*******> update.17: Author: Chuan Tian Date: Sep 9, 2019 Programs: pmemd, pmemd.MPI, pmemd.cuda and pmemd.cuda.MPI Description: This modifies the CMAP code for cpu and gpu, adds cuda test cases -------------------------------------------------------------------------------- src/pmemd/src/cmap.F90 | 1 - src/pmemd/src/gb_alltasks_setup.F90 | 8 + src/pmemd/src/master_setup.F90 | 5 +- src/pmemd/src/parallel.F90 | 3 + src/pmemd/src/pme_alltasks_setup.F90 | 9 +- src/pmemd/src/prmtop_dat.F90 | 100 +- src/pmemd/src/runfiles.F90 | 5 + src/pmemd/src/runmin.F90 | 5 + src/pmemd/src/ti.F90 | 12 +- test/cuda/Makefile | 3 +- test/cuda/ff19SB/Run | 72 + test/cuda/ff19SB/inpcrd | Bin 0 -> 7167112 bytes test/cuda/ff19SB/mdin | 13 + test/cuda/ff19SB/mdout.CPU.save | 403 + test/cuda/ff19SB/mdout.GPU_DPFP.save | 433 + test/cuda/ff19SB/mdout.GPU_SPFP.save | 439 + test/cuda/ff19SB/prmtop | 299831 ++++++++++++++++++++++++++++++++ 17 files changed, 301300 insertions(+), 42 deletions(-) diff --git src/pmemd/src/cmap.F90 src/pmemd/src/cmap.F90 index 69b53b33fc..21233b601d 100644 --- src/pmemd/src/cmap.F90 +++ src/pmemd/src/cmap.F90 @@ -455,7 +455,6 @@ subroutine generate_cmap_derivatives integer :: res,halfRes,twoRes double precision, allocatable,dimension(:) :: tmpy,tmpy2 - write(6,'(a)') '|CHARMM: Reticulating splines.' write(6,'(a)') '' do i=1,cmap_type_count diff --git src/pmemd/src/gb_alltasks_setup.F90 src/pmemd/src/gb_alltasks_setup.F90 index 76cb128704..183522414a 100644 --- src/pmemd/src/gb_alltasks_setup.F90 +++ src/pmemd/src/gb_alltasks_setup.F90 @@ -124,6 +124,8 @@ subroutine gb_alltasks_setup(num_ints, num_reals) if (charmm_active) then call angles_ub_setup(num_ints, num_reals, use_atm_map) call dihedrals_imp_setup(num_ints, num_reals, use_atm_map) + end if + if (charmm_active .or. cmap_term_count > 0) then call cmap_setup(num_ints, num_reals, use_atm_map) endif call nb14_setup(num_ints, num_reals, use_atm_map) @@ -145,6 +147,8 @@ subroutine gb_alltasks_setup(num_ints, num_reals) if (charmm_active) then call angles_ub_setup(num_ints, num_reals, use_atm_map) call dihedrals_imp_setup(num_ints, num_reals, use_atm_map) + endif + if (charmm_active .or. cmap_term_count > 0) then call cmap_setup(num_ints, num_reals, use_atm_map) endif call nb14_setup(num_ints, num_reals, use_atm_map) @@ -160,6 +164,10 @@ subroutine gb_alltasks_setup(num_ints, num_reals) num_ints = num_ints - size(gbl_cmap) * cmap_rec_ints deallocate(gbl_angle_ub, gbl_dihed_imp, gbl_cmap) + else if (cmap_term_count > 0) then + num_ints = num_ints - size(gbl_cmap) * cmap_rec_ints + + deallocate(gbl_cmap) endif deallocate(gbl_angle, gbl_dihed) diff --git src/pmemd/src/master_setup.F90 src/pmemd/src/master_setup.F90 index f615d350f9..efbebff98c 100644 --- src/pmemd/src/master_setup.F90 +++ src/pmemd/src/master_setup.F90 @@ -31,7 +31,7 @@ contains subroutine master_setup(num_ints, num_reals, new_stack_limit, terminal_flag) use axis_optimize_mod - use charmm_mod, only : charmm_active + !use charmm_mod, only : charmm_active use cmap_mod, only : generate_cmap_derivatives use cit_mod use constantph_mod, only : cnstph_read, cph_igb @@ -329,7 +329,8 @@ subroutine master_setup(num_ints, num_reals, new_stack_limit, terminal_flag) if (icnste .gt. 0) call cnste_read(num_ints, num_reals) - if (charmm_active) call generate_cmap_derivatives + !if (charmm_active) call generate_cmap_derivatives + if (cmap_term_count > 0) call generate_cmap_derivatives ! If the user has requested NMR restraints, do a cursory read of the ! restraints file(s) now to determine the amount of memory necessary diff --git src/pmemd/src/parallel.F90 src/pmemd/src/parallel.F90 index 99ae9b9557..f0bc61db5e 100644 --- src/pmemd/src/parallel.F90 +++ src/pmemd/src/parallel.F90 @@ -2177,6 +2177,7 @@ subroutine do_atm_distribution(atm_cnt, num_ints, num_reals, & use charmm_mod, only : charmm_active use mdin_ctrl_dat_mod, only : ips, nmropt use ti_mod + use prmtop_dat_mod implicit none @@ -2233,6 +2234,8 @@ subroutine do_atm_distribution(atm_cnt, num_ints, num_reals, & if (charmm_active) then call angles_ub_setup(num_ints, num_reals, use_atm_map) call dihedrals_imp_setup(num_ints, num_reals, use_atm_map) + end if + if (cmap_term_count > 0) then call cmap_setup(num_ints, num_reals, use_atm_map) endif call nb14_setup(num_ints, num_reals, use_atm_map) diff --git src/pmemd/src/pme_alltasks_setup.F90 src/pmemd/src/pme_alltasks_setup.F90 index 4340716fc8..86a329ac5c 100644 --- src/pmemd/src/pme_alltasks_setup.F90 +++ src/pmemd/src/pme_alltasks_setup.F90 @@ -260,6 +260,8 @@ endif if (charmm_active) then call angles_ub_setup(num_ints, num_reals, use_atm_map) call dihedrals_imp_setup(num_ints, num_reals, use_atm_map) + end if + if (cmap_term_count > 0) then call cmap_setup(num_ints, num_reals, use_atm_map) endif @@ -271,9 +273,12 @@ endif if (charmm_active) then num_ints = num_ints - size(gbl_angle_ub) * angle_ub_rec_ints num_ints = num_ints - size(gbl_dihed_imp) * dihed_imp_rec_ints - num_ints = num_ints - size(gbl_cmap) * cmap_rec_ints - deallocate(gbl_angle_ub, gbl_dihed_imp, gbl_cmap) + deallocate(gbl_angle_ub, gbl_dihed_imp) + endif + if (cmap_term_count > 0) then + num_ints = num_ints - size(gbl_cmap) * cmap_rec_ints + deallocate(gbl_cmap) endif deallocate(gbl_angle, gbl_dihed) diff --git src/pmemd/src/prmtop_dat.F90 src/pmemd/src/prmtop_dat.F90 index ef7f7d6255..afbf486478 100644 --- src/pmemd/src/prmtop_dat.F90 +++ src/pmemd/src/prmtop_dat.F90 @@ -412,17 +412,21 @@ if(.not. usemidpoint) & type = 'CHARMM_NUM_IMPR_TYPES' call nxtsec(prmtop, mdout, 0, fmtin, type, fmt, errcode) read(prmtop,fmt) nimprtyp - + end if !Read CHARMM CMAP fmtin = '(2I8)' type = 'CHARMM_CMAP_COUNT' call nxtsec(prmtop, mdout, 1, fmtin, type, fmt, errcode) + if (errcode .ne. 0) then + type = 'CMAP_COUNT' + call nxtsec(prmtop, mdout, 1, fmtin, type, fmt, errcode) + end if + if (errcode == 0) then read(prmtop,fmt) cmap_term_count, cmap_type_count end if - end if @@ -926,41 +930,54 @@ if(.not. usemidpoint) & read(prmtop, fmt) (gbl_imp_phase(i), i = 1, nimprtyp) + end if + !CMAP - if ( cmap_term_count > 0 ) then - fmtin = '(9I8)' - type = 'CHARMM_CMAP_INDEX' + if ( cmap_term_count > 0 ) then + fmtin = '(9I8)' + type = 'CHARMM_CMAP_INDEX' + call nxtsec(prmtop, mdout, 1, fmtin, type, fmt, errcode) + if (errcode .ne. 0) then + type = 'CMAP_INDEX' call nxtsec(prmtop, mdout, 0, fmtin, type, fmt, errcode) - read(prmtop, fmt) (gbl_cmap(i)%atm_i,& - gbl_cmap(i)%atm_j,& - gbl_cmap(i)%atm_k,& - gbl_cmap(i)%atm_l,& - gbl_cmap(i)%atm_m,& - gbl_cmap(i)%parm_idx,& - i=1, cmap_term_count) - - fmtin = '(20I4)' - type = 'CHARMM_CMAP_RESOLUTION' + end if + read(prmtop, fmt) (gbl_cmap(i)%atm_i,& + gbl_cmap(i)%atm_j,& + gbl_cmap(i)%atm_k,& + gbl_cmap(i)%atm_l,& + gbl_cmap(i)%atm_m,& + gbl_cmap(i)%parm_idx,& + i=1, cmap_term_count) + + fmtin = '(20I4)' + type = 'CHARMM_CMAP_RESOLUTION' + call nxtsec(prmtop, mdout, 1, fmtin, type, fmt, errcode) + if (errcode .ne. 0) then + type = 'CMAP_RESOLUTION' call nxtsec(prmtop, mdout, 0, fmtin, type, fmt, errcode) - read(prmtop, fmt) (gbl_cmap_res(i), i=1, cmap_type_count) - - do i=1,cmap_type_count - write(word,'(i2.2)')i - fmtin = '(8(F9.5))' - type = "CHARMM_CMAP_PARAMETER_" // word + end if + read(prmtop, fmt) (gbl_cmap_res(i), i=1, cmap_type_count) + + do i=1,cmap_type_count + write(word,'(i2.2)')i + fmtin = '(8(F9.5))' + type = "CHARMM_CMAP_PARAMETER_" // word + call nxtsec(prmtop, mdout, 1, fmtin, type, fmt, errcode) + if (errcode .ne. 0) then + type = "CMAP_PARAMETER_" // word call nxtsec(prmtop, mdout, 0, fmtin, type, fmt, errcode) + end if - read(prmtop, fmt) gbl_cmap_grid(i, 1:gbl_cmap_res(i), & - 1:gbl_cmap_res(i) ) + read(prmtop, fmt) gbl_cmap_grid(i, 1:gbl_cmap_res(i), & + 1:gbl_cmap_res(i) ) - !Also take advantage of this loop to set the step size - gbl_cmap_grid_step_size(i) = 360/gbl_cmap_res(i) - enddo + !Also take advantage of this loop to set the step size + gbl_cmap_grid_step_size(i) = 360/gbl_cmap_res(i) + enddo - end if ! cmap_term_count > 0 + end if ! cmap_term_count > 0 - end if ! Read the dihedral information: @@ -1670,7 +1687,8 @@ subroutine bcast_amber_prmtop_dat call mpi_bcast(gbl_imp_phase, nimprtyp, mpi_double_precision, 0, & pmemd_comm, err_code_mpi) end if - + end if + ! CMAP TODO: 24 is a hard coded hack if ( cmap_type_count > 0 ) then call mpi_bcast(gbl_cmap_res, cmap_type_count, mpi_integer, 0, & @@ -1691,8 +1709,6 @@ subroutine bcast_amber_prmtop_dat end if - end if - @@ -1735,6 +1751,8 @@ subroutine bcast_amber_prmtop_dat call mpi_bcast(gbl_dihed_imp, size(gbl_dihed_imp) * bytes_per_unit, mpi_byte, 0, & pmemd_comm, err_code_mpi) + end if + if (cmap_term_count > 0) then !CMAP call get_bytesize(gbl_cmap(1), gbl_cmap(2), bytes_per_unit) @@ -1813,6 +1831,10 @@ subroutine alloc_amber_prmtop_mem(num_ints, num_reals) gbl_ub_r0(nubtypes),& gbl_imp_pk(nimprtyp),& gbl_imp_phase(nimprtyp),& + stat = alloc_failed) + end if + if (cmap_term_count > 0) then + allocate( & gbl_cmap_res(cmap_type_count), & gbl_cmap_grid_step_size(cmap_type_count), & gbl_cmap_grid(cmap_type_count,24,24),& !24 is a hardcoded hack @@ -1836,7 +1858,10 @@ subroutine alloc_amber_prmtop_mem(num_ints, num_reals) if (charmm_active) then num_reals = num_reals + size(gbl_cn114) + size(gbl_cn214) + & size(gbl_ub_rk) + size(gbl_ub_r0) + size(gbl_imp_pk) + & - size(gbl_imp_phase) + & + size(gbl_imp_phase) + end if + if (cmap_term_count > 0) then + num_reals = num_reals + & size(gbl_cmap_res) + size(gbl_cmap_grid_step_size) + & size(gbl_cmap_grid) + size(gbl_cmap_dPhi) + & size(gbl_cmap_dPsi) + size(gbl_cmap_dPhi_dPsi) @@ -1871,6 +1896,12 @@ subroutine alloc_amber_prmtop_mem(num_ints, num_reals) allocate( & gbl_angle_ub(gbl_angle_ub_allocsize), & gbl_dihed_imp(gbl_dihed_imp_allocsize), & + stat = alloc_failed) + + if (alloc_failed .ne. 0) call setup_alloc_error + end if + if (charmm_active .or. cmap_term_count > 0) then + allocate( & gbl_cmap(gbl_cmap_allocsize), & stat = alloc_failed) @@ -1884,7 +1915,10 @@ subroutine alloc_amber_prmtop_mem(num_ints, num_reals) if (charmm_active) then num_ints = num_ints + size(gbl_angle_ub) * angle_ub_rec_ints + & - size(gbl_dihed_imp) * dihed_imp_rec_ints + & + size(gbl_dihed_imp) * dihed_imp_rec_ints + end if + if (cmap_term_count > 0) then + num_ints = num_ints + & size(gbl_cmap) * cmap_rec_ints end if diff --git src/pmemd/src/runfiles.F90 src/pmemd/src/runfiles.F90 index 5550532447..68200d1523 100644 --- src/pmemd/src/runfiles.F90 +++ src/pmemd/src/runfiles.F90 @@ -640,6 +640,7 @@ subroutine prntmd(nstep, total_nstlim, time, si, fac, iout7, rms, mdloop) use mdin_ewald_dat_mod use nmr_calls_mod use parallel_dat_mod + use prmtop_dat_mod, only : cmap_term_count #ifdef MPI use remd_mod, only : remd_method, remd_dimension, group_num, replica_indexes #endif @@ -702,6 +703,8 @@ subroutine prntmd(nstep, total_nstlim, time, si, fac, iout7, rms, mdloop) if (charmm_active) then write(mdout, 9040) si(si_angle_ub_ene), si(si_dihedral_imp_ene), & si(si_cmap_ene) + else if (cmap_term_count > 0) then + write(mdout, 9040) 0.d0,0.d0, si(si_cmap_ene) endif write(mdout, 9049) si(si_vdw_14_ene), si(si_elect_14_ene), si(si_vdw_ene) @@ -859,6 +862,8 @@ endif if (charmm_active) then write(mdinfo, 9040) si(si_angle_ub_ene), si(si_dihedral_imp_ene), & si(si_cmap_ene) + else if (cmap_term_count > 0) then + write(mdinfo, 9040) 0.d0,0.d0,si(si_cmap_ene) endif write(mdinfo, 9049) si(si_vdw_14_ene), si(si_elect_14_ene), si(si_vdw_ene) diff --git src/pmemd/src/runmin.F90 src/pmemd/src/runmin.F90 index 5a0401368f..479535c057 100644 --- src/pmemd/src/runmin.F90 +++ src/pmemd/src/runmin.F90 @@ -930,6 +930,7 @@ subroutine printe(nstep, dff, fdmax, si, iatmax, labmax) use state_info_mod use charmm_mod, only : charmm_active use ti_mod + use prmtop_dat_mod, only: cmap_term_count implicit none ! Formal arguments: @@ -958,6 +959,8 @@ subroutine printe(nstep, dff, fdmax, si, iatmax, labmax) if (charmm_active) then write(mdout, 9039) si(si_angle_ub_ene), si(si_dihedral_imp_ene), & si(si_cmap_ene) + else if (cmap_term_count > 0) then + write(mdout, 9039) 0.d0,0.d0,si(si_cmap_ene) endif if (using_pme_potential) then @@ -1023,6 +1026,8 @@ subroutine printe(nstep, dff, fdmax, si, iatmax, labmax) if (charmm_active) then write(mdinfo, 9039) si(si_angle_ub_ene), si(si_dihedral_imp_ene), & si(si_cmap_ene) + else if (cmap_term_count > 0) then + write(mdinfo, 9039) 0.d0,0.d0,si(si_cmap_ene) endif if (using_pme_potential) then diff --git src/pmemd/src/ti.F90 src/pmemd/src/ti.F90 index 30659a08f2..bb5eb0a7ee 100644 --- src/pmemd/src/ti.F90 +++ src/pmemd/src/ti.F90 @@ -2365,6 +2365,7 @@ subroutine ti_print_ene(iti, sc_ener_in, sc_ener_aug_in) use file_io_dat_mod use charmm_mod, only : charmm_active + use prmtop_dat_mod, only: cmap_term_count implicit none @@ -2400,9 +2401,14 @@ subroutine ti_print_ene(iti, sc_ener_in, sc_ener_aug_in) sc_ener_in(si_vdw_ene) write (mdout,1200) sc_ener_in(si_elect_ene) - if(charmm_active) write (mdout,1400) sc_ener_in(si_angle_ub_ene), & - sc_ener_in(si_dihedral_imp_ene), & - sc_ener_in(si_cmap_ene) !charmm + if(charmm_active) then + write (mdout,1400) sc_ener_in(si_angle_ub_ene), & + sc_ener_in(si_dihedral_imp_ene), & + sc_ener_in(si_cmap_ene) !charmm + else if(cmap_term_count > 0) then + write(mdout,1400) 0, 0, & + sc_ener_in(si_cmap_ene) !charmm + end if write (mdout,1300) sc_ener_aug_in(ti_rest_dist_ene), & sc_ener_aug_in(ti_rest_ang_ene), & diff --git test/cuda/Makefile test/cuda/Makefile index 284fba627e..cca5a3194e 100644 --- test/cuda/Makefile +++ test/cuda/Makefile @@ -170,7 +170,8 @@ test.pmemd.cuda.pme: -cd lipid_npt_tests/ && ./Run_npt_isotropic.lipid14 $(PREC_MODEL) $(NETCDF) -cd lipid_npt_tests/ && ./Run_npt_anisotropic.lipid14 $(PREC_MODEL) $(NETCDF) -cd lipid_npt_tests/ && ./Run_npt_semiisotropic.lipid14 $(PREC_MODEL) $(NETCDF) - +#Newly implemented CMAP test for ff19SB + -cd ff19SB/ && ./Run $(PREC_MODEL) $(NETCDF) test.pmemd.cuda.remd: -cd remd && $(MAKE) test # -cd remd/rem_gb_2rep && ./Run.rem $(PREC_MODEL) $(NETCDF) diff --git test/cuda/ff19SB/Run test/cuda/ff19SB/Run new file mode 100755 index 0000000000..746e037ccb --- /dev/null +++ test/cuda/ff19SB/Run @@ -0,0 +1,72 @@ +#!/bin/csh -f +#TEST-PROGRAM pmemd.cuda +#TEST-DESCRIP TO_BE_DEtermined +#TEST-PURPOSE regression, basic +#TEST-STATE undocumented + + +if( ! $?DO_PARALLEL ) then + setenv DO_PARALLEL " " + if( $?TESTsander ) then + set sander = $TESTsander + else + set sander = ${AMBERHOME}/bin/pmemd.cuda_$1 + endif +else + if( $?TESTsander ) then + set sander = $TESTsander + else + set sander = ${AMBERHOME}/bin/pmemd.cuda_$1.MPI + endif +endif + + +cat > mdin <