! #include "copyright.h" #include "dprec.h" #include "assert.h" !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+ main driver routine to compute energies and forces subroutine force(xx,ix,ih,ipairs,x,f,ener,vir, & fs, rborn, reff, onereff, qsetup,qpsander,do_list_update ) use genborn use poisson_boltzmann, only : pb_force use dispersion_cavity, only : npopt, np_force use trace use stack use qmmm_module, only : qmmm_nml,qmmm_struct, qm2_struct, & qmmm_mpi, qmewald use constants, only : zero use relax_mat use ew_recip use parms, only:cn1,cn2,asol,bsol,pk,rk,tk,numbnd,numang,nptra,nphb,nimprp use nblist, only:nonbond_list, & a,b,c,alpha,beta,gamma !only used in qmmm call #ifdef PSANDER use spatial, only: nb2h,ib2h,jb2h,icb2h, & nb2a,ib2a,jb2a,icb2a, & nt2h,it2h,jt2h,kt2h,ict2h, & nt2a,it2a,jt2a,kt2a,ict2a, & nphi2h,ip2h,jp2h,kp2h,lp2h,icp2h, & nphi2a,ip2a,jp2a,kp2a,lp2a,icp2a, & max_dcmp_partners,ntskfrcrcv,ntskfrcsnd, & maxf_sndsiz, maxf_rcvsiz, & known_atm,heavy_atoms,h_atoms,nmyatm,spmylist #endif use amoeba_interface,only: AM_VAL_eval,AM_NonBond_eval use amoeba_mdin, only : iamoeba #ifdef PIMD use pimd_vars, only: natomCL, nbead, equal_part, ineb, Epot_external, & Epot_spring, Epot_deriv use qmmm_module, only: qm_gb,qmmm_pimd_first_call, & pimd_den_matrix, pimd_fock_matrix, pimd_old_den_matrix use pimd_vars, only: nrg_bnd, nrg_ang, nrg_dih, nrg_ele, nrg_vdw, nrg_egb, & nrg_scf,nrg_all,nrg_tmp use pimd_vars, only: xrep, frep, vrep, mmchg, qmchg, empot, reff_rep, & onereff_rep, iqmatoms, atomtype, iqmatomicnumbers, atommask, & mm_link_mask, link_pairs #endif #if defined(MPI) use evb_amber, only: xnrg, xq, xf, morsify_initialized, evb_step use evb_parm, only: natm3, nmorse,evb_pimd_debug, ntw_evb use evb_data, only: evb_frc_type, evb_frc use ntrfc_evb, only: evb_ntrfc # if defined(PIMD) use evb_pimd, only: evb_begin, evb_end & , pie, piq, pif, nrg_bead, fpimd # endif /* PIMD */ #endif /* MPI */ implicit none character(kind=1,len=5) :: routine="force" integer ier integer ipairs(*) _REAL_ xx(*) integer ix(*) character(len=4) ih(*) _REAL_ fs(*),rborn(*),reff(*),dvdl _REAL_, intent(out) :: onereff(*) # include "def_time.h" # include "ew_frc.h" # include "ew_cntrl.h" # include "extra_pts.h" # include "parallel.h" #ifdef MPI # include "ew_parallel.h" # ifdef MPI_DOUBLE_PRECISION # undef MPI_DOUBLE_PRECISION # endif # include "mpif.h" integer ierr, gb_rad_mpistart # ifdef CRAY_PVP # define MPI_DOUBLE_PRECISION MPI_REAL8 # endif #endif logical belly,nocrst # include "md.h" # include "pb_md.h" # include "box.h" # include "nmr.h" # include "memory.h" # include "files.h" # include "extra.h" # include "tgtmd.h" # include "flocntrl.h" # include "les.h" # ifdef LES # include "rem.h" # endif # include "sgld.h" integer istart,iend _REAL_ evdwex, eelex, virex(3,3) _REAL_ escf logical qsetup,do_list_update,qpsander _REAL_ enmr(3),devdis(4),devang(4),devtor(4),entr,ecap _REAL_ x(*),f(*),ene(30),vir(4) _REAL_ ener(*) ! offsets in this ener array = offsets in runmd ener - 22 save ene integer i,m,nttyp,i3 _REAL_ virvsene,eelt,epol,esurf,edisp _REAL_ epolar,aveper,aveind,avetot,emtot,dipiter,dipole_temp integer l_r2x,l_rjx,l_tmp1,l_tmp2,l_tmp3,l_tmp4,l_tmp5 integer l_tmp6,l_tmp7,l_tmp8,l_jj,l_skipv, l_kvls,l_jvls,l_psi integer l_da,l_sumd integer mmoffset, mmoffset3, qmoffset integer newbalance save newbalance #ifdef PSANDER integer tmpene,tmpfrc,tmpf2 #else !dummy versions of arrays used in psander integer, dimension(1):: known_atm,heavy_atoms,h_atoms, & ib2h,jb2h,icb2h,ib2a,jb2a,icb2a, & it2h,jt2h,kt2h,ict2h,it2a,jt2a,kt2a,ict2a, & ip2h,jp2h,kp2h,lp2h,icp2h, & ip2a,jp2a,kp2a,lp2a,icp2a, & spmylist integer :: nb2h,nb2a,nt2h,nt2a,nphi2h,nphi2a,nmyatm ! call bond(nb2h,ib2h,jb2h,icb2h, ! call bond(nb2a,ib2a,jb2a,icb2a, ! call angl(nt2h,it2h,jt2h,kt2h,ict2h ! call angl(nt2a,it2a,jt2a,kt2a,ict2a ! call ephi(nphi2h,ip2h,jp2h,kp2h,lp2h,icp2h ! call ephi(nphi2a,ip2a,jp2a,kp2a,lp2a,icp2a #endif _REAL_ :: etot #ifdef PIMD integer j,k integer :: nrg_tag, frc_tag _REAL_ :: spring_energy _REAL_ :: evb_bead(nbead), f_tmp(3,natomCL,1) _REAL_ :: evb_nrg, evb_nrg_sum #endif #if defined(MPI) integer :: status( MPI_STATUS_SIZE ) #endif /* MPI */ integer :: root, jobs_per_node, alloc_error, dealloc_error, n, nn integer :: ii, jj, kk call trace_enter( 'force' ) call timer_start(TIME_FORCE) ene(:) = ZERO belly = ibelly > 0 nocrst = .false. nttyp = ishft(ntypes*(ntypes+1),-1) !Division by 2 #ifdef MPI if (mpi_orig) then ! Check to see if we are done yet in mpi_orig case (tec3). ! This is done by monitoring the status of an integer notdone. ! If notdone .eq. 1 then we keep going. notdone is set to zero ! when we no longer want to call force(). This perhaps is not the ! most efegbficient means to implement this check... call mpi_bcast(notdone,1,mpi_integer,0,commsander,ierr) if (notdone /= 1) return ! Send copies of xyz coords, setbox common block, vir array ! and NTNB value to all nodes from master with a broadcast. if (numtasks > 1) then call mpi_bcast(box,BC_BOXR,MPI_DOUBLE_PRECISION,0, & commsander,ierr) call mpi_bcast(ntb,BC_BOXI,mpi_integer,0,commsander,ierr) call mpi_bcast(vir,3,MPI_DOUBLE_PRECISION,0, & commsander,ierr) call mpi_bcast(xx(lcrd),3*natom,MPI_DOUBLE_PRECISION, & 0,commsander,ierr) call mpi_bcast(ntnb,1,mpi_integer,0,commsander,ierr) if (iabs(ntb) >= 2) then call mpi_bcast(xx(l45),3*natom,MPI_DOUBLE_PRECISION, & 0,commsander,ierr) end if end if end if istart = iparpt(mytaskid) + 1 iend = iparpt(mytaskid+1) #else istart = 1 iend = natom #endif #ifdef PIMD nrg_bnd(:) = 0.0d0 nrg_ang(:) = 0.0d0 nrg_dih(:) = 0.0d0 nrg_ele(:) = 0.0d0 nrg_vdw(:) = 0.0d0 nrg_egb(:) = 0.0d0 nrg_scf(:) = 0.0d0 #endif ! ----- ZERO OUT THE ENERGIES AND FORCES ----- aveper=0.d0 aveind=0.d0 avetot=0.d0 dipiter=0.d0 dvdl=0.d0 dipole_temp=0.d0 enmr(1:3) = 0.d0 vir(1:4) = 0.d0 virvsene = 0.d0 f(1:3*natom+iscale) = 0.d0 if( igb == 0 .and. iyammp == 0 ) then ! (for GB: do all nonbondeds together below) call timer_start(TIME_NONBON) call timer_start(TIME_LIST) call nonbond_list(x,ix(i04),ix(i06),ix(i08),ix(i10), & ntypes,natom,xx,ix,ipairs,ntnb, & ix(ibellygp),belly,newbalance,cn1, & xx(lvel),xx(lvel2),ntp,xx(l45), qsetup,qpsander, & do_list_update, & known_atm,heavy_atoms,h_atoms) call timer_stop(TIME_LIST) call timer_stop(TIME_NONBON) end if !----------------------------------------------------------------- ! QMMM Link atom positioning in main amber coordinate array !----------------------------------------------------------------- !We need to adjust the link atoms positions !in the main amber coordinate array, replacing the mm link pair !atom coordinates temporarily. if (qmmm_nml%ifqnt) then call timer_start(TIME_QMMM) call timer_start(TIME_QMMMCOORDSX) if( qmmm_nml%idc == 0 ) call adj_mm_link_pair_crd(x) call timer_stop(TIME_QMMMCOORDSX) call timer_stop(TIME_QMMM) end if ! ---------------------------------------------------------------- ! Do weight changes, if requested ! ---------------------------------------------------------------- if (nmropt > 0) & call nmrcal(x,f,ih(m04),ih(m02),ix(i02),xx(lwinv),enmr,devdis, & devang,devtor,temp0,tautp,cut,ntb,xx(lnmr01), & ix(inmr02),xx(l95),31,6,rk,tk,pk,cn1, & cn2,asol,bsol,xx(l15),numbnd,numang,nptra-nimprp, & nimprp,nttyp,nphb,natom,natom,ntypes,nres, & rad,wel,radhb,welhb,rwell,isftrp,tgtrmsd,temp0les,-1,'WEIT') epolar = 0.d0 ! ----------------------------------------------------------------- ! EGB if igb>0 and /=6 then we need to calculate the GB radii for this ! structure ! ----------------------------------------------------------------- ! If we are doing qm_gb=2 then we need to calculate the GB radii ! before calling qm_mm if( igb > 0 .and. igb /= 6 .and. igb /=10 .and. & ( irespa < 2 .or. mod(irespa,nrespai) == 0) ) then #ifdef MPI gb_rad_mpistart = mytaskid+1 #endif call timer_start(TIME_EGB) call timer_start(TIME_GBRAD1) !If qmmm and then this will calculate us the radii for the !link atoms and not the mm link pair atoms. The initial radii used are those !of the mm link pair's atom type though. call egb_calc_radii(igb,natom,x,fs,reff, & onereff,fsmax,rgbmax, rborn, offset,gbalpha, & gbbeta,gbgamma,rbornstat,xx(l188),xx(l189), & xx(l186),xx(l187), gbneckscale, ncopy, rdt & #ifdef MPI ,gb_rad_mpistart & #endif ) call timer_stop(TIME_GBRAD1) call timer_stop(TIME_EGB) end if ! QM/MM Contributions are now calculated before the NON-Bond info. ! ---------------------------------------------------------------- ! Calculate the qm/mm contributions ! ---------------------------------------------------------------- if(qmmm_nml%ifqnt) then ! If we are doing periodic boundaries with QM/MM PME then we need to ! do the PME calculation twice. First here to get the potential at ! each QM atom due to the PME and then again after all the QM is done ! to get the MM-MM potential and all of the gradients. if(qmmm_nml%qm_pme) then ! Ewald force will put the potential into the qm_ewald%empot array. call timer_start(TIME_EWALD) call ewald_force(x,natom,ix(i04),ix(i06),ntypes, & xx(l15),cn1,cn2,asol,bsol,eelt,epolar, & f,xx,ix,ipairs,xx(l45),virvsene, xx(lpol),qpsander,.true.) call timer_stop(TIME_EWALD) endif call timer_start(TIME_QMMM) #ifdef PIMD !======================================================== ! PIMD QMMM !======================================================== call rst_mm_link_pair_crd(x) iqmatoms = qmmm_nml%iqmatoms iqmatomicnumbers = qmmm_struct%iqm_atomic_numbers atommask = qmmm_struct%atom_mask atomtype = qmmm_struct%qm_atom_type mm_link_mask(1:natom) = qmmm_struct%mm_link_mask(1:natom) if( .not.allocated( link_pairs ) ) then allocate( link_pairs(2,qmmm_struct%nlink), stat=alloc_error ) REQUIRE( alloc_error==0) end if link_pairs(1:2,1:qmmm_struct%nlink) = qmmm_struct%link_pairs(1:2,1:qmmm_struct%nlink) if( qmmm_nml%qm_pme ) empot = qmewald%empot qmmm_struct%nquant = qmmm_struct%nquant/ncopy qmmm_struct%nlink = qmmm_struct%nlink/ncopy qmmm_struct%nquant_nlink = qmmm_struct%nquant_nlink/ncopy do j=1,qmmm_struct%nquant_nlink qmmm_nml%iqmatoms(j) = iqmatoms( ncopy * j ) / ncopy qmmm_struct%iqm_atomic_numbers(j) = iqmatomicnumbers(ncopy*j) qmmm_struct%qm_atom_type(j) = atomtype( ncopy*j ) end do do j=1,natomCL qmmm_struct%atom_mask(j) = atommask(ncopy * j) qmmm_struct%mm_link_mask(j) = mm_link_mask( ncopy * j ) end do if(qm2_struct%calc_mchg_scf .or. qmmm_nml%qm_ewald > 0 .or. qmmm_nml%printcharges) & qm2_struct%scf_mchg = 0.0 nrg_scf = 0.0 do i=1,ncopy if( (ncopy.ge.numtasks.and.mod(i-1,numtasks)==mytaskid) .or. & (ncopy.lt.numtasks.and.mod(mytaskid,ncopy)==i-1) ) then if( .not.qmmm_pimd_first_call ) then qm2_struct%den_matrix(:) = pimd_den_matrix(:,i) qm2_struct%old_den_matrix(:) = pimd_old_den_matrix(:,i) qm2_struct%fock_matrix(:) = pimd_fock_matrix(:,i) end if call dispatch3( natomCL, i, ncopy, 1,1, f, frep ) call dispatch3( natomCL, i, ncopy, 1,1, x, xrep ) call dispatch ( natomCL, i, ncopy, 1,1, & qmmm_struct%scaled_mm_charges, mmchg ) if( qmmm_nml%qm_pme ) call dispatch( qmmm_struct%nquant_nlink, i, ncopy, & 1,1, empot, qmewald%empot ) if( qmmm_nml%qmgb == 2 ) then call dispatch( natomCL, i, ncopy,1,1,reff, reff_rep ) call dispatch( natomCL, i, ncopy,1,1,onereff, onereff_rep ) end if escf = 0.0 frep = 0.0 qmchg = 0.0 mmchg = mmchg * ncopy call qm_mm( xrep, natomCL, mmchg, frep, escf, periodic, a,b,c, & alpha,beta,gamma, reff_rep, onereff_rep, & intdiel, extdiel, Arad, cut, qmchg,ih(m06) ) escf = escf / ncopy frep = frep / ncopy qmchg= qmchg/ncopy nrg_scf( i ) = escf call dispatch3( natomCL, 1,1, i, ncopy, frep, f ); if( qm2_struct%calc_mchg_scf .or. qmmm_nml%qm_ewald > 0 ) then call dispatch( qmmm_struct%nquant_nlink,1,1,i,ncopy,qmchg, & qm2_struct%scf_mchg ) end if if( .not.allocated(pimd_den_matrix ) ) & allocate( pimd_den_matrix(qm2_struct%matsize, ncopy) ) if( .not.allocated(pimd_fock_matrix) ) & allocate( pimd_fock_matrix(qm2_struct%matsize, ncopy ) ) if( .not.allocated(pimd_old_den_matrix ) ) & allocate( pimd_old_den_matrix(qm2_struct%matsize, ncopy) ) pimd_den_matrix(:,i) = qm2_struct%den_matrix(:) pimd_fock_matrix(:,i) = qm2_struct%fock_matrix(:) pimd_old_den_matrix(:,i) = qm2_struct%old_den_matrix(:) ene(25) = ene(25) +escf end if end do qmmm_pimd_first_call = .false. qmmm_struct%nquant = qmmm_struct%nquant*ncopy qmmm_struct%nlink = qmmm_struct%nlink*ncopy qmmm_struct%nquant_nlink = qmmm_struct%nquant_nlink*ncopy qmmm_nml%iqmatoms = iqmatoms qmmm_struct%iqm_atomic_numbers = iqmatomicnumbers qmmm_struct%atom_mask = atommask qmmm_struct%mm_link_mask = mm_link_mask qmmm_struct%qm_atom_type = atomtype qmmm_struct%link_pairs(1:2,1:qmmm_struct%nlink) = link_pairs(1:2,1:qmmm_struct%nlink) if( qmmm_nml%qm_pme ) qmewald%empot = empot # ifdef MPI if( qm2_struct%calc_mchg_scf .or. qmmm_nml%qm_ewald > 0 .or. qmmm_nml%printcharges) then if(.not.qmmm_mpi%commqmmm_master) qm2_struct%scf_mchg = 0.0 call mpi_allreduce( qm2_struct%scf_mchg, qmchg, & qmmm_struct%nquant_nlink, MPI_DOUBLE_PRECISION, & mpi_sum, commsander, ierr ) qm2_struct%scf_mchg = qmchg(1:qmmm_struct%nquant_nlink) end if # endif !======================================================== ! END PIMD QMMM !======================================================== #else # ifndef LES !======================================================== ! REGULAR QMMM !======================================================== if(qmmm_nml%idc>0)then call qm_div(x, ix, f, mmoffset, escf, 1, mmoffset3) else call qm_mm(x, natom,qmmm_struct%scaled_mm_charges, & f,escf,periodic,a,b,c,alpha,beta,gamma,reff,onereff, & intdiel,extdiel,Arad, cut,qm2_struct%scf_mchg,ih(m06)) endif ene(25) = escf # endif #endif /* ifdef PIMD */ !======================================================== ! END REGULAR QMMM !======================================================== call timer_stop(TIME_QMMM) end if !if(qmmm_nml%ifqnt) !--------------------------------------------------------------- !END qm/mm contributions !--------------------------------------------------------------- ! ---------------------------------------------------------------- ! Calculate the non-bonded contributions ! ---------------------------------------------------------------- call timer_start(TIME_NONBON) call timer_start(TIME_EEXIPS) if( ips > 0 ) then call eexips(evdwex,eelex,istart,iend, ntb, & ix(i04),ix(i08),ix(i10),xx(l15),cn1,cn2,f,x,virex) vir(1) = vir(1)+0.5D0*(virips+virex(1,1)) vir(2) = vir(2)+0.5D0*(virips+virex(2,2)) vir(3) = vir(3)+0.5D0*(virips+virex(3,3)) endif call timer_stop(TIME_EEXIPS) if( igb == 0 .and. iyammp == 0 ) then ! (for GB: do all nonbondeds together below) call timer_start(TIME_EWALD) if ( iamoeba == 1 )then call AM_NonBond_eval(natom,x,f,vir,xx,ipairs, & evdw,eelt,epolar,& enb14,ee14,diprms,dipiter,qpsander) else if ( induced == 1 )then call handle_induced(x,natom,ix(i04),ix(i06),ntypes, & xx(l15),cn1,cn2,asol,bsol, & eelt,epolar,f,xx,ix,ipairs,xx(lpol), & xx(l45),virvsene,ix(i02),ibgwat,nres, & aveper,aveind,avetot,emtot,diprms,dipiter,dipole_temp,dt, & scee,scnb,ntb,qpsander) else call ewald_force(x,natom,ix(i04),ix(i06),ntypes, & xx(l15),cn1,cn2,asol,bsol,eelt,epolar, & f,xx,ix,ipairs,xx(l45),virvsene, xx(lpol),qpsander,.false.) end if end if ! iamoeba == 1 call timer_stop(TIME_EWALD) # ifdef MPI if(mytaskid == 0)then # endif ene(2) = evdw ene(3) = eelt ene(4) = ehb # ifdef MPI else ! energies have already been reduced to the master ! node in ewald_force, so here we zero out elements ! on non-master nodes: ene(2) = 0.d0 ene(3) = 0.d0 ene(4) = 0.d0 end if # endif if( ips > 0 )then ene(2)=ene(2)+evdwex ene(3)=ene(3)+eelex endif end if ! ( igb == 0 .and. iyammp == 0 ) call timer_stop(TIME_NONBON) ! ---------------------------------------------------------------- ! Calculate the other contributions ! ---------------------------------------------------------------- ! +---------------------------------------------------------------+ ! | Bonds with H | ! +---------------------------------------------------------------+ call timer_start(TIME_BOND) if( ntf < 2 ) then ebdev = 0.d0 if(qpsander)then call bond(nb2h,ib2h,jb2h,icb2h,x,xx,ix,f,ene(6),qpsander) else call bond(nbonh,ix(iibh),ix(ijbh),ix(iicbh),x,xx,ix,f,ene(6),qpsander) #ifdef LES if(rem == 2) then ene(24) = ene(24) + elesb endif #endif endif end if ! +---------------------------------------------------------------+ ! | Bonds without H | ! +---------------------------------------------------------------+ if( ntf < 3 ) then if(qpsander)then call bond(nb2a,ib2a,jb2a,icb2a,x,xx,ix,f,ene(7),qpsander) else call bond(nbona+nbper,ix(iiba),ix(ijba),ix(iicba),x,xx,ix,f,ene(7), & qpsander) #ifdef LES if(rem == 2) then ene(24) = ene(24) + elesb endif #endif endif if (nbonh+nbona > 0) ebdev = sqrt( ebdev/(nbonh+nbona) ) end if ! +---------------------------------------------------------------+ ! | Angles with H | ! +---------------------------------------------------------------+ if( ntf < 4 ) then eadev = 0.d0 if(qpsander)then call angl(nt2h,it2h,jt2h,kt2h,ict2h,x,xx,ix,f, & ene(8),qpsander) else call angl(ntheth,ix(i24),ix(i26),ix(i28),ix(i30),x,xx,ix,f,ene(8), & qpsander) #ifdef LES if(rem == 2) then ene(24) = ene(24) + elesa endif #endif endif end if ! +---------------------------------------------------------------+ ! | Angles without H | ! +---------------------------------------------------------------+ if( ntf < 5 ) then if(qpsander)then call angl(nt2a,it2a,jt2a,kt2a,ict2a,x,xx,ix,f, & ene(9),qpsander) else call angl(ntheta+ngper,ix(i32),ix(i34),ix(i36),ix(i38),x,xx,ix,f, & ene(9),qpsander) #ifdef LES if(rem == 2) then ene(24) = ene(24) + elesa endif #endif endif if (ntheth+ntheta > 0) eadev = 57.296*sqrt( eadev/(ntheth+ntheta) ) end if ! +---------------------------------------------------------------+ ! | Dihedrals with H | ! +---------------------------------------------------------------+ if( ntf < 6 ) then if(qpsander)then call ephi(nphi2h, & ip2h,jp2h,kp2h,lp2h,icp2h,xx(l15), & ix(i04),x,xx,ix,f,dvdl,ene(10),ene(11),ene(12),xx(l190),qpsander) else call ephi(nphih,ix(i40),ix(i42),ix(i44),ix(i46),ix(i48), & xx(l15),ix(i04),x,xx,ix,f,dvdl,ene(10),ene(11),ene(12),xx(l190), & qpsander) #ifdef LES if(rem == 2) then ene(24) = ene(24) + elesd endif #endif endif end if ! +---------------------------------------------------------------+ ! | Dihedrals without H | ! +---------------------------------------------------------------+ if( ntf < 7 ) then if(qpsander)then call ephi(nphi2a, & ip2a,jp2a,kp2a,lp2a,icp2a,xx(l15), & ix(i04),x,xx,ix,f,dvdl,ene(13),ene(14),ene(15),xx(l190), & qpsander) else call ephi(nphia+ndper,ix(i50),ix(i52),ix(i54),ix(i56),ix(i58), & xx(l15),ix(i04),x,xx,ix,f,dvdl,ene(13),ene(14),ene(15),xx(l190), & qpsander) #ifdef LES if(rem == 2) then ene(24) = ene(24) + elesd endif #endif endif end if if ( iamoeba == 1 ) call AM_VAL_eval(x,f,vir,ene(6),ene(8),ene(10)) call timer_stop(TIME_BOND) ! --- calculate the position constraint energy --- if(natc > 0 .and. ntr==1) then ! ntr=1 (positional restraints) call xconst(natc,entr,ix(icnstrgp),x,f,xx(lcrdr),xx(l60),natom, & nmyatm,spmylist,qpsander ) ene(20) = entr end if if ( itgtmd==1 .and. (nattgtfit > 0 .or. nattgtrms > 0) ) then ! Calculate rmsd for targeted md (or minimization) if requested. ! All nodes do rms fit, could just be master then broadcast. ! All nodes need all coordinates for this. call rmsfit(xx(lcrdr),x,xx(lmass),ix(itgtfitgp), & ix(itgtrmsgp),rmsdvalue,rmsok) if (.not.rmsok) then if (master) write (6,*) 'Fatal Error calculating RMSD !' call mexit(6, 1) end if call xtgtmd(entr,ix(itgtrmsgp),x,f,xx(lcrdr),xx(lmass)) ene(20) = entr end if if(ifcap > 0) then call capwat(natom,x,f,ecap) ene(20) = ene(20) + ecap end if ! (this seems very weird: we have already done an allreduce on molvir ! in ewald_force(); this just collects it on processor 0 (with zeroes ! on all slave nodes), then later does an allreduce...) if( mytaskid == 0 .and. iamoeba == 0 ) then vir(1) = vir(1)+0.5d0*molvir(1,1) vir(2) = vir(2)+0.5d0*molvir(2,2) vir(3) = vir(3)+0.5d0*molvir(3,3) end if if( igb == 0 .and. iyammp == 0 ) then ener(21) = virvsene ener(22) = diprms ener(23) = dipiter ener(24) = dipole_temp end if ! ---- get the noesy volume penalty energy: ------ enoe = 0.d0 if( iredir(4) /= 0 ) then call timer_start(TIME_NOE) call noecalc(x,f,xx,ix) call timer_stop(TIME_NOE) end if ene(22) = enoe ! -- when igb!=0 and igb!=10, all nonbonds are done in routine egb: esurf = 0.d0 if( igb /= 0 .and. igb /= 10) then call timer_start(TIME_EGB) call egb( x,f,rborn,fs,reff,onereff,xx(l15),ix(i04),ix(i06), & ix(i08),ix(i10),xx(l190), & cut,ntypes,natom,natbel,epol,eelt,evdw, & esurf,dvdl,xx(l165),ix(i82),xx(l170),xx(l175),xx(l180), & xx(l185), xx(l186),xx(l187),xx(l188),xx(l189),ncopy ) ene(2) = evdw ene(3) = eelt ene(4) = epol ene(23) = esurf ene(21) = dvdl call timer_stop(TIME_EGB) #ifdef LES if(rem == 2) then ene(24) = ene(24) + elesp endif #endif end if ! ( igb /= 0 .and. igb /= 10) ! -- when igb==10, all nonbonds are done in routine pb_force, and ! all nonpolar interactions are done in np_force: #ifdef MPI if(mytaskid == 0)then #endif if( igb == 10 ) then call timer_start(TIME_PBFORCE) call pb_force(natom,nres,ntypes,ix(i02),ix(i04),ix(i06),ix(i10), & cn1,cn2,xx(l15),x,f,evdw,eelt,epol) if ( pbgrid ) pbgrid = .false. if ( pbinit ) pbinit = .false. ene(2) = evdw ene(3) = eelt ene(4) = epol call timer_stop(TIME_PBFORCE) call timer_start(TIME_NPFORCE) esurf = 0.0d0; edisp = 0.0d0 if ( ifcap == 0 .and. npopt /= 0 ) & call np_force(natom,nres,ntypes,ix(i02),ix(i04),ix(i06),& cn1,cn2,x,f,esurf,edisp) if ( pbprint ) pbprint = .false. ene(23) = esurf ene(26) = edisp call timer_stop(TIME_NPFORCE) end if ! ( igb == 10 ) #ifdef MPI end if #endif #ifdef MPI call timer_barrier( commsander ) call timer_start(TIME_COLLFRC) ! add force, ene, vir, copies from all nodes ! also add up newbalance for nonperiodic. # ifdef PSANDER call get_stack(tmpene,50,routine) call get_stack(tmpfrc,maxf_sndsiz*ntskfrcsnd,routine) call get_stack(tmpf2, maxf_rcvsiz*ntskfrcrcv,routine) if(.not. rstack_ok)then deallocate(r_stack) allocate(r_stack(1:lastrst),stat=alloc_ier) call reassign_rstack(routine) endif REQUIRE(rstack_ok) call spatial_fdist(f,ene,vir, & r_stack(tmpene),r_stack(tmpfrc),r_stack(tmpf2)) call free_stack(tmpf2,routine) call free_stack(tmpfrc,routine) call free_stack(tmpene,routine) # else call fdist(f,xx(lfrctmp),ene,vir,newbalance) # endif call timer_stop(TIME_COLLFRC) #endif ! ---- at this point, the parallel part of the force calculation is ! finished, and the forces have been distributed to their needed ! locations. All forces below here are computed redundantly on ! all processors, and added into the force vector. Hence, below ! is the place to put any component of the force calculation that ! has not (yet) been parallelized. ! Calculate the NMR restraint energy contributions, if requested. if (nmropt > 0) then call nmrcal(x,f,ih(m04),ih(m02),ix(i02),xx(lwinv),enmr,devdis, & devang,devtor,temp0,tautp,cut,ntb,xx(lnmr01), & ix(inmr02),xx(l95),31,6,rk,tk,pk,cn1, & cn2,asol,bsol,xx(l15),numbnd,numang,nptra-nimprp, & nimprp,nttyp,nphb,natom,natom,ntypes,nres, & rad,wel,radhb,welhb,rwell,isftrp,tgtrmsd,temp0les,-1,'CALC') end if eshf = 0.d0 epcshf = 0.d0 ealign = 0.d0 if (iredir(5) /= 0) call cshf(natom,x,f) if (iredir(7) /= 0) call pcshift(natom,x,f) if (iredir(8) /= 0) call align1(natom,x,f,xx(lmass)) enoe = ene(22) ! ----- CALCULATE TOTAL ENERGY AND GROUP THE COMPONENTS ----- #if defined(LES) || defined(PIMD) #else if( igb == 0 ) then ene(11) = enb14 ene(14) = 0.d0 ene(12) = ee14 ene(15) = 0.d0 end if #endif do m = 2,15 ene(1) = ene(1) + ene(m) end do ene(1) = ene(1) + epolar + ene(23) + ene(25) + ene(26) ene(5) = ene(6)+ene(7) ene(6) = ene(8)+ene(9) ene(7) = ene(10)+ene(13) ene(8) = ene(11)+ene(14) ene(9) = ene(12)+ene(15) ene(10) = ene(17)+ene(20)+eshf+epcshf+enoe+enmr(1)+enmr(2)+enmr(3)+ealign ene(1) = ene(1)+ene(10) ! Here is a summary of how the ene array is used. For parallel runs, ! these values get summed then rebroadcast to all nodes (via ! mpi_allreduce). ! ene(1): total energy ! ene(2): van der Waals ! ene(3): electrostatic energy ! ene(4): 10-12 (hb) energy, or GB energy when igb.gt.0 ! ene(5): bond energy ! ene(6): angle energy ! ene(7): torsion angle energy ! ene(8): 1-4 nonbonds ! ene(9): 1-4 electrostatics ! ene(10): constraint energy ! ene(11-19): used a scratch, but not needed further below ! ene(20): position constraint energy + cap energy ! ene(21): charging free energy result ! ene(22): noe volume penalty ! ene(23): surface-area dependent energy, or cavity energy ! ene(24): potential energy for a subset of atoms ! ene(25): SCF Energy when doing QMMM ! ene(26): implicit solvation dispersion energy ! ----- TRANSFER THE ENERGIES TO THE ARRAY ENER, USED IN PRNTMD ----- ener(1:10) = ene(1:10) ener(11) = epolar ener(12) = aveper ener(13) = aveind ener(14) = avetot ener(15) = ene(23) ener(17) = ene(21) ener(18) = ene(24) ener(25) = ene(25) ener(28) = ene(26) ! ----- IF BELLY IS ON THEN SET THE BELLY ATOM FORCES TO ZERO ----- ! FROZEN_FORCES if (belly .and. master) call prt_frozen_frcs(natom,ix(ibellygp),f) ! FROZEN_FORCES if (belly) call bellyf(natom,ix(ibellygp),f) ! +---------------------------------------------------------------+ ! | Interface to EVB | ! +---------------------------------------------------------------+ #if defined(MPI) if( ievb /= 0 ) then if( nmorse > 0 ) then if( .not. morsify_initialized ) call morsify_init ( ix ) call morsify ( x, ener(5), ener(1), f ) endif #if defined(PIMD) nrg_bead(:) = nrg_bnd(:)+nrg_ang(:)+nrg_dih(:)+nrg_vdw(:)+nrg_ele(:) if( evb_pimd_debug ) then write(6,'(2(/))') write(6,'( A,(8F10.5, 2X))') 'nrg_bnd = ', nrg_bnd(:), sum( nrg_bnd ) write(6,'( A,(8F10.5, 2X))') 'nrg_ang = ', nrg_ang(:), sum( nrg_ang ) write(6,'( A,(8F10.5, 2X))') 'nrg_dih = ', nrg_dih(:), sum( nrg_dih ) write(6,'( A,(8F10.5, 2X))') 'nrg_vdw = ', nrg_vdw(:), sum( nrg_vdw ) write(6,'( A,(8F10.5, 2X))') 'nrg_ele = ', nrg_ele(:), sum( nrg_ele ) write(6,'( A,(8F10.5, 2X))') 'total_E = ', nrg_bead(:), sum( nrg_bead ) fpimd(:,:,:) = reshape( f(1:3*natom), (/ 3, nbead, natomCL /) ) do n = 1, nbead write(6,'(/)') write(6,'(A/, ( 6F12.7) )' ) '<< F = ', fpimd(:,n,:) enddo endif piq(:,:,:) = reshape( x(1:3*natom), (/ 3, nbead, natomCL /) ) call mpi_allgather ( nrg_bead, nbead, MPI_DOUBLE_PRECISION, pie, nbead & , MPI_DOUBLE_PRECISION, commworld, ierr ) call mpi_allgather ( f, natm3*nbead, MPI_DOUBLE_PRECISION, pif, & natm3*nbead, MPI_DOUBLE_PRECISION, commworld, ierr ) call mpi_barrier ( commworld, ierr ) evb_nrg_sum = 0.0d0 do n = evb_begin(worldrank+1), evb_end(worldrank+1) call evb_ntrfc( pif(:,n,:,:), piq(:,n,:), pie(n,:), natm3 ) call mpi_send ( evb_frc% evb_f(:), natm3, MPI_DOUBLE_PRECISION, 0, & n, commworld, ierr ) call mpi_send ( evb_frc% evb_nrg, 1, MPI_DOUBLE_PRECISION, 0, & n+nbead, commworld, ierr ) enddo if( worldrank == 0 ) then do n = 1, nbead call mpi_recv( f_tmp, natm3, MPI_DOUBLE_PRECISION, MPI_ANY_SOURCE, & n, commworld, status, ierr ) frc_tag = status( MPI_TAG ) fpimd(:,frc_tag,:) = f_tmp(:,:,1) call mpi_recv( evb_nrg, 1, MPI_DOUBLE_PRECISION, MPI_ANY_SOURCE, & n+nbead, commworld, status, ierr ) nrg_tag = status( MPI_TAG ) - nbead evb_bead(nrg_tag) = evb_nrg enddo evb_nrg_sum = sum( evb_bead(:) ) endif call mpi_barrier ( commworld, ierr ) if( worldrank == 0 ) then f(1:natm3*nbead) = reshape( fpimd(:,:,:), (/ natm3*nbead /) ) endif if( numtasks > 0 ) then call mpi_bcast ( f, natm3*nbead, MPI_DOUBLE_PRECISION, 0, & commworld, ierr ) endif #else /* PIMD */ root = 0 !KFW call mpi_bcast( x, natm3, MPI_DOUBLE_PRECISION, root, commworld, ierr ) if( master ) then call mpi_allgather ( f, natm3, MPI_DOUBLE_PRECISION, xf, natm3, & MPI_DOUBLE_PRECISION, commmaster, ierr ) call mpi_allgather ( ener(1), 1, MPI_DOUBLE_PRECISION, xnrg, 1, & MPI_DOUBLE_PRECISION, commmaster, ierr ) xq(:) = x(1:natm3) call evb_ntrfc( xf, xq, xnrg, natm3 ) f(1:natm3) = evb_frc% evb_f(:) endif if( numtasks > 0 ) then call mpi_bcast( f, natm3, MPI_DOUBLE_PRECISION, root, commworld, ierr ) endif #endif /* PIMD */ endif #endif /*MPI*/ #ifdef PIMD if( ineb > 0 ) then nrg_all = nrg_bnd+nrg_ang+nrg_dih+nrg_vdw+nrg_ele+nrg_egb+nrg_scf ener(1:26) = ener(1:26) * ncopy # ifdef MPI call mpi_allreduce( nrg_all, nrg_tmp, ncopy, MPI_DOUBLE_PRECISION, & mpi_sum, commsander, ierr ) nrg_all = nrg_tmp # endif nrg_all = nrg_all * ncopy f(1:natom*3) = f(1:natom*3) * ncopy call pimd_neb_forces( x, f, spring_energy, xx(lvel) ) else call pimd_spring_force ( x, f, spring_energy, Epot_deriv ) etot = ener(1) + equal_part + Epot_deriv end if #endif /*PIMD*/ #ifdef MPI if( ievb /= 0 .and. worldrank == 0 ) then call out_evb ( etot ) evb_step = evb_step + 1 endif #endif call timer_stop(TIME_FORCE) call trace_exit( 'force' ) return end subroutine force ! ------------------------------------------------------------------- ! Allocate storage space for xq, xf, xnrg for EVB ! ------------------------------------------------------------------- subroutine evb_amber_alloc use evb_parm, only: nevb, natm3 use evb_amber, only: xnrg, xq, xf implicit none integer :: alloc_error allocate( xnrg(nevb), stat = alloc_error ) REQUIRE( alloc_error == 0 ) allocate( xq(natm3), stat = alloc_error ) REQUIRE( alloc_error == 0 ) allocate( xf(natm3,nevb), stat = alloc_error ) REQUIRE( alloc_error == 0 ) end subroutine evb_amber_alloc ! ------------------------------------------------------------------- ! Deallocate storage space for xnrg, xq, xf for EVB ! ------------------------------------------------------------------- subroutine evb_amber_dealloc use evb_amber, only: xnrg, xq, xf implicit none integer :: dealloc_error deallocate( xnrg, xq, xf, stat = dealloc_error ) REQUIRE( dealloc_error == 0 ) end subroutine evb_amber_dealloc