From 96ca97756fc001542d123285b4f5e7a282098bc5 Mon Sep 17 00:00:00 2001 From: Mike Wall Date: Fri, 15 May 2026 09:13:10 -0700 Subject: [PATCH 01/16] Use restart forces when annealing. --- examples/gpmdk/src/gpmdcov_mdloop.F90 | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/examples/gpmdk/src/gpmdcov_mdloop.F90 b/examples/gpmdk/src/gpmdcov_mdloop.F90 index 8ef2bad..c5f4c34 100644 --- a/examples/gpmdk/src/gpmdcov_mdloop.F90 +++ b/examples/gpmdk/src/gpmdcov_mdloop.F90 @@ -36,6 +36,7 @@ subroutine gpmdcov_MDloop() real(dp) :: ke_tensor(3,3) real(dp) :: pressure_tensor(3,3) real(dp), allocatable :: saved_velocities(:,:) + real(dp), allocatable :: saved_forces(:,:) integer :: total_steps integer :: cuda_error logical :: newnl ! Indicates new neighbor list @@ -78,6 +79,7 @@ end function cudaProfilerStop !do mdstep = -1,lt%mdsteps if(gpmdt%minimization_steps.ne.0)then saved_velocities = sy%velocity + saved_forces = sy%force sy%velocity = 0.0_dp endif ! Compute box volume @@ -728,7 +730,9 @@ end function cudaProfilerStop if(mdstep.eq.gpmdt%minimization_steps)then if(gpmdt%temp0.gt.1.0E-10.or.gpmdt%restartfromdump)then sy%velocity = saved_velocities + sy%force = saved_forces deallocate(saved_velocities) + deallocate(saved_forces) endif endif From 11565e72a1fde2c3c48e5623ec03efa2df0f2d61 Mon Sep 17 00:00:00 2001 From: Mike Wall Date: Mon, 1 Jun 2026 12:13:04 -0700 Subject: [PATCH 02/16] Set initial force to zero during annealing --- examples/gpmdk/src/gpmdcov_mdloop.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/examples/gpmdk/src/gpmdcov_mdloop.F90 b/examples/gpmdk/src/gpmdcov_mdloop.F90 index c5f4c34..05142d3 100644 --- a/examples/gpmdk/src/gpmdcov_mdloop.F90 +++ b/examples/gpmdk/src/gpmdcov_mdloop.F90 @@ -81,6 +81,7 @@ end function cudaProfilerStop saved_velocities = sy%velocity saved_forces = sy%force sy%velocity = 0.0_dp + sy%force = 0.0_dp endif ! Compute box volume call gpmdcov_get_vol(sy%lattice_vector,sy%volr) From 4b86196c7c366292ae5305992cc20aada959e988 Mon Sep 17 00:00:00 2001 From: Sue Mniszewski Date: Wed, 15 Apr 2026 13:13:56 -0600 Subject: [PATCH 03/16] Added water input file for testing. --- examples/gpmdk/run/water/my_waterInput.in | 130 ++++++++++++++++++++++ 1 file changed, 130 insertions(+) create mode 100644 examples/gpmdk/run/water/my_waterInput.in diff --git a/examples/gpmdk/run/water/my_waterInput.in b/examples/gpmdk/run/water/my_waterInput.in new file mode 100644 index 0000000..a511072 --- /dev/null +++ b/examples/gpmdk/run/water/my_waterInput.in @@ -0,0 +1,130 @@ +INPUT FILE FOR THE GPMD PROGRAM +=============================== + +#LATTE parameters +Latte{ + JobName= GPMD + #BMLType= Ellpack + BMLType= Dense + #Method= GSP2 + #Method= SP2 + #Method= Diag + #Method= DiagEf + Method= DiagEfFull + MDim= -1 + #Threshold= 1.0d-5 + Threshold= 0.0 + #Verbose= 2 #Verbosity levels: Basic info(0), 1(Basic routines info), 2(Print Physics data), 3(Print Relevant Matrices), 5(Print auxiliary matrices), 10(Print all) + Verbose= 10 #Verbosity levels: Basic info(0), 1(Basic routines info), 2(Print Physics data), 3(Print Relevant Matrices), 5(Print auxiliary matrices), 10(Print all) + #SCF variables# + #StopAt= "gpmdcov_Energ" + #StopAt= "gpmdcov_DM_Min" + #StopAt= "gpmdcov_FirstCharges" + MPulay= 10 + #ZMat= ZSP + ZMat= Diag + PulayCoeff= 0.1 + #MixCoeff= 0.6 #VALID FOR WAT + MixCoeff= 0.2 + SCFTol= 1.0d-8 + MaxSCFIter= 500 + CoulAcc= 1.0d-5 + TimeRatio= 10.0 + #TimeStep= 0.2 + TimeStep= 0.4 + #TimeStep= 0.00 + MDSteps= 2000 + ParamPath= "../../tests/latteTBparams" + CoordsFile= coords_300.dat + #CoordsFile= coords_2088.dat + NlistEach= 10 + MuCalcType= FromParts + EFermi= -0.0 + #kBT= 0.04308695 + #kBT= 0.025 + kBT= 0.2 + Entropy= T + DoKernel= F +} + +#SP2 Solver +SP2{ + MinSP2Iter= 10 + MaxSP2Iter= 200 + SP2Tol= 1.0d-5 + SP2Conv= Rel +} + +#Graph-based SP2 parameters +GSP2{ + + BMLType= Ellpack + GraphElement= Atom + #PartitionType= Box + PartitionType= Block + #PartitionType= Sedacs + NLGraphCut= 4.5 + CovGraphFact= 4.5 + NodesPerPart= 300 + PartitionCount= 1 + PartitionCountX= 1 + PartitionCountY= 1 + PartitionCountZ= 1 + GraphThreshold= 0.00001 + ErrLimit= 1.0e-12 + PartEach= 1000 + SmallSubgraphs= T + Alpha= 10 + Mdim= -1 +} + + +#Sparse propagation of the inverse overlap +ZSP{ + Verbose= 1 + NFirst= 8 + NRefI= 3 + NRefF= 1 + Int= .true. + NumthreshI= 1.0d-8 + NumthreshF= 1.0d-5 +} + +#Extended Lagrangian parameters +XLBO{ + JobName= XLBO + Verbose= 1 + Mprg_init= 2 + MaxSCFIter= 0 + MaxSCFInitIter= 50 + NumThresh= 0.0 +} + + +KERNEL{ + XLBOLevel1= T + ScaledDelta= T + ScaledDeltaConstant= 0.2 + KernelType= ByParts + #KernelType= Full + #KernelType= ByBlocks + BuildAlways= F + RankNUpdate= 2 + KernelMixing= T + InitialMixingWith= DIIS + UpdateEach= 1 + UpdateAfterBuild= T + Verbose= 1 +} + +GPMD{ + DoVelocityRescale= F + #VRFactor= 1.0 + WriteTrajectory= F + WriteCoordsEach= 10 + LangevinMethod= Siva + LangevinDynamics= F + LangevinGamma= 0.01 + InitialTemperature= 300.0 + SymmetrizeGraph= T +} From da80ee0f176c8b4686bd343a1c100f35b3c4143e Mon Sep 17 00:00:00 2001 From: Sue Mniszewski Date: Wed, 15 Apr 2026 15:28:27 -0600 Subject: [PATCH 04/16] Added diff changes that did not use si. Ran for t2, t3, t34, t4. --- examples/gpmdk/src/gpmdcov_mdloop.F90 | 29 +++++++++++++++++++++------ 1 file changed, 23 insertions(+), 6 deletions(-) diff --git a/examples/gpmdk/src/gpmdcov_mdloop.F90 b/examples/gpmdk/src/gpmdcov_mdloop.F90 index 05142d3..0953062 100644 --- a/examples/gpmdk/src/gpmdcov_mdloop.F90 +++ b/examples/gpmdk/src/gpmdcov_mdloop.F90 @@ -37,6 +37,9 @@ subroutine gpmdcov_MDloop() real(dp) :: pressure_tensor(3,3) real(dp), allocatable :: saved_velocities(:,:) real(dp), allocatable :: saved_forces(:,:) + real(dp) :: user_timestep,this_maxdisp + real(dp), parameter :: maxdist = 0.02 + integer :: si,num_substeps integer :: total_steps integer :: cuda_error logical :: newnl ! Indicates new neighbor list @@ -75,7 +78,7 @@ end function cudaProfilerStop endif call gpmdcov_msI("gpmdcov_MDloop","In gpmdcov_MDloop ...",lt%verbose,myRank) - savets = lt%timestep + !savets = lt%timestep !do mdstep = -1,lt%mdsteps if(gpmdt%minimization_steps.ne.0)then saved_velocities = sy%velocity @@ -92,6 +95,8 @@ end function cudaProfilerStop call freeze(gpmdt%freezef,freeze_list,sy%velocity) endif + user_timestep = lt%timestep + do mdstep = 1,total_steps ! if(mdstep < 0)then ! savets = lt%timestep @@ -125,6 +130,12 @@ end function cudaProfilerStop write(*,*)"" endif + this_maxdisp = maxval(user_timestep*sy%velocity) + num_substeps = 1 + int(this_maxdisp/maxdist) ! If max displacement is above 0.02 Angstroms, divide into substeps + if(num_substeps > 1)write(*,*)"Performing ",num_substeps," substeps for mdstep ",mdstep + lt%timestep = user_timestep/num_substeps + !do si = 1,num_substeps + maxv_atom_axis = MAXLOC(ABS(sy%velocity)) call gpmdcov_msI("gpmdcov_MDloop","Maximum Velocity "//to_string(MAXVAL(ABS(sy%velocity)))//" & &for (atom,axis) = ("//to_string(maxv_atom_axis(2))//","//to_string(maxv_atom_axis(1))//")",lt%verbose,myRank) @@ -148,7 +159,8 @@ end function cudaProfilerStop !! Total Energy in eV Energy = EKIN + EPOT; !! Time in fs - Time = mdstep*lt%timestep; + !Time = mdstep*lt%timestep; + Time = mdstep*user_timestep; !! Statistical pressure do i = 1,3 @@ -359,7 +371,8 @@ end function cudaProfilerStop !> Update neighbor list (Actialized every nlisteach times steps) mls_md1 = mls() - if(mod(mdstep,lt%nlisteach) == 0 .or. mdstep == 0 .or. mdstep == 1)then + !if(mod(mdstep,lt%nlisteach) == 0 .or. mdstep == 0 .or. mdstep == 1)then + if((mod(mdstep,lt%nlisteach) == 0 .or. mdstep == 0 .or. mdstep == 1))then call gpmdcov_msMemGPU("mdloop","Before NeighborList",lt%verbose,myRank) call gpmdcov_msMem("gpmdcov_mdloop", "Before build_nlist_int",lt%verbose,myRank) !call gpmdcov_destroy_nlist(nl,lt%verbose) @@ -373,8 +386,9 @@ end function cudaProfilerStop #ifdef USE_OFFLOAD call gpmdcov_build_nlist_sedacs(sy%coordinate,sy%lattice_vector,coulcut,nl,lt%verbose,myRank,numRanks) #else - call gpmdcov_build_nlist_sedacs(sy%coordinate,sy%lattice_vector,coulcut,nl,lt%verbose,myRank,numRanks) + !call gpmdcov_build_nlist_sedacs(sy%coordinate,sy%lattice_vector,coulcut,nl,lt%verbose,myRank,numRanks) !call gpmdcov_build_nlist_sparse_v2(sy%coordinate,sy%lattice_vector,coulcut,nl,lt%verbose,myRank,numRanks) + call gpmdcov_build_nlist_sparse_v2(sy%coordinate,sy%lattice_vector,coulcut,nl,lt%verbose,myRank,numRanks) #endif ! if(any(nl2%nrnnstruct.ne.nl%nrnnstruct))then ! write(*,*)"DEBUG: nrnnstruct not equal" @@ -460,7 +474,8 @@ end function cudaProfilerStop mls_md1 = mls() resnorm = 0.0_dp - if((mdstep >= 2) .and. (.not. (kernel%xlbolevel1.and.lt%doKernel))) resnorm = norm2(sy%net_charge - n)/sqrt(dble(sy%nats)) + !if((mdstep >= 2) .and. (.not. (kernel%xlbolevel1.and.lt%doKernel))) resnorm = norm2(sy%net_charge - n)/sqrt(dble(sy%nats)) + if((mdstep >= 2) .and. (.not. kernel%xlbolevel1)) resnorm = norm2(sy%net_charge - n)/sqrt(dble(sy%nats)) Nr_SCF_It = xl%maxscfiter; !> Use SCF the first MD steps @@ -513,7 +528,8 @@ end function cudaProfilerStop #ifdef USE_NVTX call gpmdEndRange #endif - if(kernel%xlbolevel1.and.lt%doKernel)then + if(kernel%xlbolevel1.and.lt%doKernel)then + !if(kernel%xlbolevel1)then allocate(n1(sy%nats)) if(mdstep > 1)then !sy%net_charge = n @@ -576,6 +592,7 @@ end function cudaProfilerStop call gpmdcov_msMem("gpmdcov_mdloop", "Before gpmdcov_EnergAndForces",lt%verbose,myRank) if(kernel%xlbolevel1.and.lt%doKernel)then + !if(kernel%xlbolevel1)then if(mdstep <= 1) n1 = n call gpmdcov_EnergAndForces(n1) deallocate(n1) From 86b4c9ee14bc299333e34203e0ad494cccd545e1 Mon Sep 17 00:00:00 2001 From: Sue Mniszewski Date: Mon, 20 Apr 2026 13:11:08 -0600 Subject: [PATCH 05/16] Added lines for substeps (si) in gpmdcov_mdloop. --- examples/gpmdk/src/gpmdcov_mdloop.F90 | 36 +++++++++++++++++++++------ 1 file changed, 29 insertions(+), 7 deletions(-) diff --git a/examples/gpmdk/src/gpmdcov_mdloop.F90 b/examples/gpmdk/src/gpmdcov_mdloop.F90 index 0953062..fa65c40 100644 --- a/examples/gpmdk/src/gpmdcov_mdloop.F90 +++ b/examples/gpmdk/src/gpmdcov_mdloop.F90 @@ -134,7 +134,7 @@ end function cudaProfilerStop num_substeps = 1 + int(this_maxdisp/maxdist) ! If max displacement is above 0.02 Angstroms, divide into substeps if(num_substeps > 1)write(*,*)"Performing ",num_substeps," substeps for mdstep ",mdstep lt%timestep = user_timestep/num_substeps - !do si = 1,num_substeps + do si = 1,num_substeps maxv_atom_axis = MAXLOC(ABS(sy%velocity)) call gpmdcov_msI("gpmdcov_MDloop","Maximum Velocity "//to_string(MAXVAL(ABS(sy%velocity)))//" & @@ -172,7 +172,8 @@ end function cudaProfilerStop pressure_tensor = EVOVERV2P*(ke_tensor + virial)/sy%volr - if(myRank == 1)then + !if(myRank == 1)then + if((myRank == 1).and.(si.eq.1))then write(*,*)"Time [fs] = ",Time write(*,*)"Energy Kinetic [eV] = ",EKIN write(*,*)"Energy Potential [eV] = ",EPOT @@ -199,6 +200,9 @@ end function cudaProfilerStop write(*,*)i,sy%velocity(1,i),sy%velocity(2,i),sy%velocity(3,i) enddo endif + + if(si.eq.num_substeps)then + !> Update positions call gpmdcov_msMem("gpmdcov_mdloop", "Before updatecoords",lt%verbose,myRank) if(myRank == 1 .and. lt%verbose >= 1) call prg_timer_start(dyn_timer,"Update positions") @@ -446,13 +450,28 @@ end function cudaProfilerStop call gpmdStartRange("Part",4) #endif - call gpmdcov_Part(2) + !call gpmdcov_Part(2) + if (num_substeps.eq.1) then + call gpmdcov_Part(2) + elseif (si.eq.1) then + call gpmdcov_Part(4) + else + call gpmdcov_Part(3) + endif + ! if(si.eq.num_steps)then + ! call gpmdcov_Part(3) + ! else + ! call gpmdcov_Part(2) + ! endif + #ifdef USE_NVTX call gpmdEndRange #endif call gpmdcov_msMem("gpmdcov_mdloop", "After gpmdcov_Part",lt%verbose,myRank) call gpmdcov_msI("gpmdcov_MDloop","Time for gpmdcov_Part & &"//to_string(mls() - mls_i)//" ms",lt%verbose,myRank) + + endif ! if (si.eq.substeps) !> Reprg_initialize parts. mls_i = mls() call gpmdcov_msMem("gpmdcov_mdloop", "Before gpmdcov_InitParts",lt%verbose,myRank) @@ -570,7 +589,9 @@ end function cudaProfilerStop mls_md1 = mls() call gpmdcov_msI("gpmdcov_MDloop","ResNorm = "//to_string(resnorm),lt%verbose,myRank) - if(myRank == 1)then + + !if(myRank == 1)then + if(myRank == 1.and.si.eq.1)then if(mdstep.le.gpmdt%minimization_steps)then if(.not.gpmdt%anneal_graph)then write(*,'(A35,I15,A1,F18.5,A1,ES12.5,A1,ES12.5,A1,ES12.5)')"Minstep, Energy, Egap, Resnorm, Temp", & @@ -580,10 +601,9 @@ end function cudaProfilerStop &mdstep," ", Energy," ", egap_glob," ", resnorm," ", Temp endif else - write(*,'(A35,I15,A1,F18.5,A1,ES12.5,A1,ES12.5,A1,ES12.5)')"Mdstep, Energy, Egap, Resnorm, Temp", & - &mdstep-gpmdt%minimization_steps," ", Energy," ", egap_glob," ", resnorm," ", Temp + write(*,'(A35,I15,A1,I15,A1,F18.5,A1,ES12.5,A1,ES12.5,A1,ES12.5)')"Mdstep, Substeps, Energy, Egap, Resnorm, Temp", & + &mdstep-gpmdt%minimization_steps," ", num_substeps, " ", Energy," ", egap_glob," ", resnorm," ", Temp endif - !write(*,*)"Step, Energy, EGap, Resnorm", mdstep, Energy, egap_glob, resnorm endif #ifdef USE_NVTX call gpmdStartRange("EnergAndForces",7) @@ -713,6 +733,8 @@ end function cudaProfilerStop sy%velocity = 0.0_dp endif + enddo ! end of si loop + #ifdef USE_NVTX call gpmdStartRange("Write trajectory",3) #endif From cb25db2b9406861e803621bb6ec75aea8cb8ca34 Mon Sep 17 00:00:00 2001 From: Sue Mniszewski Date: Wed, 22 Apr 2026 12:37:06 -0600 Subject: [PATCH 06/16] Added "si" lines for testing. --- examples/gpmdk/run/water/my_waterInput.in | 4 ++-- examples/gpmdk/src/gpmdcov_mdloop.F90 | 19 ++++++++++++++++--- 2 files changed, 18 insertions(+), 5 deletions(-) diff --git a/examples/gpmdk/run/water/my_waterInput.in b/examples/gpmdk/run/water/my_waterInput.in index a511072..3c6b67c 100644 --- a/examples/gpmdk/run/water/my_waterInput.in +++ b/examples/gpmdk/run/water/my_waterInput.in @@ -14,8 +14,8 @@ Latte{ MDim= -1 #Threshold= 1.0d-5 Threshold= 0.0 - #Verbose= 2 #Verbosity levels: Basic info(0), 1(Basic routines info), 2(Print Physics data), 3(Print Relevant Matrices), 5(Print auxiliary matrices), 10(Print all) - Verbose= 10 #Verbosity levels: Basic info(0), 1(Basic routines info), 2(Print Physics data), 3(Print Relevant Matrices), 5(Print auxiliary matrices), 10(Print all) + Verbose= 2 #Verbosity levels: Basic info(0), 1(Basic routines info), 2(Print Physics data), 3(Print Relevant Matrices), 5(Print auxiliary matrices), 10(Print all) + #Verbose= 10 #Verbosity levels: Basic info(0), 1(Basic routines info), 2(Print Physics data), 3(Print Relevant Matrices), 5(Print auxiliary matrices), 10(Print all) #SCF variables# #StopAt= "gpmdcov_Energ" #StopAt= "gpmdcov_DM_Min" diff --git a/examples/gpmdk/src/gpmdcov_mdloop.F90 b/examples/gpmdk/src/gpmdcov_mdloop.F90 index fa65c40..61c645d 100644 --- a/examples/gpmdk/src/gpmdcov_mdloop.F90 +++ b/examples/gpmdk/src/gpmdcov_mdloop.F90 @@ -134,6 +134,9 @@ end function cudaProfilerStop num_substeps = 1 + int(this_maxdisp/maxdist) ! If max displacement is above 0.02 Angstroms, divide into substeps if(num_substeps > 1)write(*,*)"Performing ",num_substeps," substeps for mdstep ",mdstep lt%timestep = user_timestep/num_substeps + write(*,*)"timestep = ", lt%timestep + + ! Substeps loop do si = 1,num_substeps maxv_atom_axis = MAXLOC(ABS(sy%velocity)) @@ -173,6 +176,7 @@ end function cudaProfilerStop pressure_tensor = EVOVERV2P*(ke_tensor + virial)/sy%volr !if(myRank == 1)then + ! Change? if((myRank == 1).and.(si.eq.1))then write(*,*)"Time [fs] = ",Time write(*,*)"Energy Kinetic [eV] = ",EKIN @@ -201,6 +205,7 @@ end function cudaProfilerStop enddo endif + ! This should be OK if(si.eq.num_substeps)then !> Update positions @@ -450,6 +455,7 @@ end function cudaProfilerStop call gpmdStartRange("Part",4) #endif + ! Check what si should be !call gpmdcov_Part(2) if (num_substeps.eq.1) then call gpmdcov_Part(2) @@ -458,7 +464,8 @@ end function cudaProfilerStop else call gpmdcov_Part(3) endif - ! if(si.eq.num_steps)then + + ! if(si.eq.num_steps)then ! call gpmdcov_Part(3) ! else ! call gpmdcov_Part(2) @@ -471,7 +478,9 @@ end function cudaProfilerStop call gpmdcov_msI("gpmdcov_MDloop","Time for gpmdcov_Part & &"//to_string(mls() - mls_i)//" ms",lt%verbose,myRank) - endif ! if (si.eq.substeps) + ! Should be OK + endif ! if (si.eq.num_substeps) + !> Reprg_initialize parts. mls_i = mls() call gpmdcov_msMem("gpmdcov_mdloop", "Before gpmdcov_InitParts",lt%verbose,myRank) @@ -590,8 +599,12 @@ end function cudaProfilerStop mls_md1 = mls() call gpmdcov_msI("gpmdcov_MDloop","ResNorm = "//to_string(resnorm),lt%verbose,myRank) + ! Which is right? !if(myRank == 1)then - if(myRank == 1.and.si.eq.1)then + !if(si.eq.num_substeps) + ! Check num_substeps == 1 or (num_substeps== 2 and si = 2)? + !if(myRank == 1.and.si.eq.1)then + if(myRank == 1.and.si.eq.num_substeps)then if(mdstep.le.gpmdt%minimization_steps)then if(.not.gpmdt%anneal_graph)then write(*,'(A35,I15,A1,F18.5,A1,ES12.5,A1,ES12.5,A1,ES12.5)')"Minstep, Energy, Egap, Resnorm, Temp", & From a54b157bcc6b1a867a7288c049e9714fabb1ecd6 Mon Sep 17 00:00:00 2001 From: Sue Mniszewski Date: Thu, 23 Apr 2026 14:52:47 -0600 Subject: [PATCH 07/16] Commented out "if" and "endif" inside of si loop. --- examples/gpmdk/run/water/my_waterInput.in | 2 +- examples/gpmdk/src/gpmdcov_mdloop.F90 | 11 ++++++++--- 2 files changed, 9 insertions(+), 4 deletions(-) diff --git a/examples/gpmdk/run/water/my_waterInput.in b/examples/gpmdk/run/water/my_waterInput.in index 3c6b67c..c89034b 100644 --- a/examples/gpmdk/run/water/my_waterInput.in +++ b/examples/gpmdk/run/water/my_waterInput.in @@ -31,7 +31,7 @@ Latte{ CoulAcc= 1.0d-5 TimeRatio= 10.0 #TimeStep= 0.2 - TimeStep= 0.4 + TimeStep= 0.2 #TimeStep= 0.00 MDSteps= 2000 ParamPath= "../../tests/latteTBparams" diff --git a/examples/gpmdk/src/gpmdcov_mdloop.F90 b/examples/gpmdk/src/gpmdcov_mdloop.F90 index 61c645d..a946980 100644 --- a/examples/gpmdk/src/gpmdcov_mdloop.F90 +++ b/examples/gpmdk/src/gpmdcov_mdloop.F90 @@ -132,6 +132,7 @@ end function cudaProfilerStop this_maxdisp = maxval(user_timestep*sy%velocity) num_substeps = 1 + int(this_maxdisp/maxdist) ! If max displacement is above 0.02 Angstroms, divide into substeps + if (num_substeps > 2) num_substeps=2 if(num_substeps > 1)write(*,*)"Performing ",num_substeps," substeps for mdstep ",mdstep lt%timestep = user_timestep/num_substeps write(*,*)"timestep = ", lt%timestep @@ -186,6 +187,8 @@ end function cudaProfilerStop write(*,*)"Pressure [bar] = ",pressure_tensor(1,1)+pressure_tensor(2,2)+pressure_tensor(3,3) endif + !if(si.eq.num_substeps)then + call gpmdcov_msI("gpmdcov_MDloop","Time for Preliminaries "//to_string(mls() - mls_md1)//" ms",lt%verbose,myRank) mls_md1 = mls() @@ -205,8 +208,7 @@ end function cudaProfilerStop enddo endif - ! This should be OK - if(si.eq.num_substeps)then + !if(si.eq.num_substeps)then !> Update positions call gpmdcov_msMem("gpmdcov_mdloop", "Before updatecoords",lt%verbose,myRank) @@ -378,6 +380,9 @@ end function cudaProfilerStop endif call gpmdcov_msI("gpmdcov_MDloop","Time for prg_xlbo_nint "//to_string(mls() - mls_md1)//" ms",lt%verbose,myRank) + ! causes problem + !if(si.eq.num_substeps)then + !> Update neighbor list (Actialized every nlisteach times steps) mls_md1 = mls() !if(mod(mdstep,lt%nlisteach) == 0 .or. mdstep == 0 .or. mdstep == 1)then @@ -479,7 +484,7 @@ end function cudaProfilerStop &"//to_string(mls() - mls_i)//" ms",lt%verbose,myRank) ! Should be OK - endif ! if (si.eq.num_substeps) + !endif ! if (si.eq.num_substeps) !> Reprg_initialize parts. mls_i = mls() From 87f3e1d7d5ee0479287608991cd94c69a8be2d75 Mon Sep 17 00:00:00 2001 From: Sue Mniszewski Date: Wed, 29 Apr 2026 10:02:51 -0600 Subject: [PATCH 08/16] REmoved si code. Added new logic to top of mdstep loop. --- examples/gpmdk/src/gpmdcov_mdloop.F90 | 81 +++++++++++---------------- 1 file changed, 33 insertions(+), 48 deletions(-) diff --git a/examples/gpmdk/src/gpmdcov_mdloop.F90 b/examples/gpmdk/src/gpmdcov_mdloop.F90 index a946980..2aceb0d 100644 --- a/examples/gpmdk/src/gpmdcov_mdloop.F90 +++ b/examples/gpmdk/src/gpmdcov_mdloop.F90 @@ -39,7 +39,7 @@ subroutine gpmdcov_MDloop() real(dp), allocatable :: saved_forces(:,:) real(dp) :: user_timestep,this_maxdisp real(dp), parameter :: maxdist = 0.02 - integer :: si,num_substeps + logical :: first_substep_taken integer :: total_steps integer :: cuda_error logical :: newnl ! Indicates new neighbor list @@ -96,6 +96,7 @@ end function cudaProfilerStop endif user_timestep = lt%timestep + first_substep_taken = .false. do mdstep = 1,total_steps ! if(mdstep < 0)then @@ -130,15 +131,34 @@ end function cudaProfilerStop write(*,*)"" endif + if (first_substep_taken .eqv. .true.) then + write(*,*) "for mdstep ", mdstep, "first_substep_taken is TRUE" + else + write(*,*) "for mdstep ", mdstep, "first_substep_taken is FALSE" + endif this_maxdisp = maxval(user_timestep*sy%velocity) - num_substeps = 1 + int(this_maxdisp/maxdist) ! If max displacement is above 0.02 Angstroms, divide into substeps - if (num_substeps > 2) num_substeps=2 - if(num_substeps > 1)write(*,*)"Performing ",num_substeps," substeps for mdstep ",mdstep - lt%timestep = user_timestep/num_substeps - write(*,*)"timestep = ", lt%timestep + write(*,*)"for mdstep ", mdstep, "this_maxdisp = ", this_maxdisp + if (first_substep_taken .or.(this_maxdisp > maxdist)) then + write(*,*)"Splitting mdstep ", mdstep + lt%timestep = user_timestep/2.0 + write(*,*)"for mdstep ", mdstep, "reduced timestep = ", lt%timestep + + if (first_substep_taken) then + first_substep_taken = .false. + else + first_substep_taken = .true. + endif - ! Substeps loop - do si = 1,num_substeps + else + lt%timestep = user_timestep + endif + + !Performing ",num_substeps," substeps for mdstep ",mdstep + !num_substeps = 1 + int(this_maxdisp/maxdist) ! If max displacement is above 0.02 Angstroms, divide into substeps + !if (num_substeps > 2) num_substeps=2 + !if(num_substeps > 1)write(*,*)"Performing ",num_substeps," substeps for mdstep ",mdstep + !lt%timestep = user_timestep/num_substeps + !write(*,*)"timestep = ", lt%timestep maxv_atom_axis = MAXLOC(ABS(sy%velocity)) call gpmdcov_msI("gpmdcov_MDloop","Maximum Velocity "//to_string(MAXVAL(ABS(sy%velocity)))//" & @@ -176,9 +196,7 @@ end function cudaProfilerStop pressure_tensor = EVOVERV2P*(ke_tensor + virial)/sy%volr - !if(myRank == 1)then - ! Change? - if((myRank == 1).and.(si.eq.1))then + if(myRank == 1)then write(*,*)"Time [fs] = ",Time write(*,*)"Energy Kinetic [eV] = ",EKIN write(*,*)"Energy Potential [eV] = ",EPOT @@ -187,8 +205,6 @@ end function cudaProfilerStop write(*,*)"Pressure [bar] = ",pressure_tensor(1,1)+pressure_tensor(2,2)+pressure_tensor(3,3) endif - !if(si.eq.num_substeps)then - call gpmdcov_msI("gpmdcov_MDloop","Time for Preliminaries "//to_string(mls() - mls_md1)//" ms",lt%verbose,myRank) mls_md1 = mls() @@ -208,8 +224,6 @@ end function cudaProfilerStop enddo endif - !if(si.eq.num_substeps)then - !> Update positions call gpmdcov_msMem("gpmdcov_mdloop", "Before updatecoords",lt%verbose,myRank) if(myRank == 1 .and. lt%verbose >= 1) call prg_timer_start(dyn_timer,"Update positions") @@ -380,9 +394,6 @@ end function cudaProfilerStop endif call gpmdcov_msI("gpmdcov_MDloop","Time for prg_xlbo_nint "//to_string(mls() - mls_md1)//" ms",lt%verbose,myRank) - ! causes problem - !if(si.eq.num_substeps)then - !> Update neighbor list (Actialized every nlisteach times steps) mls_md1 = mls() !if(mod(mdstep,lt%nlisteach) == 0 .or. mdstep == 0 .or. mdstep == 1)then @@ -460,21 +471,7 @@ end function cudaProfilerStop call gpmdStartRange("Part",4) #endif - ! Check what si should be - !call gpmdcov_Part(2) - if (num_substeps.eq.1) then - call gpmdcov_Part(2) - elseif (si.eq.1) then - call gpmdcov_Part(4) - else - call gpmdcov_Part(3) - endif - - ! if(si.eq.num_steps)then - ! call gpmdcov_Part(3) - ! else - ! call gpmdcov_Part(2) - ! endif + call gpmdcov_Part(2) #ifdef USE_NVTX call gpmdEndRange @@ -483,9 +480,6 @@ end function cudaProfilerStop call gpmdcov_msI("gpmdcov_MDloop","Time for gpmdcov_Part & &"//to_string(mls() - mls_i)//" ms",lt%verbose,myRank) - ! Should be OK - !endif ! if (si.eq.num_substeps) - !> Reprg_initialize parts. mls_i = mls() call gpmdcov_msMem("gpmdcov_mdloop", "Before gpmdcov_InitParts",lt%verbose,myRank) @@ -562,7 +556,6 @@ end function cudaProfilerStop call gpmdEndRange #endif if(kernel%xlbolevel1.and.lt%doKernel)then - !if(kernel%xlbolevel1)then allocate(n1(sy%nats)) if(mdstep > 1)then !sy%net_charge = n @@ -604,12 +597,7 @@ end function cudaProfilerStop mls_md1 = mls() call gpmdcov_msI("gpmdcov_MDloop","ResNorm = "//to_string(resnorm),lt%verbose,myRank) - ! Which is right? - !if(myRank == 1)then - !if(si.eq.num_substeps) - ! Check num_substeps == 1 or (num_substeps== 2 and si = 2)? - !if(myRank == 1.and.si.eq.1)then - if(myRank == 1.and.si.eq.num_substeps)then + if(myRank == 1)then if(mdstep.le.gpmdt%minimization_steps)then if(.not.gpmdt%anneal_graph)then write(*,'(A35,I15,A1,F18.5,A1,ES12.5,A1,ES12.5,A1,ES12.5)')"Minstep, Energy, Egap, Resnorm, Temp", & @@ -619,8 +607,8 @@ end function cudaProfilerStop &mdstep," ", Energy," ", egap_glob," ", resnorm," ", Temp endif else - write(*,'(A35,I15,A1,I15,A1,F18.5,A1,ES12.5,A1,ES12.5,A1,ES12.5)')"Mdstep, Substeps, Energy, Egap, Resnorm, Temp", & - &mdstep-gpmdt%minimization_steps," ", num_substeps, " ", Energy," ", egap_glob," ", resnorm," ", Temp + write(*,'(A35,I15,A1,F5.2,A1,F18.5,A1,ES12.5,A1,ES12.5,A1,ES12.5)')"Mdstep, timestep, Energy, Egap, Resnorm, Temp", & + &mdstep-gpmdt%minimization_steps," ", lt%timestep, " ", Energy," ", egap_glob," ", resnorm," ", Temp endif endif #ifdef USE_NVTX @@ -630,7 +618,6 @@ end function cudaProfilerStop call gpmdcov_msMem("gpmdcov_mdloop", "Before gpmdcov_EnergAndForces",lt%verbose,myRank) if(kernel%xlbolevel1.and.lt%doKernel)then - !if(kernel%xlbolevel1)then if(mdstep <= 1) n1 = n call gpmdcov_EnergAndForces(n1) deallocate(n1) @@ -751,8 +738,6 @@ end function cudaProfilerStop sy%velocity = 0.0_dp endif - enddo ! end of si loop - #ifdef USE_NVTX call gpmdStartRange("Write trajectory",3) #endif From 8f085fe9be36ff34b4dd26dca5d386f6cfbebad9 Mon Sep 17 00:00:00 2001 From: Sue Mniszewski Date: Mon, 4 May 2026 09:05:14 -0600 Subject: [PATCH 09/16] Added new print output for mdsteps. --- examples/gpmdk/src/gpmdcov_mdloop.F90 | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/examples/gpmdk/src/gpmdcov_mdloop.F90 b/examples/gpmdk/src/gpmdcov_mdloop.F90 index 2aceb0d..573b20c 100644 --- a/examples/gpmdk/src/gpmdcov_mdloop.F90 +++ b/examples/gpmdk/src/gpmdcov_mdloop.F90 @@ -40,7 +40,7 @@ subroutine gpmdcov_MDloop() real(dp) :: user_timestep,this_maxdisp real(dp), parameter :: maxdist = 0.02 logical :: first_substep_taken - integer :: total_steps + integer :: total_steps, print_mdstep integer :: cuda_error logical :: newnl ! Indicates new neighbor list type(neighlist_type) :: nl2 @@ -97,6 +97,7 @@ end function cudaProfilerStop user_timestep = lt%timestep first_substep_taken = .false. + print_mdstep = 0 do mdstep = 1,total_steps ! if(mdstep < 0)then @@ -607,8 +608,13 @@ end function cudaProfilerStop &mdstep," ", Energy," ", egap_glob," ", resnorm," ", Temp endif else - write(*,'(A35,I15,A1,F5.2,A1,F18.5,A1,ES12.5,A1,ES12.5,A1,ES12.5)')"Mdstep, timestep, Energy, Egap, Resnorm, Temp", & - &mdstep-gpmdt%minimization_steps," ", lt%timestep, " ", Energy," ", egap_glob," ", resnorm," ", Temp + ! Skip first substeps when writing output + if (.not.first_substep_taken)then + print_mdstep = print_mdstep + 1 + write(*,'(A35,I15,A1,F5.2,A1,F18.5,A1,ES12.5,A1,ES12.5,A1,ES12.5)')"Mdstep, timestep, Energy, Egap, Resnorm, Temp", & + &print_mdstep," ", lt%timestep, " ", Energy," ", egap_glob," ", resnorm," ", Temp + !&mdstep-gpmdt%minimization_steps," ", lt%timestep, " ", Energy," ", egap_glob," ", resnorm," ", Temp + endif endif endif #ifdef USE_NVTX From fd8dd31852af18551b41fa3da32e250d31e83493 Mon Sep 17 00:00:00 2001 From: Sue Mniszewski Date: Tue, 5 May 2026 15:03:19 -0600 Subject: [PATCH 10/16] More print changes. --- examples/gpmdk/src/gpmdcov_mdloop.F90 | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/examples/gpmdk/src/gpmdcov_mdloop.F90 b/examples/gpmdk/src/gpmdcov_mdloop.F90 index 573b20c..9272c40 100644 --- a/examples/gpmdk/src/gpmdcov_mdloop.F90 +++ b/examples/gpmdk/src/gpmdcov_mdloop.F90 @@ -611,9 +611,10 @@ end function cudaProfilerStop ! Skip first substeps when writing output if (.not.first_substep_taken)then print_mdstep = print_mdstep + 1 - write(*,'(A35,I15,A1,F5.2,A1,F18.5,A1,ES12.5,A1,ES12.5,A1,ES12.5)')"Mdstep, timestep, Energy, Egap, Resnorm, Temp", & - &print_mdstep," ", lt%timestep, " ", Energy," ", egap_glob," ", resnorm," ", Temp - !&mdstep-gpmdt%minimization_steps," ", lt%timestep, " ", Energy," ", egap_glob," ", resnorm," ", Temp + write(*,'(A35,I15,A1,F18.5,A1,ES12.5,A1,ES12.5,A1,ES12.5)')"Mdstep, Energy, Egap, Resnorm, Temp", & + &print_mdstep," ", Energy," ", egap_glob," ", resnorm," ", Temp + !&mdstep-gpmdt%minimization_steps," ", Energy," ", egap_glob," ", resnorm," ", Temp + write(*,*) "Mdstep ", print_mdstep, " was performed using two half timesteps" endif endif endif From be1897d86fce7c8a128bc793d56f6125c60177b6 Mon Sep 17 00:00:00 2001 From: Mike Wall Date: Tue, 5 May 2026 13:28:12 -0600 Subject: [PATCH 11/16] Make adaptive time step compatible with minimization --- examples/gpmdk/src/gpmdcov_mdloop.F90 | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/examples/gpmdk/src/gpmdcov_mdloop.F90 b/examples/gpmdk/src/gpmdcov_mdloop.F90 index 9272c40..49abfbd 100644 --- a/examples/gpmdk/src/gpmdcov_mdloop.F90 +++ b/examples/gpmdk/src/gpmdcov_mdloop.F90 @@ -611,10 +611,8 @@ end function cudaProfilerStop ! Skip first substeps when writing output if (.not.first_substep_taken)then print_mdstep = print_mdstep + 1 - write(*,'(A35,I15,A1,F18.5,A1,ES12.5,A1,ES12.5,A1,ES12.5)')"Mdstep, Energy, Egap, Resnorm, Temp", & - &print_mdstep," ", Energy," ", egap_glob," ", resnorm," ", Temp - !&mdstep-gpmdt%minimization_steps," ", Energy," ", egap_glob," ", resnorm," ", Temp - write(*,*) "Mdstep ", print_mdstep, " was performed using two half timesteps" + write(*,'(A35,I15,A1,F5.2,A1,F18.5,A1,ES12.5,A1,ES12.5,A1,ES12.5)')"Mdstep, Energy, Egap, Resnorm, Temp", & + &print_mdstep-gpmdt%minimization_steps," ", lt%timestep, " ", Energy," ", egap_glob," ", resnorm," ", Temp endif endif endif @@ -748,15 +746,15 @@ end function cudaProfilerStop #ifdef USE_NVTX call gpmdStartRange("Write trajectory",3) #endif - if(gpmdt%writetraj .and. myRank == 1 .and. mdstep.ge.gpmdt%minimization_steps)then + if(gpmdt%writetraj .and. myRank == 1 .and. mdstep.ge.gpmdt%minimization_steps .and. first_substep_taken .eqv. .false.)then if((gpmdt%traj_format .eq. "XYZ").and. & - (mod(mdstep-gpmdt%minimization_steps,gpmdt%writetreach).eq.0.or. & - (mdstep-gpmdt%minimization_steps).eq.1))then - call prg_write_trajectory(sy,mdstep-gpmdt%minimization_steps,gpmdt%writetreach,& + (mod(print_mdstep-gpmdt%minimization_steps,gpmdt%writetreach).eq.0.or. & + (print_mdstep-gpmdt%minimization_steps).eq.1))then + call prg_write_trajectory(sy,print_mdstep-gpmdt%minimization_steps,gpmdt%writetreach,& <%timestep,adjustl(trim(lt%jobname))//"_trajectory","xyz") call prg_write_system(sy,adjustl(trim(lt%jobname))//"_latest","pdb") else - call prg_write_trajectory(sy,mdstep-gpmdt%minimization_steps,gpmdt%writetreach,& + call prg_write_trajectory(sy,print_mdstep-gpmdt%minimization_steps,gpmdt%writetreach,& <%timestep,adjustl(trim(lt%jobname))//"_trajectory","pdb") endif endif @@ -774,7 +772,7 @@ end function cudaProfilerStop ! Save MD state each 120 steps if(gpmdt%dumpeach .gt. 0)then - if(mod(mdstep-gpmdt%minimization_steps,gpmdt%dumpeach) == 0)call gpmdcov_dump() + if(mod(print_mdstep-gpmdt%minimization_steps,gpmdt%dumpeach) == 0)call gpmdcov_dump() endif if(mdstep.eq.gpmdt%minimization_steps)then From 0dc37ba78535113cfb478e4f93e0ec022099f539 Mon Sep 17 00:00:00 2001 From: Sue Mniszewski Date: Wed, 6 May 2026 14:56:11 -0600 Subject: [PATCH 12/16] Added message to output when 2 half steps are taken. Added input file for usinf kernel and LangeDynamics. --- examples/gpmdk/run/water/my_waterInput_mod.in | 130 ++++++++++++++++++ examples/gpmdk/src/gpmdcov_mdloop.F90 | 24 +++- 2 files changed, 149 insertions(+), 5 deletions(-) create mode 100644 examples/gpmdk/run/water/my_waterInput_mod.in diff --git a/examples/gpmdk/run/water/my_waterInput_mod.in b/examples/gpmdk/run/water/my_waterInput_mod.in new file mode 100644 index 0000000..2ea8ad0 --- /dev/null +++ b/examples/gpmdk/run/water/my_waterInput_mod.in @@ -0,0 +1,130 @@ +INPUT FILE FOR THE GPMD PROGRAM +=============================== + +#LATTE parameters +Latte{ + JobName= GPMD + #BMLType= Ellpack + BMLType= Dense + #Method= GSP2 + #Method= SP2 + #Method= Diag + #Method= DiagEf + Method= DiagEfFull + MDim= -1 + #Threshold= 1.0d-5 + Threshold= 0.0 + Verbose= 2 #Verbosity levels: Basic info(0), 1(Basic routines info), 2(Print Physics data), 3(Print Relevant Matrices), 5(Print auxiliary matrices), 10(Print all) + #Verbose= 10 #Verbosity levels: Basic info(0), 1(Basic routines info), 2(Print Physics data), 3(Print Relevant Matrices), 5(Print auxiliary matrices), 10(Print all) + #SCF variables# + #StopAt= "gpmdcov_Energ" + #StopAt= "gpmdcov_DM_Min" + #StopAt= "gpmdcov_FirstCharges" + MPulay= 10 + #ZMat= ZSP + ZMat= Diag + PulayCoeff= 0.1 + #MixCoeff= 0.6 #VALID FOR WAT + MixCoeff= 0.2 + SCFTol= 1.0d-8 + MaxSCFIter= 500 + CoulAcc= 1.0d-5 + TimeRatio= 10.0 + #TimeStep= 0.2 + TimeStep= 0.4 + #TimeStep= 0.00 + MDSteps= 2000 + ParamPath= "../../tests/latteTBparams" + CoordsFile= coords_300.dat + #CoordsFile= coords_2088.dat + NlistEach= 10 + MuCalcType= FromParts + EFermi= -0.0 + #kBT= 0.04308695 + #kBT= 0.025 + kBT= 0.2 + Entropy= T + DoKernel= T +} + +#SP2 Solver +SP2{ + MinSP2Iter= 10 + MaxSP2Iter= 200 + SP2Tol= 1.0d-5 + SP2Conv= Rel +} + +#Graph-based SP2 parameters +GSP2{ + + BMLType= Ellpack + GraphElement= Atom + #PartitionType= Box + PartitionType= Block + #PartitionType= Sedacs + NLGraphCut= 4.5 + CovGraphFact= 4.5 + NodesPerPart= 300 + PartitionCount= 1 + PartitionCountX= 1 + PartitionCountY= 1 + PartitionCountZ= 1 + GraphThreshold= 0.00001 + ErrLimit= 1.0e-12 + PartEach= 1000 + SmallSubgraphs= T + Alpha= 10 + Mdim= -1 +} + + +#Sparse propagation of the inverse overlap +ZSP{ + Verbose= 1 + NFirst= 8 + NRefI= 3 + NRefF= 1 + Int= .true. + NumthreshI= 1.0d-8 + NumthreshF= 1.0d-5 +} + +#Extended Lagrangian parameters +XLBO{ + JobName= XLBO + Verbose= 1 + Mprg_init= 2 + MaxSCFIter= 0 + MaxSCFInitIter= 50 + NumThresh= 0.0 +} + + +KERNEL{ + XLBOLevel1= T + ScaledDelta= T + ScaledDeltaConstant= 0.2 + KernelType= ByParts + #KernelType= Full + #KernelType= ByBlocks + BuildAlways= F + RankNUpdate= 2 + KernelMixing= T + InitialMixingWith= DIIS + UpdateEach= 1 + UpdateAfterBuild= T + Verbose= 1 +} + +GPMD{ + DoVelocityRescale= F + #VRFactor= 1.0 + WriteTrajectory= F + WriteCoordsEach= 10 + LangevinMethod= Siva + LangevinDynamics= T + LangevinGamma= 0.01 + InitialTemperature= 300.0 + SymmetrizeGraph= T +} diff --git a/examples/gpmdk/src/gpmdcov_mdloop.F90 b/examples/gpmdk/src/gpmdcov_mdloop.F90 index 49abfbd..0985ac9 100644 --- a/examples/gpmdk/src/gpmdcov_mdloop.F90 +++ b/examples/gpmdk/src/gpmdcov_mdloop.F90 @@ -37,9 +37,9 @@ subroutine gpmdcov_MDloop() real(dp) :: pressure_tensor(3,3) real(dp), allocatable :: saved_velocities(:,:) real(dp), allocatable :: saved_forces(:,:) - real(dp) :: user_timestep,this_maxdisp + real(dp) :: user_timestep,this_maxdisp,user_half_timestep real(dp), parameter :: maxdist = 0.02 - logical :: first_substep_taken + logical :: first_substep_taken,half_timestep_flag integer :: total_steps, print_mdstep integer :: cuda_error logical :: newnl ! Indicates new neighbor list @@ -95,8 +95,17 @@ end function cudaProfilerStop call freeze(gpmdt%freezef,freeze_list,sy%velocity) endif + ! user_timestep is a full timestep + ! user_half_timestep is a half timestep + ! first_substep_taken indicates that the first of 2 half timesteps was taken + ! half_timestep_flag indicates that 2 half timesteps were used + ! an output message is printed after the mdsteps line + ! print_mdstep is the mdstep used for output + ! user_timestep = lt%timestep + user_half_timestep = lt%timestep/2.0 first_substep_taken = .false. + half_timestep_flag = .false. print_mdstep = 0 do mdstep = 1,total_steps @@ -141,7 +150,8 @@ end function cudaProfilerStop write(*,*)"for mdstep ", mdstep, "this_maxdisp = ", this_maxdisp if (first_substep_taken .or.(this_maxdisp > maxdist)) then write(*,*)"Splitting mdstep ", mdstep - lt%timestep = user_timestep/2.0 + lt%timestep = user_half_timestep + half_timestep_flag = .true. write(*,*)"for mdstep ", mdstep, "reduced timestep = ", lt%timestep if (first_substep_taken) then @@ -152,6 +162,7 @@ end function cudaProfilerStop else lt%timestep = user_timestep + half_timestep_flag = .false. endif !Performing ",num_substeps," substeps for mdstep ",mdstep @@ -611,8 +622,11 @@ end function cudaProfilerStop ! Skip first substeps when writing output if (.not.first_substep_taken)then print_mdstep = print_mdstep + 1 - write(*,'(A35,I15,A1,F5.2,A1,F18.5,A1,ES12.5,A1,ES12.5,A1,ES12.5)')"Mdstep, Energy, Egap, Resnorm, Temp", & - &print_mdstep-gpmdt%minimization_steps," ", lt%timestep, " ", Energy," ", egap_glob," ", resnorm," ", Temp + write(*,'(A35,I15,A1,F18.5,A1,ES12.5,A1,ES12.5,A1,ES12.5)')"Mdstep, Energy, Egap, Resnorm, Temp", & + &print_mdstep-gpmdt%minimization_steps," ", Energy," ", egap_glob," ", resnorm," ", Temp + if (half_timestep_flag)then + write(*,*) "Mdstep ", print_mdstep, " was performed using two half timesteps" + endif endif endif endif From 077356179650d477d0090dae99f92ed1d5d6df85 Mon Sep 17 00:00:00 2001 From: Sue Mniszewski Date: Tue, 19 May 2026 09:09:17 -0600 Subject: [PATCH 13/16] Changed message for 2 half steps.Checking criteria for mdsteps greater than minimization steps. --- examples/gpmdk/src/gpmdcov_mdloop.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/gpmdk/src/gpmdcov_mdloop.F90 b/examples/gpmdk/src/gpmdcov_mdloop.F90 index 0985ac9..405b796 100644 --- a/examples/gpmdk/src/gpmdcov_mdloop.F90 +++ b/examples/gpmdk/src/gpmdcov_mdloop.F90 @@ -148,7 +148,7 @@ end function cudaProfilerStop endif this_maxdisp = maxval(user_timestep*sy%velocity) write(*,*)"for mdstep ", mdstep, "this_maxdisp = ", this_maxdisp - if (first_substep_taken .or.(this_maxdisp > maxdist)) then + if ((first_substep_taken .or.(this_maxdisp > maxdist)) .and. mdstep.gt.gpmdt%minimization_steps) then write(*,*)"Splitting mdstep ", mdstep lt%timestep = user_half_timestep half_timestep_flag = .true. @@ -625,7 +625,7 @@ end function cudaProfilerStop write(*,'(A35,I15,A1,F18.5,A1,ES12.5,A1,ES12.5,A1,ES12.5)')"Mdstep, Energy, Egap, Resnorm, Temp", & &print_mdstep-gpmdt%minimization_steps," ", Energy," ", egap_glob," ", resnorm," ", Temp if (half_timestep_flag)then - write(*,*) "Mdstep ", print_mdstep, " was performed using two half timesteps" + write(*,*) "WARNING: Two half timesteps were performed for step ", print_mdstep endif endif endif From 096cf95266f742a1e978d94a922a3569c834113f Mon Sep 17 00:00:00 2001 From: Mike Wall Date: Thu, 21 May 2026 14:10:05 -0600 Subject: [PATCH 14/16] Use sedacs neighborlist --- examples/gpmdk/src/gpmdcov_mdloop.F90 | 8 +------- 1 file changed, 1 insertion(+), 7 deletions(-) diff --git a/examples/gpmdk/src/gpmdcov_mdloop.F90 b/examples/gpmdk/src/gpmdcov_mdloop.F90 index 405b796..885dd7e 100644 --- a/examples/gpmdk/src/gpmdcov_mdloop.F90 +++ b/examples/gpmdk/src/gpmdcov_mdloop.F90 @@ -419,14 +419,8 @@ end function cudaProfilerStop #ifdef USE_NVTX call gpmdStartRange("build_nlist_sparse_sedacs",3) #endif - !call gpmdcov_destroy_nlist(nl2,lt%verbose) -#ifdef USE_OFFLOAD call gpmdcov_build_nlist_sedacs(sy%coordinate,sy%lattice_vector,coulcut,nl,lt%verbose,myRank,numRanks) -#else - !call gpmdcov_build_nlist_sedacs(sy%coordinate,sy%lattice_vector,coulcut,nl,lt%verbose,myRank,numRanks) - !call gpmdcov_build_nlist_sparse_v2(sy%coordinate,sy%lattice_vector,coulcut,nl,lt%verbose,myRank,numRanks) - call gpmdcov_build_nlist_sparse_v2(sy%coordinate,sy%lattice_vector,coulcut,nl,lt%verbose,myRank,numRanks) -#endif + ! if(any(nl2%nrnnstruct.ne.nl%nrnnstruct))then ! write(*,*)"DEBUG: nrnnstruct not equal" ! do k = 1,size(nl%nrnnstruct) From 987307871f43d80413cd7f382d3e6f13c63e9614 Mon Sep 17 00:00:00 2001 From: Mike Wall Date: Mon, 1 Jun 2026 12:01:17 -0700 Subject: [PATCH 15/16] Debug print_mdstep --- examples/gpmdk/src/gpmdcov_mdloop.F90 | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/examples/gpmdk/src/gpmdcov_mdloop.F90 b/examples/gpmdk/src/gpmdcov_mdloop.F90 index 885dd7e..0f56c05 100644 --- a/examples/gpmdk/src/gpmdcov_mdloop.F90 +++ b/examples/gpmdk/src/gpmdcov_mdloop.F90 @@ -617,7 +617,7 @@ end function cudaProfilerStop if (.not.first_substep_taken)then print_mdstep = print_mdstep + 1 write(*,'(A35,I15,A1,F18.5,A1,ES12.5,A1,ES12.5,A1,ES12.5)')"Mdstep, Energy, Egap, Resnorm, Temp", & - &print_mdstep-gpmdt%minimization_steps," ", Energy," ", egap_glob," ", resnorm," ", Temp + &print_mdstep," ", Energy," ", egap_glob," ", resnorm," ", Temp if (half_timestep_flag)then write(*,*) "WARNING: Two half timesteps were performed for step ", print_mdstep endif @@ -756,13 +756,13 @@ end function cudaProfilerStop #endif if(gpmdt%writetraj .and. myRank == 1 .and. mdstep.ge.gpmdt%minimization_steps .and. first_substep_taken .eqv. .false.)then if((gpmdt%traj_format .eq. "XYZ").and. & - (mod(print_mdstep-gpmdt%minimization_steps,gpmdt%writetreach).eq.0.or. & - (print_mdstep-gpmdt%minimization_steps).eq.1))then - call prg_write_trajectory(sy,print_mdstep-gpmdt%minimization_steps,gpmdt%writetreach,& + (mod(print_mdstep,gpmdt%writetreach).eq.0.or. & + (print_mdstep).eq.1))then + call prg_write_trajectory(sy,print_mdstep,gpmdt%writetreach,& <%timestep,adjustl(trim(lt%jobname))//"_trajectory","xyz") call prg_write_system(sy,adjustl(trim(lt%jobname))//"_latest","pdb") else - call prg_write_trajectory(sy,print_mdstep-gpmdt%minimization_steps,gpmdt%writetreach,& + call prg_write_trajectory(sy,print_mdstep,gpmdt%writetreach,& <%timestep,adjustl(trim(lt%jobname))//"_trajectory","pdb") endif endif @@ -780,7 +780,7 @@ end function cudaProfilerStop ! Save MD state each 120 steps if(gpmdt%dumpeach .gt. 0)then - if(mod(print_mdstep-gpmdt%minimization_steps,gpmdt%dumpeach) == 0)call gpmdcov_dump() + if(mod(print_mdstep,gpmdt%dumpeach) == 0)call gpmdcov_dump() endif if(mdstep.eq.gpmdt%minimization_steps)then From 47ba959392b284d0e2904b37374bee30075b88ab Mon Sep 17 00:00:00 2001 From: Mike Wall Date: Mon, 1 Jun 2026 12:28:27 -0700 Subject: [PATCH 16/16] Fix timestep for .pdb trajectory output --- examples/gpmdk/src/gpmdcov_mdloop.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/gpmdk/src/gpmdcov_mdloop.F90 b/examples/gpmdk/src/gpmdcov_mdloop.F90 index 0f56c05..0f4f49d 100644 --- a/examples/gpmdk/src/gpmdcov_mdloop.F90 +++ b/examples/gpmdk/src/gpmdcov_mdloop.F90 @@ -759,11 +759,11 @@ end function cudaProfilerStop (mod(print_mdstep,gpmdt%writetreach).eq.0.or. & (print_mdstep).eq.1))then call prg_write_trajectory(sy,print_mdstep,gpmdt%writetreach,& - <%timestep,adjustl(trim(lt%jobname))//"_trajectory","xyz") + &user_timestep,adjustl(trim(lt%jobname))//"_trajectory","xyz") call prg_write_system(sy,adjustl(trim(lt%jobname))//"_latest","pdb") else call prg_write_trajectory(sy,print_mdstep,gpmdt%writetreach,& - <%timestep,adjustl(trim(lt%jobname))//"_trajectory","pdb") + &user_timestep,adjustl(trim(lt%jobname))//"_trajectory","pdb") endif endif #ifdef USE_NVTX