*******> update.9 Author: Jinan Wang, Yinglong Miao Date: Jan 5, 2021 Programs: pmemd.cuda Description: Fixes a bug in LiGaMD & Pep-GaMD in pmemd.cuda -------------------------------------------------------------------------------- src/pmemd/src/cuda/gti_f95.cpp | 107 ++++++++++++++++++++++++++++++++++++++ src/pmemd/src/cuda/gti_f95.h | 4 ++ src/pmemd/src/extra_pnts_nb14.F90 | 16 +++++- src/pmemd/src/gamd.F90 | 25 +++++++++ src/pmemd/src/mdin_ctrl_dat.F90 | 12 +++-- src/pmemd/src/runmd.F90 | 99 ++++++++++++++++++++++++++--------- 6 files changed, 230 insertions(+), 33 deletions(-) diff --git src/pmemd/src/cuda/gti_f95.cpp src/pmemd/src/cuda/gti_f95.cpp index ef6d62b506..4712190a4d 100644 --- src/pmemd/src/cuda/gti_f95.cpp +++ src/pmemd/src/cuda/gti_f95.cpp @@ -295,6 +295,113 @@ extern "C" void gti_ti_nb_setup_(int* pNtypes, int iac[], int ico[], int* pNb14_ gpu->sim.pTINb142 = gpu->pbTINb142->_pDevData; gpu->sim.pTINb14ID = gpu->pbTINb14ID->_pDevData; } +//--------------------------------------------------------------------------------------------- +// gti_ti_nb_gamd_setup_: +// +// Arguments: +// pNtypes: +// iac: +// ico: +// pNb14_cnt: +// cit_nb14: +// gbl_one_scee: +// gbl_one_scnb: +// cn114: +// cn214: +//--------------------------------------------------------------------------------------------- +extern "C" void gti_ti_nb_gamd_setup_(int* pNtypes, int iac[], int ico[], int* pNb14_cnt, + int cit_nb14[][3], double gbl_one_scee[], + double gbl_one_scnb[], double cn114[], double cn214[]) +{ + PRINTMETHOD("gti_ti_nb_gamd_setup"); + + gpuContext gpu = theGPUContext::GetPointer(); + unsigned int ntypes = *pNtypes; + unsigned maxNumberTINBEntries = gpu->sim.numberTIAtoms * 3000; + gpu->pbTINBList = std::unique_ptr< GpuBuffer > (new GpuBuffer(maxNumberTINBEntries)); + + for (unsigned i = 0; i < maxNumberTINBEntries; i++) { + gpu->pbTINBList->_pSysData[i].x = -1; + gpu->pbTINBList->_pSysData[i].y = -1; + gpu->pbTINBList->_pSysData[i].z = -1; + gpu->pbTINBList->_pSysData[i].w = -1; + } + gpu->pbTINBList->Upload(); + gpu->sim.pTINBList = gpu->pbTINBList->_pDevData; + gpu->pbNumberTINBEntries = std::unique_ptr< GpuBuffer > (new GpuBuffer(1)); + gpu->pbNumberTINBEntries->_pSysData[0] = 0; + gpu->pbNumberTINBEntries->Upload(); + gpu->sim.pNumberTINBEntries = gpu->pbNumberTINBEntries->_pDevData; + + // 1:4 NB part + unsigned int nb14s = (gpu->ntf < 8) ? *pNb14_cnt : 0; + int(*temp)[4] = new int[nb14s][4]; + bool inRegion[2][2]; + unsigned counter = 0; + for (unsigned int i = 0; i < nb14s; i++) { + for (unsigned int j = 0; j < 2; j++) { + unsigned atom = cit_nb14[i][j] - 1; + inRegion[0][j] = (gpu->pbTIList->_pSysData[atom])>0; + inRegion[1][j] = (gpu->pbTIList->_pSysData[atom + gpu->sim.stride])>0; + } + temp[counter][3] = -1; + if (inRegion[0][0] || inRegion[0][1]) { + temp[counter][3] = 0; + } + if (inRegion[1][0] || inRegion[1][1]) { + temp[counter][3] = 1; + } + if (temp[counter][3] >= 0) { + for (unsigned int j = 0; j < 3; j++) { + temp[counter][j] = cit_nb14[i][j]; + } +// cit_nb14[i][2]=-1; + counter++; + } + } + nb14s = counter; + + // Copy 1-4 interactions + gpu->pbTINb141 = std::unique_ptr< GpuBuffer > (new GpuBuffer(nb14s)); + gpu->pbTINb142 = std::unique_ptr< GpuBuffer > (new GpuBuffer(nb14s)); + gpu->pbTINb14ID = std::unique_ptr< GpuBuffer > (new GpuBuffer(nb14s)); + if (gpu->ntf < 8) { + for (unsigned int i = 0; i < nb14s; i++) { + int parm_idx = temp[i][2] - 1; + unsigned atom0 = abs(temp[i][0]) - 1; + unsigned atom1 = abs(temp[i][1]) - 1; + + gpu->pbTINb141->_pSysData[i].x = gbl_one_scee[parm_idx]; + gpu->pbTINb141->_pSysData[i].y = gbl_one_scnb[parm_idx]; + int tt = ntypes * (iac[abs(temp[i][0]) - 1] - 1) + (iac[abs(temp[i][1]) - 1] - 1); + int nbtype = tt >= 0 ? ico[tt] - 1 : -1; + if (nbtype >= 0) { + gpu->pbTINb142->_pSysData[i].x = cn114[nbtype]; + gpu->pbTINb142->_pSysData[i].y = cn214[nbtype]; + gpu->pbTINb142->_pSysData[i].z = cn214[nbtype]/cn114[nbtype]; // B/A + } + else { + gpu->pbTINb142->_pSysData[i].x = ZeroF; + gpu->pbTINb142->_pSysData[i].y = ZeroF; + gpu->pbTINb142->_pSysData[i].z = ZeroF; + } + gpu->pbTINb14ID->_pSysData[i].x = atom0; + gpu->pbTINb14ID->_pSysData[i].y = atom1; + gpu->pbTINb14ID->_pSysData[i].z = temp[i][3]; + } + } + delete[] temp; + gpu->pbTINb141->Upload(); + gpu->pbTINb142->Upload(); + gpu->pbTINb14ID->Upload(); + + // Set constants + gpu->sim.numberTI14NBEntries = nb14s; + gpu->sim.pTINb141 = gpu->pbTINb141->_pDevData; + gpu->sim.pTINb142 = gpu->pbTINb142->_pDevData; + gpu->sim.pTINb14ID = gpu->pbTINb14ID->_pDevData; +} + //--------------------------------------------------------------------------------------------- // gti_lj1264_nb_setup_: diff --git src/pmemd/src/cuda/gti_f95.h src/pmemd/src/cuda/gti_f95.h index 95f9d86203..be07d82e15 100644 --- src/pmemd/src/cuda/gti_f95.h +++ src/pmemd/src/cuda/gti_f95.h @@ -71,6 +71,10 @@ extern "C" void gti_ti_nb_setup_(int* pNtypes, int iac[], int ico[], int* pNb14_ int cit_nb14[][3], double gbl_one_scee[], double gbl_one_scnb[], double cn114[], double cn214[]); +extern "C" void gti_ti_nb_gamd_setup_(int* pNtypes, int iac[], int ico[], int* pNb14_cnt, + int cit_nb14[][3], double gbl_one_scee[], + double gbl_one_scnb[], double cn114[], double cn214[]); + extern "C" void gti_lj1264_nb_setup_(char atm_igraph[][4], int* ntypes, int atm_iac[], int ico[], double gbl_cn6[]); diff --git src/pmemd/src/extra_pnts_nb14.F90 src/pmemd/src/extra_pnts_nb14.F90 index 093f279cb6..6469d060a5 100644 --- src/pmemd/src/extra_pnts_nb14.F90 +++ src/pmemd/src/extra_pnts_nb14.F90 @@ -493,6 +493,7 @@ subroutine nb14_setup(num_ints, num_reals, use_atm_map) do atm_i=1, natom if (ti_vdw_mask_list(atm_i).ne.0) iac_ti(atm_i)=0 end do + if(.not.(igamd.ge.6.and.igamd.le.11))then if (charmm_active) then call gti_ti_nb_setup(ntypes, iac_ti, typ_ico, & cit_nb14_cnt, cit_nb14, gbl_one_scee, gbl_one_scnb, & @@ -502,6 +503,17 @@ subroutine nb14_setup(num_ints, num_reals, use_atm_map) cit_nb14_cnt, cit_nb14, gbl_one_scee, gbl_one_scnb, & gbl_cn1, gbl_cn2) end if + else + if (charmm_active) then + call gti_ti_nb_gamd_setup(ntypes, iac_ti, typ_ico, & + cit_nb14_cnt, cit_nb14, gbl_one_scee, gbl_one_scnb, & + gbl_cn114, gbl_cn214) + else + call gti_ti_nb_gamd_setup(ntypes, iac_ti, typ_ico, & + cit_nb14_cnt, cit_nb14, gbl_one_scee, gbl_one_scnb, & + gbl_cn1, gbl_cn2) + end if + end if ! if(.not.(igamd.eq.12.or.igamd.eq.13.or.igamd.eq.14.or.igamd.eq.15.or.igamd.eq.16.or.igamd.eq.17 & ! .or.igamd.eq.18.or.igamd.eq.19.or.igamd.eq.21.or.(igamd.ge.110.and.igamd.le.120)))then do atm_i=1, natom @@ -513,11 +525,11 @@ subroutine nb14_setup(num_ints, num_reals, use_atm_map) if(igamd.eq.12.or.igamd.eq.13.or.igamd.eq.14.or.igamd.eq.15.or.igamd.eq.16.or.igamd.eq.17 & .or.igamd.eq.18.or.igamd.eq.19.or.igamd.eq.21.or.(igamd.ge.110.and.igamd.le.120))then if (charmm_active) then - call gti_ti_nb_setup(ntypes, iac_ti, typ_ico, & + call gti_ti_nb_gamd_setup(ntypes, iac_ti, typ_ico, & cit_nb14_cnt, cit_nb14, gbl_one_scee, gbl_one_scnb, & gbl_cn114, gbl_cn214) else - call gti_ti_nb_setup(ntypes, iac_ti, typ_ico, & + call gti_ti_nb_gamd_setup(ntypes, iac_ti, typ_ico, & cit_nb14_cnt, cit_nb14, gbl_one_scee, gbl_one_scnb, & gbl_cn1, gbl_cn2) end if diff --git src/pmemd/src/gamd.F90 src/pmemd/src/gamd.F90 index d15d01a5ce..ffc4089706 100644 --- src/pmemd/src/gamd.F90 +++ src/pmemd/src/gamd.F90 @@ -603,6 +603,16 @@ subroutine write_gamd_weights(ntwx,total_nstep) end subroutine write_gamd_weights +pure double precision function r2(v) + + implicit none + + double precision, intent(in) :: v(3) + + r2 = v(1)**2 + v(2)**2 + v(3)**2 + +end function r2 + pure double precision function norm3(v) implicit none @@ -704,4 +714,19 @@ function PBC_distance(r1, r2) result(value) end function PBC_distance +double precision function msd(lcrd0, lcrd1) + + implicit none + + integer :: i + double precision, intent(in) :: lcrd0(:,:), lcrd1(:,:) + + msd = 0.0d0 + do i = 1, size(lcrd0,2) + msd = msd + r2(lcrd1(:,i) - lcrd0(:,i)) + end do + msd = msd/size(lcrd0,2) + +end function msd + end module gamd_mod diff --git src/pmemd/src/mdin_ctrl_dat.F90 src/pmemd/src/mdin_ctrl_dat.F90 index 29eddc2f0e..8a2aab9f9b 100644 --- src/pmemd/src/mdin_ctrl_dat.F90 +++ src/pmemd/src/mdin_ctrl_dat.F90 @@ -37,7 +37,7 @@ use file_io_dat_mod ntcnste, ntrelaxe, ntwf, & iphmd, & !PHMD igamd,irest_gamd,iE,iEP,iED,iVm,igamdlag,ntcmd,nteb, & - ntcmdprep,ntebprep,ibblig,nlig,atom_p,atom_l,bgpro2atm,edpro2atm,& + ntcmdprep,ntebprep,ibblig,nlig,atom_p,atom_l,ntmsd,nftau,bgpro2atm,edpro2atm,& plumed, & ! PLUMED w_amd, emil_do_calc,isgld,isgsta,isgend,fixcom, & tishake, emil_sc,iemap, lj1264, efn,mcwat,mcint,& @@ -104,13 +104,13 @@ use file_io_dat_mod ifmbar_net, nebfreq , ramdint, ramdboostfreq, & !145 plumed, & !146 PLUMED reweight, iextpot, & !148 - ionstepvelocities,iEP,iED,ibblig,nlig,atom_p,atom_l !155 + ionstepvelocities,iEP,iED,ibblig,nlig,atom_p,atom_l,ntmsd,nftau !155 ! This variable must be set to the total number of integers in the ! mdin_ctrl_int common block. ! NOTE: ihwtnm adds 2 to the mdin_ctrl_int_cnt tally: add 1 to tally ! in comments above - integer, parameter :: mdin_ctrl_int_cnt = 156 + integer, parameter :: mdin_ctrl_int_cnt = 158 save :: / mdin_ctrl_int / @@ -353,7 +353,7 @@ use file_io_dat_mod EthreshD_w, alphaD_w, EthreshP_w, alphaP_w, emil_sc, lj1264, & igamd,irest_gamd,iE,iED,iEP,iVm,igamdlag,kD,kP,ntcmd,nteb,ntcmdprep,ntebprep, & VmaxD,VminD,VavgD,k0D,VmaxP,VminP,VavgP,k0P,ibblig,nlig,atom_p,atom_l,dblig,bgpro2atm,edpro2atm, & - sigma0D,sigma0P,sigmaVD,sigmaVP,cutoffD,cutoffP,fswitch, & + sigma0D,sigma0P,sigmaVD,sigmaVP,cutoffD,cutoffP,fswitch,ntmsd,nftau, & efx,efy,efz,efn,effreq,efphase,infe,irandom, & mbar_lambda, mbar_states, mask_from_ref, & gti_eq_nstep, gti_init_vel, gti_syn_mass, gremd_acyc, & @@ -471,7 +471,9 @@ subroutine init_mdin_ctrl_dat(remd_method, rremd_type) nlig = 0 atom_p = 0 atom_l = 0 - dblig = 4.0 + dblig = 1.0d99 + ntmsd = 0 + nftau = 0 ! IPS defaults ips = 0 diff --git src/pmemd/src/runmd.F90 src/pmemd/src/runmd.F90 index da6345de2b..531170df5b 100644 --- src/pmemd/src/runmd.F90 +++ src/pmemd/src/runmd.F90 @@ -275,9 +275,10 @@ subroutine runmd(atm_cnt, crd, mass, frc, vel, last_vel) logical :: update_gamd double precision :: VmaxDt,VminDt,VavgDt,sigmaVDt double precision :: VmaxPt,VminPt,VavgPt,sigmaVPt - double precision :: pot_ene_ll, pot_ene_l_pe, pot_ene_notl - integer :: blig,blig0,atm_i,atm_j - double precision :: dlig,pos_tmp(3) + double precision :: pot_ene_ll, pot_ene_l_pe, pot_ene_notl ! not used + integer :: blig,blig_min,blig0,atm_i,atm_j,nfmsd,ntwin,ift + double precision :: dlig,dmin,pos_tmp(3),fr(3) + double precision, allocatable :: lcrd0(:,:,:),lcrd1(:,:,:) integer,save :: counts=0 @@ -503,6 +504,19 @@ subroutine runmd(atm_cnt, crd, mass, frc, vel, last_vel) VavgPt = 0.0d0 sigmaVPt = 0.d0 +! GaMD settings + if (igamd.gt.0 .and. ibblig.gt.0) then + blig0 = 1 + blig = 1 + open(unit=gamdlig,file=gamdlig_name,action='write') + write(gamdlig,'(a)')'#step ligand distance' + ! allocate dynamic arrays to calculate ligand MSD + if (ibblig.eq.2) then + nfmsd = ntmsd/ntwx + ntwin = ntmsd+nftau*ntwx + allocate(lcrd0(3,nftau,nlig), lcrd1(3,nftau,nlig)) + end if + end if if ((igamd.gt.0) .and. (ntcmd.gt.0)) then igamd = 0 if (master) write(mdout,'(/,a,i10,/)') & @@ -511,12 +525,6 @@ subroutine runmd(atm_cnt, crd, mass, frc, vel, last_vel) call gpu_igamd_update(igamd) #endif end if - if (igamd.gt.0 .and. ibblig.gt.0) then - blig0 = 1 - blig = 1 - open(unit=gamdlig,file=gamdlig_name,action='write') - write(gamdlig,'(a)')'#step ligand_bound distance' - end if #if defined(MPI) && defined(CUDA) if (ineb .gt. 0 .and. (beadid .eq. 1 .or. beadid .eq. neb_nbead)) then @@ -4736,27 +4744,62 @@ DBG_ARRAYS_TIME_STEP(nstep) end if ! swap atomic coordinates of TI selected ligand residue with those of the bound ligand - if ( (mod(total_nstep, ntwx).eq.0) .and. (total_nstep.gt.(ntcmd+nteb)) .and. onstep) then - if (ibblig.gt.0 .and. igamd.eq.11) then + if (ibblig.gt.0 .and. onstep) then + if ( (mod(total_nstep, ntwx).eq.0) ) then #ifdef CUDA call gpu_download_crd(crd) #endif - do i = 1, nlig - atm_i = ti_atm_lst(1,1) + ti_ti_atm_cnt(1)*(i-1) + atom_l-1 - if (ntb.eq.0) then - dlig = distance(crd(1:3,atom_p), crd(1:3,atm_i)) - else - dlig = PBC_distance(crd(1:3,atom_p), crd(1:3,atm_i)) - end if - if (dlig.le.dblig) then - blig = i - write(mdout,'(a,2I10,F10.3)')'| (step, ligand_bound, distance): ', total_nstep, blig, dlig - write(gamdlig,'(2I10,F10.3)') total_nstep, blig, dlig - end if - end do + if ( ibblig.eq.2 ) then + do i = 1, nlig + atm_i = ti_atm_lst(1,1) + ti_ti_atm_cnt(1)*(i-1) + atom_l-1 + pos_tmp = crd(1:3,atm_i) + ! use PBC wrapped atomic coordinate + ! call PBC_r2f(crd(1:3,atm_i), fr) + ! call PBC_f2r(fr, pos_tmp, 0, 0, 0) + + ift = mod(total_nstep, ntwin)/ntwx + if (ift.lt.nftau) lcrd0(1:3,ift+1,i) = pos_tmp + if (ift.ge.nfmsd) lcrd1(1:3,ift-nfmsd+1,i) = pos_tmp + end do + end if + + ! calculate protein-ligand atom distance + dmin = 1.0d99 + if (ibblig.eq.1) then + do i = 1, nlig + atm_i = ti_atm_lst(1,1) + ti_ti_atm_cnt(1)*(i-1) + atom_l-1 + if (ntb.eq.0) then + dlig = distance(crd(1:3,atom_p), crd(1:3,atm_i)) + else + dlig = PBC_distance(crd(1:3,atom_p), crd(1:3,atm_i)) + end if + if (dmin.gt.dlig) then + dmin = dlig + blig_min = i + end if + end do + if (dmin.le.dblig) then + blig = blig_min + write(mdout,'(a,2I10,F14.6)')'| GaMD step, ligand, distance_min: ', total_nstep, blig, dmin + write(gamdlig,'(2I10,F14.6)') total_nstep, blig, dmin + end if + ! calculate ligand MSD + else if (ibblig.eq.2 .and. (mod(total_nstep, ntwin).eq.0)) then + do i = 1, nlig + dlig = msd(lcrd0(:,:,i), lcrd1(:,:,i)) + if (dmin.gt.dlig) then + dmin = dlig + blig_min = i + end if + end do + if (dmin.le.dblig) then + blig = blig_min + write(mdout,'(a,2I10,F14.6)')'| GaMD step, ligand, MSD_min: ', total_nstep, blig, dmin + write(gamdlig,'(2I10,F14.6)') total_nstep, blig, dmin + end if + end if if (blig.ne.blig0) then - ! swap ligand coordinates - write(mdout,'(a,2I10)')'| swap coordinates, forces and velocities of ligands: ', blig0, blig + write(mdout,'(a,2I10)')'| GaMD swap coordinates, forces and velocities of ligands: ', blig0, blig #ifdef CUDA call gpu_download_vel(vel) call gpu_download_frc(frc) @@ -5052,6 +5095,10 @@ DBG_ARRAYS_TIME_STEP(nstep) if (ti_mode .ne. 0) call ti_cleanup if (igamd.gt.0 .and. ibblig.gt.0) close(gamdlig) + if (igamd.gt.0 .and. ibblig.eq.2) then + if (allocated(lcrd0)) deallocate(lcrd0) + if (allocated(lcrd1)) deallocate(lcrd1) + end if call update_time(runmd_time)