diff --git a/examples/gpmdk/run/water/my_waterInput.in b/examples/gpmdk/run/water/my_waterInput.in new file mode 100644 index 00000000..c89034bc --- /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.2 + #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 +} 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 00000000..2ea8ad00 --- /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 8ef2bad9..0f4f49dd 100644 --- a/examples/gpmdk/src/gpmdcov_mdloop.F90 +++ b/examples/gpmdk/src/gpmdcov_mdloop.F90 @@ -36,7 +36,11 @@ subroutine gpmdcov_MDloop() real(dp) :: ke_tensor(3,3) real(dp) :: pressure_tensor(3,3) real(dp), allocatable :: saved_velocities(:,:) - integer :: total_steps + real(dp), allocatable :: saved_forces(:,:) + real(dp) :: user_timestep,this_maxdisp,user_half_timestep + real(dp), parameter :: maxdist = 0.02 + logical :: first_substep_taken,half_timestep_flag + integer :: total_steps, print_mdstep integer :: cuda_error logical :: newnl ! Indicates new neighbor list type(neighlist_type) :: nl2 @@ -74,11 +78,13 @@ 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 + 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) @@ -89,6 +95,19 @@ 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 ! if(mdstep < 0)then ! savets = lt%timestep @@ -122,6 +141,37 @@ 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) + write(*,*)"for mdstep ", mdstep, "this_maxdisp = ", this_maxdisp + 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. + write(*,*)"for mdstep ", mdstep, "reduced timestep = ", lt%timestep + + if (first_substep_taken) then + first_substep_taken = .false. + else + first_substep_taken = .true. + endif + + else + lt%timestep = user_timestep + half_timestep_flag = .false. + 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)))//" & &for (atom,axis) = ("//to_string(maxv_atom_axis(2))//","//to_string(maxv_atom_axis(1))//")",lt%verbose,myRank) @@ -145,7 +195,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 @@ -184,6 +235,7 @@ end function cudaProfilerStop write(*,*)i,sy%velocity(1,i),sy%velocity(2,i),sy%velocity(3,i) enddo endif + !> 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") @@ -356,7 +408,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) @@ -366,13 +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) -#endif + ! if(any(nl2%nrnnstruct.ne.nl%nrnnstruct))then ! write(*,*)"DEBUG: nrnnstruct not equal" ! do k = 1,size(nl%nrnnstruct) @@ -430,12 +478,14 @@ end function cudaProfilerStop #endif call gpmdcov_Part(2) + #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) + !> Reprg_initialize parts. mls_i = mls() call gpmdcov_msMem("gpmdcov_mdloop", "Before gpmdcov_InitParts",lt%verbose,myRank) @@ -457,7 +507,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 @@ -510,7 +561,7 @@ end function cudaProfilerStop #ifdef USE_NVTX call gpmdEndRange #endif - if(kernel%xlbolevel1.and.lt%doKernel)then + if(kernel%xlbolevel1.and.lt%doKernel)then allocate(n1(sy%nats)) if(mdstep > 1)then !sy%net_charge = n @@ -551,6 +602,7 @@ end function cudaProfilerStop mls_md1 = mls() call gpmdcov_msI("gpmdcov_MDloop","ResNorm = "//to_string(resnorm),lt%verbose,myRank) + if(myRank == 1)then if(mdstep.le.gpmdt%minimization_steps)then if(.not.gpmdt%anneal_graph)then @@ -561,10 +613,16 @@ 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 + ! 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 + if (half_timestep_flag)then + write(*,*) "WARNING: Two half timesteps were performed for step ", print_mdstep + endif + endif endif - !write(*,*)"Step, Energy, EGap, Resnorm", mdstep, Energy, egap_glob, resnorm endif #ifdef USE_NVTX call gpmdStartRange("EnergAndForces",7) @@ -696,16 +754,16 @@ 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,& - <%timestep,adjustl(trim(lt%jobname))//"_trajectory","xyz") + (mod(print_mdstep,gpmdt%writetreach).eq.0.or. & + (print_mdstep).eq.1))then + call prg_write_trajectory(sy,print_mdstep,gpmdt%writetreach,& + &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,mdstep-gpmdt%minimization_steps,gpmdt%writetreach,& - <%timestep,adjustl(trim(lt%jobname))//"_trajectory","pdb") + call prg_write_trajectory(sy,print_mdstep,gpmdt%writetreach,& + &user_timestep,adjustl(trim(lt%jobname))//"_trajectory","pdb") endif endif #ifdef USE_NVTX @@ -722,13 +780,15 @@ 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%dumpeach) == 0)call gpmdcov_dump() endif 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