#include "copyright.h" #include "../include/dprec.fh" #include "../include/assert.fh" #include "ncsu-config.h" !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+ The Molecular Dynamics/NMR Refinement/Modeling Module of the AMBER !----------------------------------------------------------------------- ! --- SANDER --- !----------------------------------------------------------------------- subroutine sander() use state #ifndef DISABLE_NCSU use ncsu_sander_hooks, only : & ncsu_on_sander_init => on_sander_init, & ncsu_on_sander_exit => on_sander_exit #endif /* DISABLE_NCSU */ use lmod_driver use constants, only : INV_AMBER_ELECTROSTATIC ! The main qmmm_struct contains all the QMMM variables and arrays use qmmm_module, only : qmmm_nml, qmmm_struct, deallocate_qmmm, qmmm_mpi, & #ifdef MPI qmmm_mpi_setup, & #endif qm2_struct, qmewald, qm_gb, qmmm_vsolv, qm2_params use qmmm_vsolv_module, only: qmmm_vsolv_store_parameters, new use qm2_extern_module, only: qm2_extern_finalize use sebomd_module, only : sebomd_obj, & sebomd_open_files, sebomd_close_files, & #ifdef MPI sebomd_bcast_obj, & #endif sebomd_setup use sebomd_arrays, only : init_sebomd_arrays, cleanup_sebomd_arrays use genborn use decomp, only : allocate_int_decomp, allocate_real_decomp, & deallocate_int_decomp, deallocate_real_decomp, & #ifdef MPI synchronize_dec, build_dec_mask, decmask, & #endif nat, nrs, jgroup, indx use fastwt use relax_mat use nmr, only: nmrrad, impnum use ew_recip, only: deallocate_m1m2m3,first_pme use parms use molecule, only : mol_info, n_iwrap_mask_atoms, iwrap_mask_atoms, & allocate_molecule, deallocate_molecule use nblist, only:cutoffnb,skinnb,nblist_allocate,nblist_deallocate, & nblist_allreal,nblist_allint, num_calls_nblist, first_list_flag use stack use amoeba_runmd, only : AM_RUNMD_get_coords,AM_RUNMD use amoeba_mdin, only : beeman_integrator,iamoeba,am_nbead use amoeba_interface, only : AMOEBA_deallocate,AMOEBA_readparm #ifdef RISMSANDER use sander_rism_interface, only: rismprm,rism_setparam, rism_init,rism_finalize #endif #ifdef PUPIL_SUPPORT use pupildata #endif /* PUPIL */ #ifdef APBS use apbs #endif /* APBS */ #ifdef _XRAY use xray_interface_module, only: xray_init, xray_read_parm, xray_fini #endif #ifdef MPI /* SOFT CORE */ use softcore, only: setup_sc, cleanup_sc, ifsc, extra_atoms, sc_sync_x, & summarize_ti_changes, sc_check_perturbed_molecules, ti_check_neutral, tishake use mbar, only: setup_mbar, cleanup_mbar, ifmbar #endif ! for LIE calculations use linear_response, only: ilrt, setup_linear_response, & cleanup_linear_response #if defined(MPI) use evb_parm, only: xch_type # if defined(LES) use evb_pimd, only: evb_pimd_init, PE_slice, master_worldrank, jobs_per_node # endif ! REMD use remd, only : rem, mdloop, repnum, remd1d_setup, remd_exchange, & remd_cleanup, hremd_exchange, ph_remd_exchange, & multid_remd_setup, multid_remd_exchange use bintraj, only: setup_remd_indices #else # define rem 0 #endif /* MPI */ use pimd_vars, only: ipimd use neb_vars, only: ineb use trajenemod, only: trajene !RCW+MJW CHARMM SUPPORT use charmm_mod, only : charmm_active, charmm_deallocate_arrays, & charmm_filter_out_qm_atoms use ff11_mod, only : cmap_active, deallocate_cmap_arrays use memory_module, only: x, ix, ih, memory_init ! Self-Guided molecular/Langevin Dynamics (SGLD) use sgld, only : isgld,psgld use nbips, only: ipssys,ips use crg_reloc, only: ifcr, cr_backup_charge, cr_cleanup, cr_allocate, & cr_read_input, cr_check_input, cr_print_info use emap,only: temap,pemap,qemap use file_io_dat use constantph, only : cnstph_finalize use barostats, only : mcbar_setup !AMD use amd_mod !scaledMD use scaledMD_mod use abfqmmm_module ! lam81 implicit none logical belly, erstop integer ier,ifind,jn,ncalls,xmin_iter character(len=4) itest logical ok logical newstyle # include "../include/memory.h" # include "nmr.h" # include "box.h" # include "../include/md.h" # include "extra.h" # include "tgtmd.h" # include "multitmd.h" # include "les.h" # include "parallel.h" #ifdef MPI ! =========================== AMBER/MPI =========================== # ifdef MPI_DOUBLE_PRECISION # undef MPI_DOUBLE_PRECISION # endif include 'mpif.h' # ifdef CRAY_PVP # define MPI_DOUBLE_PRECISION MPI_REAL8 # endif # ifdef MPI_BUFFER_SIZE integer*4 mpibuf(mpi_buffer_size) # endif ! REMD: loop is the current exchange. runmd is called numexchg times. integer loop integer nrank, istat _REAL_ ener(30),vir(4) integer ierr integer partner ! ========================= END AMBER/MPI ========================= #endif # include "ew_pme_recip.h" # include "ew_frc.h" # include "ew_erfc_spline.h" # include "ew_parallel.h" # include "ew_mpole.h" # include "ew_cntrl.h" # include "def_time.h" type(state_rec) :: ene integer native,nr3,nr ! nmrcal vars _REAL_ f,enmr,devdis,devang,devtor,devplpt,devpln,devgendis,ag,bg,cg ! Updated 9/2007 by Matthew Seetin to enable plane-point and plane-plane ! restraints _REAL_ emtmd integer numphi,nhb ! runmin/trajene var _REAL_ carrms ! dipole moment stuff integer ngrp character(len=8) initial_date, setup_end_date, final_date character(len=10) initial_time, setup_end_time, final_time integer nstlim_total ! for final time printout _REAL_ time0, time1 integer idiff,i,j,istop,index,ierror,itemp integer, dimension(:), allocatable :: ipairs integer :: n_force_calls logical qsetup logical :: do_list_update=.false. _REAL_ :: box_center(3) ! lam81 #ifdef MPI_DEBUGGER integer, volatile :: release_debug ! So only the master master thread has release_debug = 0 release_debug = worldrank ! Lock us into an infinite loop while release_debug == 0 on any thread (only ! the master here). This allows you to connect a debugger to any running ! process without having to 'race' program execution. A debugger MUST be ! attached to the master thread (typically the thread with the lowest PID), ! and have release_debug set to NOT 0 (e.g., via "set release_debug=1"). do if (release_debug .ne. 0) exit end do ! Prevent any other threads from progressing past this point until all ! threads you want to watch with a debugger are watched and all those ! threads are continued. call mpi_barrier(mpi_comm_world, ierr) #endif ! ---- HERE BEGIN THE EXECUTABLE STATEMENTS ---- ! Initialize the cpu timer. Needed for machines where returned cpu times ! are relative. call date_and_time( initial_date, initial_time ) call wallclock( time0 ) call init_timers() ! Initialize the printing of ongoing time and performance summaries. ! call print_ongoing_time_summary(0,0,0.0d0,7) call print_ongoing_time_summary(0,0,0.0d0,7,time0) ! BPR - original location of PUPIL interface. I moved it further down ! because, if it's here, it can't print stuff; write(6,...) statements ! assume mdread1() has already been invoked. However, moving this down ! may break other things. ! ==== Flag to tell list builder to print size of list on first call ======= first_list_flag = .true. ! ==== Flag to tell recip space routines to allocate on first call ======= first_pme = .true. ! ==== Initialise first_call flags for QMMM ==== qmmm_struct%qm_mm_first_call = .true. qmmm_struct%fock_first_call = .true. qmmm_struct%fock2_2atm_first_call = .true. qmmm_struct%qm2_allocate_e_repul_first_call = .true. qmmm_struct%qm2_calc_rij_eqns_first_call = .true. qmmm_struct%qm2_scf_first_call = .true. qmmm_struct%zero_link_charges_first_call = .true. qmmm_struct%adj_mm_link_pair_crd_first_call = .true. qmmm_struct%num_qmmm_calls = 0 #ifdef MPI ! Parallel initialization (setup is done in multisander.F90). ! Make PE 0 the master master = mytaskid == 0 master_master = masterrank == 0 if ( master .and. numtasks > MPI_MAX_PROCESSORS ) then write(0, '(a,i4,a,i4)') & 'Error: the number of processors must not be greater than ', & MPI_MAX_PROCESSORS, ', but is ', numtasks call mexit(6,1) end if # ifdef MPI_BUFFER_SIZE call mpi_buffer_attach(mpibuf, mpi_buffer_size*4, ierr) # endif #else /* not MPI follows */ ! In the single-threaded version, the one process is master master = .true. #endif /* MPI */ erstop = .false. qsetup = .true. ! --- generic packing scheme --- nwdvar = 1 native = 32 #ifdef ISTAR2 ! --- Int*2 packing scheme --- nwdvar = 2 #endif /*ISTAR2*/ numpk = nwdvar nbit = native/numpk ! ----- Only the master node (only node when single-process) ! performs the initial setup and reading/writing ----- call timer_start(TIME_TOTAL) call abfqmmm_init_param() ! lam81 do while ( (abfqmmm_param%qmstep <= abfqmmm_param%maxqmstep) & .or. (abfqmmm_param%maxqmstep == 0 .and. abfqmmm_param%system == 2) ) ! lam81 masterwork: if (master) then if (abfqmmm_param%abfqmmm == 0) then ! lam81 ! ---- first, initial reads to determine memory sizes: call mdread1() call amopen(8,parm,'O','F','R') call rdparm1(8) if (mtmd /= 'mtmd' .or. itgtmd == 2) call mtmdlx(natom) ! --- now, we can allocate memory: call locmem() write(6,'(/,a,5x,a)') '|','Memory Use Allocated' write(6,'(a,5x,a,i14)') '|', 'Real ', lastr write(6,'(a,5x,a,i14)') '|', 'Hollerith ', lasth write(6,'(a,5x,a,i14)') '|', 'Integer ', lasti write(6,'(a,5x,a,i14)') '|', 'Max Pairs ', lastpr ! --- dynamic memory allocation: ! GMS: ! Allocate space for module molecule ! in the master node mol_info%natom = natom mol_info%nres = nres call allocate_molecule() ! Allocate all global arrays allocate( x(lastr), ix(lasti), ipairs(lastpr), ih(lasth), stat = ier ) REQUIRE( ier == 0 ) ix(1:lasti) = 0 ! This sets up pointer arrays in MEMORY_MODULE to match array-offsets into ! the shared X, IX, and IH arrays. Eventually, LOCMEM code should be ! merged with MEMORY_MODULE to allocate individual allocatable arrays, but ! that will also require updating the MPI code to handle individual ! arrays. call memory_init() ! Allocate the parm arrays call allocate_parms() if ((igb /= 0 .and. igb /= 10 .and. ipb == 0) & .or.hybridgb>0.or.icnstph.gt.1) & call allocate_gb( natom, ncopy ) if( idecomp > 0 ) then #ifdef MPI if (ifsc > 0) then call synchronize_dec(natom, nres) else nat = natom nrs = nres end if #else nat = natom nrs = nres #endif call allocate_int_decomp(natom, nres) else call allocate_int_decomp(1, 1) end if write(6,'(a,5x,a,i14)' ) '|', 'nblistReal', nblist_allreal write(6,'(a,5x,a,i14)' ) '|', 'nblist Int', nblist_allint write(6,'(a,5x,a,i14,a)') '|', ' Total ', & (8*(lastr+lastrst+nblist_allreal) & + 4*(lasth+lasti+lastpr+lastist+nblist_allint))/1024, & ' kbytes' ! --- finish reading the prmtop file and other user input: call rdparm2(x,ix,ih,ipairs,8) call AMOEBA_readparm(8,ntf,ntc,natom,x(lmass))! ntf,ntc get reset if amoeba prmtop #ifdef _XRAY call xray_read_parm(8,6) #endif end if ! lam81 if (qmmm_nml%ifqnt .or. abfqmmm_param%abfqmmm == 1) then ! lam81 if(abfqmmm_param%abfqmmm == 0) then ! lam81 call sebomd_setup call read_qmmm_nm_and_alloc(igb, ih, ix, x, cut, use_pme, ntb, 0) ! lam81 if (qmmm_nml%qmtheory%SEBOMD) then ! don't do QM/MM qmmm_nml%ifqnt= .false. sebomd_obj%do_sebomd = .true. end if end if ! lam81 if(qmmm_struct%abfqmmm == 1 .and. abfqmmm_param%abfqmmm == 0) then ! lam81 call abfqmmm_setup(natom,nres,ix(i02),ih(m04),ih(m02),x(lmass), & ! lam81 nbonh,nbona,ix(iibh),ix(ijbh),ix(iiba),ix(ijba)) ! lam81 nr=natom ! lam81 call AMOEBA_check_newstyle_inpcrd(inpcrd,newstyle) ! lam81 if (newstyle) then ! lam81 call AM_RUNMD_get_coords(natom,t,irest,ntb,x(lcrd),x(lvel)) ! lam81 else ! lam81 call getcor(nr,x(lcrd),x(lvel),x(lforce),ntx,box,irest,t,temp0,.FALSE.) ! lam81 end if ! lam81 abfqmmm_param%maxqmstep = nstlim ! lam81 end if ! lam81 if(abfqmmm_param%abfqmmm == 1) then ! lam81 if(abfqmmm_param%system == 1) then ! lam81 call abfqmmm_update_qmatoms(x(lcrd)) ! lam81 if(abfqmmm_param%ntwpdb < 0) then ! lam81 call abfqmmm_write_pdb(x(lcrd),ix(i70)) ! lam81 close(6) ! lam81 call mexit(6,1) ! lam81 end if ! lam81 end if ! lma81 call abfqmmm_select_system_qmatoms(natom) ! lam81 if(qmmm_nml%ifqnt) then ! lam81 call read_qmmm_nm_and_alloc(igb,ih,ix,x,cut,use_pme,ntb,abfqmmm_param%qmstep, & ! lam81 abfqmmm_param%isqm,abfqmmm_param%abfcharge) ! lam81 end if ! lam81 endif ! lam81 end if ! lam81 if (abfqmmm_param%qmstep == 1 .and. abfqmmm_param%system == 1) then ! lam81 call mdread2(x,ix,ih,ipairs) endif ! lam81 #if defined(RISMSANDER) call rism_setparam(mdin,& commsander,& natom,ntypes,x(L15:L15+natom-1),& x(LMASS:LMASS+natom-1),cn1,cn2,& ix(i04:i04+ntypes**2-1), ix(i06:i06+natom-1)) #endif /*RISMSANDER*/ if ( ifcr /= 0 ) then call cr_read_input(natom) call cr_check_input( ips ) !call cr_print_info(6) call cr_backup_charge( x(l15), natom ) end if ! --- alloc memory for decomp module that needs info from mdread2 if (idecomp == 1 .or. idecomp == 2) then call allocate_real_decomp(nrs) #ifdef MPI ! -- ti decomp if (ifsc > 0) then ! following lines don't really seem to make sense(?) ! partner = ieor(masterrank,1) ! if (nat == natom) then ! nrank = masterrank ! else ! nrank = partner ! end if call mpi_bcast(jgroup, nat, MPI_INTEGER, 0, commmaster, ierr) end if #endif else if( idecomp == 3 .or. idecomp == 4 ) then call allocate_real_decomp(npdec*npdec) end if ! ----- EVALUATE SOME CONSTANTS FROM MDREAD SETTINGS ----- nr = nrp nr3 = 3*nr belly = ibelly > 0 ! ========================= PUPIL INTERFACE ========================= #ifdef PUPIL_SUPPORT ! I moved the PUPIL interface down here so that write() statements work ! as advertised. BPR 9/7/09 ! Initialise the CORBA interface puperror = 0 call fixport() call inicorbaintfcmd(puperror) if (puperror .ne. 0) then write(6,*) 'Error creating PUPIL CORBA interface.' call mexit(6,1) end if pupactive = .true. write(6,*) 'PUPIL CORBA interface initialized.' ! Allocation of memory and initialization pupStep = 0 puperror = 0 allocate (qcell (12 ),stat=puperror) allocate (pupmask (natom ),stat=puperror) allocate (pupqlist(natom ),stat=puperror) allocate (pupatm (natom ),stat=puperror) allocate (pupchg (natom ),stat=puperror) allocate (qfpup (natom*3),stat=puperror) allocate (qcdata (natom*9),stat=puperror) allocate (keyMM (natom ),stat=puperror) allocate (pupres (nres ),stat=puperror) allocate (keyres (nres ),stat=puperror) if (puperror /= 0) then write(6,*) 'Error allocating PUPIL interface memory.' call mexit(6,1) end if ! Initialise the "atomic numbers" and "quantum forces" vectors pupqatoms = 0 iresPup = 1 pupres(1) = 1 do iPup=1,natom bs1 = (iPup-1)*3 call get_atomic_number_pupil(ih(iPup+m06-1),x(lmass+iPup-1),pupatm(iPup)) if (iresPup .lt. nres) then if (iPup .ge. ix(iresPup+i02)) then iresPup = iresPup + 1 pupres(iresPup) = iPup end if end if write (strAux,"(A4,'.',A4)") trim(ih(iresPup+m02-1)),adjustl(ih(iPup+m04-1)) keyres(iresPup) = trim(ih(iresPup+m02-1)) keyMM(iPup) = trim(strAux) ! Retrieve the initial charges pupchg(iPup) = x(L15+iPup-1) !write(6,*) 'Atom num.',iPup,'Label,Mass,Atomic Num.', keyMM(iPup),x(lmass+iPup-1),pupatm(iPup), 'Charge', pupchg(iPup) do jPup=1,3 qfpup(bs1+jPup) = 0.0d0 end do end do write(6,*) 'Got all atomic numbers.' ! Initialise the PUPIL cell do iPup=1,12 qcell(iPup) = 0.0d0 end do ! Submit the KeyMM particles and their respective atomic numbers to PUPIL puperror = 0 call putatomtypes(natom,puperror,pupatm,keyMM) if (puperror .ne. 0) then write(6,*) 'Error sending MM atom types to PUPIL.' call mexit(6,1) end if ! Submit the Residue Pointer vector to PUPIL write(6,"(a20,1x,i6,3x,a17,1x,i6)") 'Number of residues =', nres, 'Number of atoms =', natom !do iPup=1,nres ! write(6,*) 'Residue ',iPup,keyres(iPup),pupres(iPup) !end do puperror = 0 call putresiduetypes(nres,puperror,pupres,keyres) if (puperror .ne. 0) then write(6,*) 'Error sending MM residue types to PUPIL.' call mexit(6,1) end if write(6,*) 'Sent system data to PUPIL.' write(*,*) 'PUPIL structure initialized.' #endif ! ========================= PUPIL INTERFACE ========================= ! --- seed the random number generator --- ! DAN ROE: Note master node only here if (abfqmmm_param%qmstep == 1 .and. abfqmmm_param%system == 1) then ! lam81 #ifdef MPI if (rem == 0) then call amrset(ig) else ! carlos: set random seed different for different replicas ! but keep same seed for cpus in the same replica since ! we want data from diff # cpus to match call amrset(ig + (17 * nodeid)) ! nodeid is in md.h and is repnum - 1 end if #else call amrset(ig) #endif if (nbit < 32 .and. nr > 32767) then write(6, *) ' Too many atoms for 16 bit pairlist -' write(6, *) ' Recompile without ISTAR2' call mexit(6, 1) end if if (ntp > 0.and.iabs(ntb) /= 2) then write(6,*) 'Input of NTP/NTB inconsistent' call mexit(6, 1) end if end if ! lam81 ! ----- READ COORDINATES AND VELOCITIES ----- if (abfqmmm_param%qmstep == 1 .and. abfqmmm_param%system == 1) then ! lam81 <<< QMSTEP=1 BLOCK >>> call timer_start(TIME_RDCRD) #ifdef LES call getcor(nr,x(lcrd),x(lvel),x(lforce),ntx,box,irest,t,temp0les,.TRUE.,solvph) #else call AMOEBA_check_newstyle_inpcrd(inpcrd,newstyle) if (newstyle) then call AM_RUNMD_get_coords(natom,t,irest,ntb,x(lcrd),x(lvel)) else call getcor(nr,x(lcrd),x(lvel),x(lforce),ntx,box,irest,t,temp0,.TRUE.,solvph) end if #endif if (iamoeba > 0) then natom = natom*am_nbead nrp = nrp*am_nbead nr = nr*am_nbead nr3 = nr3*am_nbead ncopy = am_nbead end if ! M-WJ !if( igb == 0 .and. induced == 1 ) call get_dips(x,nr) ! WJM if is a polarizable model, reading input dipole information if (igb == 0 .and. ipb == 0 .and. induced > 0) call get_dips(x,nr) #ifdef APBS ! APBS initialization if (mdin_apbs) then ! in: natom, coords, charge and radii (from prmtop) ! out: pb charges and pb radii (via apbs_vars module) call apbs_init(natom, x(lcrd), x(l15), x(l97)) end if #endif /* APBS */ #ifdef _XRAY call xray_init() #endif ! ----- SET THE INITIAL VELOCITIES ----- if (ntx <= 3) then call setvel(nr,x(lvel),x(lwinv),tempi,init,iscale,scalm) ! random numbers may have been "used up" in setting the intial ! velocities; re-set the generator so that all nodes are back in ! sync ! DAN ROE: Note master node only here #ifdef MPI if (rem == 0) call amrset(ig) #else call amrset(ig) #endif end if if (belly) call bellyf(natom,ix(ibellygp),x(lvel)) call timer_stop(TIME_RDCRD) if(abfqmmm_param%abfqmmm == 1 .and. ntb > 0) then ! lam81 call iwrap2(abfqmmm_param%n_user_qm, abfqmmm_param%user_qm, x(lcrd), box_center) ! lam81 end if ! lam81 ! --- If we are reading NMR restraints/weight changes, ! read them now: if (nmropt >= 1) then call nmrcal(x(lcrd),f,ih(m04),ih(m02),ix(i02),x(lwinv),enmr, & devdis,devang,devtor,devplpt,devpln,devgendis,temp0,tautp,cut,ntb,x(lnmr01), & ix(inmr02),x(l95),5,6,rk,tk,pk,cn1,cn2, & ag,bg,cg,numbnd,numang,numphi,nimprp, & nhb,natom,natom,ntypes,nres,rad,wel,radhb, & welhb,rwell,isftrp,tgtrmsd,temp0les,-1,'READ') ! Updated 9/2007 by Matthew Seetin to enable plane-point and plane-plane restraints ! --- Determine how many of the torsional parameters ! are impropers call impnum(ix(i46),ix(i56),ix(i48),ix(i58),nphih,nphia, & 0,nptra,nimprp) end if ! -- Set up info related to weight changes for the non-bonds: call nmrrad(rad,wel,cn1,cn2,ntypes,0,0.0d0) call decnvh(asol,bsol,nphb,radhb,welhb) if (iredir(4) > 0) call noeread(x,ix,ih) if (iredir(8) > 0) call alignread(natom, x(lcrd)) if (iredir(9) > 0) call csaread end if ! lam81 <<< QMSTEP=1 BLOCK >>> !--------------------------------------------------------------- ! --- Call FASTWAT, which will tag those bonds which are part ! of 3-point water molecules. Constraints will be effected ! for these waters using a fast analytic routine -- dap. call timer_start(TIME_FASTWT) call fastwat(ih(m04),nres,ix(i02),ih(m02), & nbonh,nbona,ix(iibh),ix(ijbh),ibelly,ix(ibellygp), & iwtnm,iowtnm,ihwtnm,jfastw,ix(iifstwt), & ix(iifstwr),ibgwat,ienwat,ibgion,ienion,iorwat, & 6,natom) call timer_stop(TIME_FASTWT) call getwds(ih(m04), nres , ix(i02) , ih(m02) , & nbonh , nbona , 0 , ix(iibh) , & ix(ijbh) , iwtnm , iowtnm , ihwtnm , & jfastw , ix(iicbh) , req , x(lwinv) , & rbtarg , ibelly , ix(ibellygp), 6 ) ! Assign link atoms between quantum mechanical and molecular mechanical ! atoms if quantum atoms are present. ! After assigning the link atoms, delete all connectivity between the ! QM atoms. if(qmmm_nml%ifqnt) then call identify_link_atoms(nbona,ix(iiba),ix(ijba)) ! Variable QM solvent: ! Store the original bond parameters since we will need to rebuild ! the QM region (delete bonded terms etc) repeatedly if ( qmmm_nml%vsolv > 0 ) then call new(qmmm_vsolv, nbonh, nbona, ntheth, ntheta, nphih, nphia) call qmmm_vsolv_store_parameters(qmmm_vsolv, numbnd, & ix(iibh), ix(ijbh), ix(iicbh), & ix(iiba), ix(ijba), ix(iicba), & ix(i24), ix(i26), ix(i28), ix(i30), & ix(i32), ix(i34), ix(i36), ix(i38), & ix(i40), ix(i42), ix(i44), ix(i46), ix(i48), & ix(i50), ix(i52), ix(i54), ix(i56), ix(i58)) end if if( abfqmmm_param%abfqmmm == 1 ) then ! lam81 if(abfqmmm_param%qmstep == 1 .and. abfqmmm_param%system == 1) then ! lam81 call abfqmmm_allocate_arrays_of_parameters(numbnd, nbonh, nbona, ntheth, ntheta, nphih, nphia) ! lam81 call abfqmmm_store_parameters(ix(iibh), ix(ijbh), ix(iicbh), & ! lam81 ix(iiba), ix(ijba), ix(iicba), & ! lam81 ix(i24), ix(i26), ix(i28), ix(i30), & ! lam81 ix(i32), ix(i34), ix(i36), ix(i38), & ! lam81 ix(i40), ix(i42), ix(i44), ix(i46), ix(i48), & ! lam81 ix(i50), ix(i52), ix(i54), ix(i56), ix(i58), & ! lam81 x(l15), rk, req) ! lam81 else ! lam81 call abfqmmm_set_parameters(numbnd, nbonh, nbona, ntheth, ntheta, nphih, nphia, & ! lam81 ix(iibh), ix(ijbh), ix(iicbh), & ! lam81 ix(iiba), ix(ijba), ix(iicba), & ! lam81 ix(i24), ix(i26), ix(i28), ix(i30), & ! lam81 ix(i32), ix(i34), ix(i36), ix(i38), & ! lam81 ix(i40), ix(i42), ix(i44), ix(i46), ix(i48), & ! lam81 ix(i50), ix(i52), ix(i54), ix(i56), ix(i58), & ! lam81 x(l15), rk, req) ! lam81 call init_extra_pts(ix(iibh),ix(ijbh),ix(iicbh), & ! lam81 ix(iiba),ix(ijba),ix(iicba), & ! lam81 ix(i24),ix(i26),ix(i28),ix(i30), & ! lam81 ix(i32),ix(i34),ix(i36),ix(i38), & ! lam81 ix(i40),ix(i42),ix(i44),ix(i46),ix(i48), & ! lam81 ix(i50),ix(i52),ix(i54),ix(i56),ix(i58), & ! lam81 ih(m06),ix,x,ix(i08),ix(i10), & ! lam81 nspm,ix(i70),x(l75),tmass,tmassinv,x(lmass),x(lwinv),req) ! lam81 end if ! lam81 end if ! lam81 ! Remove bonds between QM atoms from list (Hydrogen) if (nbonh .gt. 0) call setbon(nbonh,ix(iibh),ix(ijbh),ix(iicbh), & ix(ibellygp)) ! Remove bonds between QM atoms from list (Heavy) if (nbona .gt. 0) call setbon(nbona,ix(iiba),ix(ijba),ix(iicba), & ix(ibellygp)) ! Remove angles between QM atoms from list (Hydrogen) if (ntheth .gt. 0) call setang(ntheth,ix(i24),ix(i26),ix(i28),ix(i30), & ix(ibellygp)) ! Remove angles between QM atoms from list (Heavy) if (ntheta .gt. 0) call setang(ntheta,ix(i32),ix(i34),ix(i36),ix(i38),& ix(ibellygp)) ! Remove dihedrals between QM atoms from list (Hydrogen) if (nphih .gt. 0) call setdih(nphih,ix(i40),ix(i42),ix(i44),ix(i46), & ix(i48), ix(ibellygp)) ! Remove dihedrals between QM atoms from list (Heavy) if (nphia .gt. 0) call setdih(nphia,ix(i50),ix(i52),ix(i54),ix(i56), & ix(i58), ix(ibellygp)) ! Remove CHARMM energy terms from QM region call charmm_filter_out_qm_atoms() ! Now we should work out the type of each quantum atom present. ! This is used for our arrays of pre-computed parameters. It is ! essentially a re-basing of the atomic numbers and is done to save ! memory. Note: qm_assign_atom_types will allocate the qm_atom_type ! array for us. Only the master calls this routine. All other ! threads get this allocated and broadcast to them by the mpi setup ! routine. call qm_assign_atom_types ! Set default QMMM MPI parameters - for single cpu operation. ! These will get overwritten by qmmm_mpi_setup if MPI is on. qmmm_mpi%commqmmm_master = master qmmm_mpi%numthreads = 1 qmmm_mpi%mytaskid = 0 qmmm_mpi%natom_start = 1 qmmm_mpi%natom_end = natom qmmm_mpi%nquant_nlink_start = 1 qmmm_mpi%nquant_nlink_end = qmmm_struct%nquant_nlink ! Now we know how many link atoms we can allocate the scf_mchg array... allocate(qm2_struct%scf_mchg(qmmm_struct%nquant_nlink), stat = ier) REQUIRE(ier == 0) ! Deallocated in deallocate qmmm ! We can also allocate ewald_memory if (qmmm_nml%qm_ewald > 0 ) then call allocate_qmewald(natom) end if if (qmmm_nml%qmgb == 2 ) then ! lam81 call allocate_qmgb(qmmm_struct%nquant_nlink) end if ! lam81 allocate( qmmm_struct%dxyzqm(3, qmmm_struct%nquant_nlink), stat = ier ) REQUIRE(ier == 0) ! Deallocated in deallocate qmmm else if(abfqmmm_param%abfqmmm == 1) then ! lam81 call abfqmmm_set_parameters(numbnd, nbonh, nbona, ntheth, ntheta, nphih, nphia, & ! lam81 ix(iibh), ix(ijbh), ix(iicbh), & ! lam81 ix(iiba), ix(ijba), ix(iicba), & ! lam81 ix(i24), ix(i26), ix(i28), ix(i30), & ! lam81 ix(i32), ix(i34), ix(i36), ix(i38), & ! lam81 ix(i40), ix(i42), ix(i44), ix(i46), ix(i48), & ! lam81 ix(i50), ix(i52), ix(i54), ix(i56), ix(i58), & ! lam81 x(l15), rk, req) ! lam81 call init_extra_pts(ix(iibh),ix(ijbh),ix(iicbh), & ! lam81 ix(iiba),ix(ijba),ix(iicba), & ! lam81 ix(i24),ix(i26),ix(i28),ix(i30), & ! lam81 ix(i32),ix(i34),ix(i36),ix(i38), & ! lam81 ix(i40),ix(i42),ix(i44),ix(i46),ix(i48), & ! lam81 ix(i50),ix(i52),ix(i54),ix(i56),ix(i58), & ! lam81 ih(m06),ix,x,ix(i08),ix(i10), & ! lam81 nspm,ix(i70),x(l75),tmass,tmassinv,x(lmass),x(lwinv),req) ! lam81 end if !if (qmmm_nml%ifqnt) ! --- Open the data dumping files and position it depending ! on the type of run: if (abfqmmm_param%qmstep == 1 .and. abfqmmm_param%system == 1) then ! lam81 #ifdef MPI ! adaptive QM/MM (qmmm_nml%vsolv=2) via multisander ! all groups have identical coords and velocities ! only master of first group needs to dump results if ( (qmmm_nml%vsolv < 2) .or. (worldrank == 0) ) then #endif call open_dump_files #ifdef MPI end if #endif end if ! lam81 if (master) call amflsh(6) ! --- end of master process setup --- end if masterwork ! (master) # if defined(RISMSANDER) call rism_init(commsander) # endif /* RISMSANDER */ #ifdef MPI call mpi_barrier(commsander,ierr) ! =========================== AMBER/MPI =========================== ! NOTE: in the current AMBER/MPI implementation, two means of ! running in parallel within sander are supported. The value ! of mpi_orig determines which approach is used. ! This is turned on when minimization (imin .ne. 0) is requested, ! and is otherwise off. ! When running the mpi_orig case, a variable notdone is now ! set by the master and determines when to exit the force() ! loop. When the master has finished calling force, the ! master changes notdone to 0 and broadcasts the data one more ! time to signal end of the loop. force() is modified so that ! in the mpi_orig case, an initial broadcast is done to receive ! the value from the master to decide whether to do the work or ! simply exit. ! ...set up initial data and send all needed data to other nodes, ! now that the master has it ! First, broadcast parameters in memory.h, so that all processors ! will know how much memory to allocate: call mpi_bcast(natom,BC_MEMORY,mpi_integer,0,commsander,ierr) ! -- ti decomp call mpi_bcast(idecomp,1,mpi_integer,0,commsander,ierr) call mpi_bcast(nat,1,mpi_integer,0,commsander,ierr) call mpi_bcast(nrs,1,mpi_integer,0,commsander,ierr) ! ----- Set up integer stack initial size -------------- call mpi_bcast(lastist,1,mpi_integer,0,commsander,ierr) call mpi_bcast(lastrst,1,mpi_integer,0,commsander,ierr) call mpi_bcast(qmmm_struct%abfqmmm,1,mpi_integer,0,commsander,ierr) ! lam81 call mpi_bcast(abfqmmm_param%abfqmmm,1,mpi_integer,0,commsander,ierr) ! lam81 if(abfqmmm_param%abfqmmm == 1) then ! lam81 call mpi_bcast(abfqmmm_param%natom,1,mpi_integer,0,commsander,ierr) ! lam81 call mpi_bcast(abfqmmm_param%qmstep,1,mpi_integer,0,commsander,ierr) ! lam81 call mpi_bcast(abfqmmm_param%maxqmstep,1,mpi_integer,0,commsander,ierr) ! lam81 call mpi_bcast(abfqmmm_param%system,1,mpi_integer,0,commsander,ierr) ! lam81 end if ! lam81 if(abfqmmm_param%abfqmmm == 1) then ! lam81 if(abfqmmm_param%qmstep == 1 .and. abfqmmm_param%system == 1) then ! lam81 call mpi_bcast(abfqmmm_param%r_core_in,1,mpi_double_precision,0,commsander,ierr) ! lam81 call mpi_bcast(abfqmmm_param%r_core_out,1,mpi_double_precision,0,commsander,ierr) ! lam81 call mpi_bcast(abfqmmm_param%r_qm_in,1,mpi_double_precision,0,commsander,ierr) ! lam81 call mpi_bcast(abfqmmm_param%r_qm_out,1,mpi_double_precision,0,commsander,ierr) ! lam81 call mpi_bcast(abfqmmm_param%r_buffer_in,1,mpi_double_precision,0,commsander,ierr) ! lam81 call mpi_bcast(abfqmmm_param%r_buffer_out,1,mpi_double_precision,0,commsander,ierr) ! lam81 call mpi_bcast(abfqmmm_param%mom_cons_type,1,mpi_integer,0,commsander,ierr) ! lam81 call mpi_bcast(abfqmmm_param%mom_cons_region,1,mpi_integer,0,commsander,ierr) ! lam81 if (.not. master) then ! lam81 allocate(abfqmmm_param%x(3*natom), stat=ier) ! lam81 REQUIRE(ier==0) ! lam81 allocate(abfqmmm_param%id(natom), stat=ier) ! lam81 REQUIRE(ier==0) ! lam81 allocate(abfqmmm_param%mass(natom), stat=ier) ! lam81 REQUIRE(ier==0) ! lam81 end if ! lam81 allocate(abfqmmm_param%v(3*natom+iscale), stat=ier) ! lam81 REQUIRE(ier==0) ! lam81 allocate(abfqmmm_param%f(3*natom+iscale), stat=ier) ! lam81 REQUIRE(ier==0) ! lam81 allocate(abfqmmm_param%f1(3*natom+iscale), stat=ier) ! lam81 REQUIRE(ier==0) ! lam81 allocate(abfqmmm_param%f2(3*natom+iscale), stat=ier) ! lam81 REQUIRE(ier==0) ! lam81 end if ! lam81 call mpi_bcast(abfqmmm_param%x,3*natom,mpi_double_precision,0,commsander,ierr) ! lam81 call mpi_bcast(abfqmmm_param%id,natom,mpi_integer,0,commsander,ierr) ! lam81 call mpi_bcast(abfqmmm_param%mass,natom,mpi_double_precision,0,commsander,ierr) ! lam81 end if ! lam81 if (abfqmmm_param%qmstep == 1 .and. abfqmmm_param%system == 1) then ! lam81 call stack_setup() else ! lam81 call deallocate_stacks ! lam81 call stack_setup() ! lam81 end if ! lam81 ! GMS: Broadcast parameters from module 'molecule' call mpi_bcast(mol_info%natom,1,mpi_integer,0,commsander,ierr) call mpi_bcast(mol_info%nres,1,mpi_integer,0,commsander,ierr) call mpi_barrier(commsander,ierr) ! ---allocate memory on the non-master nodes: if (.not. master) then if (abfqmmm_param%qmstep /= 1 .or. abfqmmm_param%system /= 1) then ! lam81 deallocate(x,ix,ipairs,ih) ! lam81 end if ! lam81 allocate( x(1:lastr), stat = ier ) REQUIRE( ier == 0 ) allocate( ix(1:lasti), stat = ier ) REQUIRE( ier == 0 ) allocate( ipairs(1:lastpr), stat = ier ) REQUIRE( ier == 0 ) allocate( ih(1:lasth), stat = ier ) REQUIRE( ier == 0 ) ! GMS: ! Allocate space for molecule module ! arrays in the other nodes if (abfqmmm_param%qmstep == 1 .and. abfqmmm_param%system == 1) then ! lam81 call allocate_molecule() else ! lam81 call deallocate_molecule() ! lam81 call allocate_molecule() ! lam81 end if ! lam81 ! -------------------------- ! -- ti decomp if (idecomp > 0) then call allocate_int_decomp(natom, nres) if (idecomp == 1 .or. idecomp == 2) then call allocate_real_decomp(nrs) else if( idecomp == 3 .or. idecomp == 4 ) then call allocate_real_decomp(npdec*npdec) end if end if end if ! ( .not.master ) !GMS: Broadcast arrays from module 'molecule' call mpi_bcast(mol_info%natom_res,mol_info%nres,MPI_INTEGER,0,commsander,ierr) call mpi_bcast(mol_info%atom_to_resid_map,mol_info%natom,MPI_INTEGER,0,commsander,ierr) call mpi_bcast(mol_info%atom_mass,mol_info%natom, MPI_DOUBLE_PRECISION, 0, commsander, ier) if(idecomp == 1 .or. idecomp == 2) then call mpi_bcast(jgroup,nat,MPI_INTEGER,0,commsander,ierr) end if if(icfe == 0 .and. (idecomp ==3 .or. idecomp == 4)) then call mpi_bcast(jgroup,natom,MPI_INTEGER,0,commsander,ierr) call mpi_bcast(indx,nres,MPI_INTEGER,0,commsander,ierr) end if call startup_groups(ierr) call startup(x,ix,ih) ! +---------------------------------------------------------------------------+ ! | Broadcast EVB/PIMD inputs/parameters to all PEs | ! +:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::+ ! | Note: The masters have all required EVB/PIMD inputs/parameters via call | ! | to mdread2 ( evb_input, evb_init, evb_pimd_init ). For EVB/PIMD, all | ! | PEs need the inputs/parameters ... so we perform this initialization | ! | again for all PEs besides the masters. The alternative is to use | ! | MPI_BCAST. | ! +---------------------------------------------------------------------------+ call mpi_bcast( ievb , 1, MPI_INTEGER, 0, commworld, ierr ) call mpi_bcast( ipimd, 1, MPI_INTEGER, 0, commworld, ierr ) ! KFW if (ievb /= 0) then call mpi_bcast(evbin, MAX_FN_LEN, MPI_CHARACTER, 0, commworld, ierr) if (.not. master) then call evb_input call evb_init end if end if ! KFW # if defined(LES) call mpi_bcast ( ncopy, 1, MPI_INTEGER, 0, commworld, ierr ) call mpi_bcast ( cnum(1:natom), natom, MPI_INTEGER, 0, commworld, ierr ) call mpi_bcast ( evbin, MAX_FN_LEN, MPI_CHARACTER, 0, commworld, ierr ) # endif /* LES */ if (ievb /= 0) then # if defined(LES) !call mpi_bcast ( mastersize, 1, MPI_INTEGER, 0, commworld, ierr ) !call mpi_bcast ( jobs_per_node, 1, MPI_INTEGER, 0, commworld, ierr ) !call mpi_bcast ( nsize, 1, MPI_INTEGER, 0, commworld, ierr ) if (ipimd > 0 .and. .not. master) then call evb_input call evb_init !call evb_pimd_init end if call evb_pimd_init !call mpi_bcast ( master_worldrank, mastersize, MPI_INTEGER, 0, commworld, ierr ) !call mpi_bcast ( PE_slice, jobs_per_node*nsize, MPI_INTEGER, 0, commworld, ierr ) # endif /* LES */ call evb_bcast call evb_alloc end if ! +---------------------------------------------------------------------------+ ! | Obtain B vector for Schlegel's distributed Gaussian method | ! +---------------------------------------------------------------------------+ if (trim(adjustl(xch_type)) == "dist_gauss") call schlegel_dg if (ifsc /= 0) then ! multi-CPU minimization does not work with soft core ! if (imin > 0 .and. numtasks > 1) then call sander_bomb('imin > 0 and numtasks > 1', & 'TI minimizations cannot be performed with > 2 CPUs','') end if call setup_sc(natom, nres, ih(m04), ih(m06), & ix(i02), ih(m02), x(lcrd), ntypes, clambda, nstlim) if (ntp > 0 .and. master) then ! Check which molecules are perturbed in NPT runs call sc_check_perturbed_molecules(nspm, ix(i70)) end if ! -- ti decomp if (idecomp > 0) then if (sanderrank == 0) call build_dec_mask(natom) call mpi_bcast(decmask, natom, MPI_INTEGER, 0, commsander, ierr) end if ! Make sure all common atoms have the same v (that of V0) in TI runs if ( ifsc /= 2 ) then if (master) call sc_sync_x(x(lvel),nr3) !call mpi_barrier(commsander,ierr) if (numtasks > 1) then call mpi_bcast(nr3,1,MPI_INTEGER,0,commsander,ierr) call mpi_bcast(x(lvel),nr3,MPI_DOUBLE_PRECISION,0,commsander,ierr) end if end if if (tishake /= 0) call setnoshake_sc(ix,ntc,num_noshake,master) else extra_atoms=0 end if if (ifmbar /= 0) then call setup_mbar(nstlim) end if if (.not. master .and. igb == 0 .and. ipb == 0) then if (abfqmmm_param%qmstep == 1 .and. abfqmmm_param%system == 1) then ! lam81 call nblist_allocate(natom,ntypes,num_direct,numtasks) else ! lam81 call nblist_deallocate ! lam81 call nblist_allocate(natom,ntypes,num_direct,numtasks) ! lam81 end if ! lam81 end if ! --- Allocate memory for GB on the non-master nodes: if( .not.master ) then if ((igb /= 0 .and. igb /= 10 .and. ipb == 0) .or. hybridgb>0 .or. icnstph.gt.1) then ! lam81 if (abfqmmm_param%qmstep == 1 .and. abfqmmm_param%system == 1) then ! lam81 call allocate_gb( natom, ncopy ) else ! lam81 call deallocate_gb ! lam81 call allocate_gb( natom, ncopy ) ! lam81 end if ! lam81 end if ! lam81 end if ! ( .not.master ) nr = nrp nr3 = 3*nr belly = ibelly > 0 ! ---Do setup for QMMM in parallel if ifqnt is on - this is essentially ! everything in the QMMM namelist and a few other bits and pieces. ! ---Note, currently only master node knows if qmmm_nml%ifqnt is ! on so we need to broadcast this first and then make decisions based ! on this. call mpi_bcast(qmmm_nml%ifqnt, 1, mpi_logical, 0, commsander, ierror) if (qmmm_nml%ifqnt) then ! Broadcast all of the stuff in qmmm_nml and allocate the relevant ! arrays on all processors. This will also set up information for ! openmp on the master processor if it is in use. call qmmm_mpi_setup( master, natom ) if (qmmm_nml%qm_ewald > 0 .and. .not. master) then call allocate_qmewald(natom) end if if (qmmm_nml%qmgb==2 .and. .not. master) then call allocate_qmgb(qmmm_struct%nquant_nlink) end if if (.not. master) then allocate( qmmm_struct%dxyzqm(3, qmmm_struct%nquant_nlink), stat = ier ) REQUIRE(ier == 0) !Deallocated in deallocate qmmm end if end if ! --- END QMMM MPI SETUP --- ! DAN ROE: Note: all nodes are calling this. amrset(ig) has already been called ! by the master node (twice if initial velocities are set, ntx<=3). ! DAN ROE: REMD: Now need to call amrset on all child threads, masters have called ! it above before initial coord read if (abfqmmm_param%qmstep == 1 .and. abfqmmm_param%system == 1) then ! lam81 if(rem == 0) then call amrset(ig+1) else if (.not. master) call amrset(ig + 17 * nodeid) endif if (nmropt >= 1) & call nmrcal(x(lcrd),f,ih(m04),ih(m02),ix(i02),x(lwinv),enmr, & devdis,devang,devtor,devplpt,devpln,devgendis,temp0,tautp,cut,ntb,x(lnmr01), & ix(inmr02),x(l95),5,6,rk,tk,pk,cn1,cn2, & ag,bg,cg,numbnd,numang,numphi,nimprp, & nhb,natom,natom,ntypes,nres,rad,wel,radhb, & welhb,rwell,isftrp,tgtrmsd,temp0les,-1,'MPI ') end if ! lam81 ! Updated 9/2007 by Matthew Seetin to enable plane-point and plane-plane restraints call mpi_bcast(lmtmd01,1,mpi_integer,0, commsander,ierr) call mpi_bcast(imtmd02,1,mpi_integer,0, commsander,ierr) if (itgtmd == 2) call mtmdcall(emtmd,x(lmtmd01),ix(imtmd02),x(lcrd),x(lforce),& ih(m04),ih(m02),ix(i02),ih(m06),x(lmass),natom,nres,'MPI ') ! ---------------- Check system is neutral and print warning message ------ ! ---------------- adjust charges for roundoff error. ------ if( igb == 0 .and. ipb == 0 .and. iyammp == 0 ) then if ( icfe == 0 ) then call check_neutral(x(l15),natom) else call ti_check_neutral(x(l15),natom) end if end if ! tell all nodes if this is a SEBOMD run call mpi_bcast(sebomd_obj%do_sebomd, 1, MPI_LOGICAL, 0, commsander, ierr) if (sebomd_obj%do_sebomd) then ! transfer SEBOMD info to the nodes call sebomd_bcast_obj() ! open necessary files if (master) then call sebomd_open_files end if ! initialize SEBOMD arrays call init_sebomd_arrays(natom) end if ! ---------------- Old parallel for minimization ---------------------- if (imin /= 0) then mpi_orig = .true. notdone = 1 else mpi_orig = .false. end if if (mpi_orig .and. .not. master) then ! ...all nodes only do the force calculations (JV) ! Minimisation so only only master gets past the loop below ! hence need to zero QM charges on non-master threads here. if (qmmm_nml%ifqnt) then ! Apply charge correction if required. if (qmmm_nml%adjust_q>0) then call qmmm_adjust_q(qmmm_nml%adjust_q, natom, qmmm_struct%nquant, qmmm_struct%nquant_nlink, & qmmm_struct%nlink, x(L15), & qmmm_struct%iqmatoms, qmmm_nml%qmcharge, qmmm_struct%atom_mask, & qmmm_struct%mm_link_mask, master,x(LCRD), qmmm_nml%vsolv) end if ! At this point we can also fill the qmmm_struct%scaled_mm_charges ! array - we only need to do this once as the charges are constant ! during a run. Having a separate array of scaled charges saves us ! having to do it on every qmmm routine call. Do this BEFORE zeroing ! the QM charges since that routine take care of these values as well. do i = 1, natom qmmm_struct%scaled_mm_charges(i) = x(L15+(i-1)) * INV_AMBER_ELECTROSTATIC & * qmmm_nml%chg_lambda ! charge scaling factor for FEP end do ! Zero out the charges on the quantum mechanical atoms call qm_zero_charges(x(L15),qmmm_struct%scaled_mm_charges,.true.) if (qmmm_struct%nlink > 0 ) then ! We need to exclude all electrostatic ! interactions with MM link pairs, both QM-MM and MM-MM. Do this by ! zeroing the MM link pair charges in the main charge array. ! These charges are stored in qmmm_struct%mm_link_pair_resp_charges in case ! they are later needed. call qm_zero_mm_link_pair_main_chg(qmmm_struct%nlink,qmmm_struct%link_pairs,x(L15), & qmmm_struct%scaled_mm_charges,.true.) end if end if if (igb == 7 .or. igb == 8) call igb7_init(natom, ncopy, x(l97)) !x(l97) is rborn() !Hai Nguyen: add igb == 8 here if ( ifcr /= 0 ) call cr_allocate( master, natom ) n_force_calls = 0 do while( notdone == 1 ) n_force_calls = n_force_calls+1 call force(x,ix,ih,ipairs,x(lcrd),x(lforce),ene,vir, & x(l96), x(l97), x(l98), x(l99), qsetup, & do_list_update,n_force_calls) end do ! Deallocate and return goto 999 end if ! ---------------------------------------------------------------------- if (master) then if(abfqmmm_param%abfqmmm /=1) write(6, '(a,i4,a,/)') '| Running AMBER/MPI version on ',numtasks, ' nodes' ! lam81 !BTREE is selected by default if cpu is a power of two. !The number of processes is required to be a power of two for Btree !Print a warning about inefficiency with cpus not being a power of 2. if (numtasks > 1 .and. logtwo(numtasks) <= 0) then if(abfqmmm_param%abfqmmm /=1) write(6, '(a,i4,a,/)') '| WARNING: The number of processors is not a power of 2' ! lam81 if(abfqmmm_param%abfqmmm /=1) write(6, '(a,i4,a,/)') '| this may be inefficient on some systems.' ! lam81 end if end if if (master .and. numgroup > 1) write(6, '(a,i4,a,i4,a,i4,a)') & '| MULTISANDER: ', numgroup, ' groups. ', & numtasks, ' processors out of ', worldsize, ' total.' if (master) call amflsh(6) ! ========================= END AMBER/MPI ========================= #else if(abfqmmm_param%abfqmmm == 1) then ! lam81 if(abfqmmm_param%qmstep == 1 .and. abfqmmm_param%system == 1) then ! lam81 allocate(abfqmmm_param%v(3*natom+iscale), stat=ier) ! lam81 REQUIRE(ier==0) ! lam81 allocate(abfqmmm_param%f(3*natom+iscale), stat=ier) ! lam81 REQUIRE(ier==0) ! lam81 allocate(abfqmmm_param%f1(3*natom+iscale), stat=ier) ! lam81 REQUIRE(ier==0) ! lam81 allocate(abfqmmm_param%f2(3*natom+iscale), stat=ier) ! lam81 REQUIRE(ier==0) ! lam81 end if ! lam81 end if ! lam81 ! debug needs to copy charges at start and they can't change later ! ---------------- Check system is neutral and print warning message ------ ! ---------------- adjust charges for roundoff error. ------ if( igb == 0 .and. ipb == 0 .and. iyammp == 0 ) call check_neutral(x(l15),natom) if (abfqmmm_param%qmstep == 1 .and. abfqmmm_param%system == 1) then ! lam81 call amrset(ig+1) end if ! lam81 if (abfqmmm_param%qmstep == 1 .and. abfqmmm_param%system == 1) then ! lam81 call stack_setup() else ! lam81 call deallocate_stacks ! lam81 call stack_setup() ! lam81 end if ! lam81 if (sebomd_obj%do_sebomd) then ! open necessary files call sebomd_open_files ! initialize SEBOMD arrays call init_sebomd_arrays(natom) endif #endif /* MPI */ #ifdef OPENMP ! If -openmp was specified to configure_amber then -DOPENMP is defined and the ! threaded version of MKL will have been linked in. It is important here that ! we set the default number of openmp threads for MKL to be 1 to stop conflicts ! with threaded vectorization routines when running in parallel etc. ! Individual calls to MKL from routines that know what they are doing - e.g. ! QMMM calls to diagonalizers etc can increase this limit as long as they ! put it back afterwards. call omp_set_num_threads(1) ! If we are using openmp for matrix diagonalization print some information. if (qmmm_nml%ifqnt .and. master) call qm_print_omp_info() #endif ! allocate memory for crg relocation if (ifcr /= 0) call cr_allocate( master, natom ) ! initialize LIE module if used if ( ilrt /= 0 ) then call setup_linear_response(natom,nres,ih(m04),ih(m06),ix(i02),ih(m02),x(lcrd),x(l15), & ntypes, ix(i04), ix(i06), cn1, cn2, master) end if call date_and_time( setup_end_date, setup_end_time ) ! Initialize the printing of ongoing time and performance summaries. ! We do this quite late here to avoid including all the setup time. if (master) call print_ongoing_time_summary(0,0,0.0d0,7) ! ---------------------------------------------------------------------- ! Now do the dynamics or minimization. ! ---------------------------------------------------------------------- if (igb == 7 .or. igb == 8 ) call igb7_init(natom, ncopy, x(l97)) !x(l97) is rborn() !Hai Nguyen: add igb ==8 here if (qmmm_nml%ifqnt) then ! Apply charge correction if required. if (qmmm_nml%adjust_q>0) then call qmmm_adjust_q(qmmm_nml%adjust_q, natom, qmmm_struct%nquant, qmmm_struct%nquant_nlink, & qmmm_struct%nlink, x(L15), & qmmm_struct%iqmatoms, qmmm_nml%qmcharge, qmmm_struct%atom_mask, & qmmm_struct%mm_link_mask, master,x(LCRD), qmmm_nml%vsolv) end if ! At this point we can also fill the qmmm_struct%scaled_mm_charges ! array - we only need to do this once as the charges are constant ! during a run. Having a separate array of scaled charges saves us ! having to do it on every qmmm routine call. Do this BEFORE zeroing ! the QM charges since that routine take care of these values as well. do i = 1, natom qmmm_struct%scaled_mm_charges(i) = x(L15+(i-1)) * INV_AMBER_ELECTROSTATIC & * qmmm_nml%chg_lambda ! charge scaling factor for FEP end do ! Zeroing of QM charges MUST be done AFTER call to check_neutral. ! Zero out the charges on the quantum mechanical atoms. call qm_zero_charges(x(L15),qmmm_struct%scaled_mm_charges,.true.) if (qmmm_struct%nlink > 0) then ! We need to exclude all electrostatic ! interactions with MM link pairs, both QM-MM and MM-MM. Do this by ! zeroing the MM link pair charges in the main charge array. ! These charges are stored in qmmm_struct%mm_link_pair_resp_charges in case ! they are later needed. call qm_zero_mm_link_pair_main_chg(qmmm_struct%nlink,qmmm_struct%link_pairs,x(L15), & qmmm_struct%scaled_mm_charges,.true.) end if end if ! Use the debugf namelist to activate call debug_frc(x,ix,ih,ipairs,x(lcrd),x(lforce),cn1,cn2,qsetup) ! Prepare for SGLD simulation if (isgld > 0) call psgld(natom,x(lmass),x(lvel), rem) ! Prepare for EMAP constraints if (temap) call pemap(dt,temp0,x,ix,ih) ! Prepare for Isotropic periodic sum of nonbonded interaction if (ips .gt. 0) call ipssys(natom,ntypes,ntb,x(l15), & cut,cn1,cn2,ix(i04),ix(i06),x(lcrd)) if (master .and. (.not. qmmm_nml%ifqnt) .and. (abfqmmm_param%abfqmmm /= 1)) & ! lam81 write(6,'(/80(''-'')/,'' 4. RESULTS'',/80(''-'')/)') ! Set up the MC barostat if requested if (ntp > 0 .and. barostat == 2) call mcbar_setup(ig) ! Input flag imin determines the type of calculation: MD, minimization, ... select case (imin) case (0) ! --- Dynamics: call timer_start(TIME_RUNMD) ! Set up AMD if (iamd .gt. 0) then call amd_setup(ntwx) endif ! Set up scaledMD if (scaledMD .gt. 0) then call scaledMD_setup(ntwx) endif if (ipimd > 0) then call pimd_init(natom,x(lmass),x(lwinv),x(lvel),ipimd) end if if(ineb>0) call neb_init() #ifdef MPI ! ----===== REMD =====---- ! If this is not a REMD run, runmd is called only once. ! If this is a REMD run, runmd is called 0 to numexchg times, ! where the 0th runmd is just for getting initial PEs (no dynamics). if (rem == 0) then ! Not a REMD run. runmd will be called once. loop = 0 else if (rem > 0) then ! This is a REMD run. runmd will be called numexchg times. loop = numexchg ! Set up temptable, open remlog, etc. call remd1d_setup(numexchg, hybridgb, numwatkeep, & temp0, mxvar, natom, ig, solvph) else ! Multi-D REMD run loop = numexchg call multid_remd_setup(numexchg, hybridgb, numwatkeep, & temp0, mxvar, natom, ig, solvph, irest) ! Now set up REMD indices for traj/restart writes. Only do this on ! master since only master handles writes. if (master) call setup_remd_indices end if ! Replica run setup ! Loop over REMD exchanges do mdloop = 0, loop ! ----===== REMD EXCHANGE HANDLING =====---- ! Note: mdloop==0 is just used to get initial energies for the ! first exchange. if (rem < 0 .and. mdloop > 0) then call multid_remd_exchange(x, ix, ih, ipairs, qsetup, & do_list_update, temp0, solvph) else if ((rem == 1 .or. rem == 2) .and. mdloop > 0) then call remd_exchange(1, rem, x(lcrd), x(lvel), x(lmass), & nr3, natom, nr, temp0) else if (rem == 3 .and. mdloop > 0) then call hremd_exchange(1, x, ix, ih, ipairs, qsetup, do_list_update) ! force was called inside hremd_exchange, so call nmrdcp ! to decrement the NMR counter, since this should not count ! as a real step. This is OK, since the counter got ! incremented at the _very_ end of nmrcal, so we haven't already ! printed an unwanted value (JMS 2/12) if (nmropt /= 0) call nmrdcp else if (rem == 4 .and. mdloop > 0) then call ph_remd_exchange(1, solvph) end if ! rem>0 and mdloop>0 ! ----===== END REMD EXCHANGE HANDLING =====---- #ifdef VERBOSE_REMD if (rem > 0 .and. mdloop .eq. 0 .and. master) & write (6,'(a,i4)') "REMD: Getting initial energy for replica ",repnum #endif /* VERBOSE_REMD */ #endif /* MPI */ #ifndef DISABLE_NCSU if(abfqmmm_param%qmstep == 1 .and. abfqmmm_param%system == 1) then ! lam81 call ncsu_on_sander_init(ih, x(lmass), x(lcrd), rem) end if ! lam81 #endif /* DISABLE_NCSU */ if (beeman_integrator == 1) then call AM_RUNMD(ix,ih,ipairs, & x(lwinv),x(lmass),x, & x(lcrd),x(lvel),x(lforce),qsetup) else call runmd(x,ix,ih,ipairs, & x(lcrd),x(lwinv),x(lmass),x(lforce), & x(lvel),x(lvel2),x(l45),x(lcrdr), & x(l50),x(l95),ix(i70),x(l75), & erstop,qsetup) end if !beeman_integrator == 1 #ifndef DISABLE_NCSU if(abfqmmm_param%qmstep == 1 .and. abfqmmm_param%system == 1) then ! lam81 call ncsu_on_sander_exit() end if ! lam81 #endif /* DISABLE_NCSU */ ! ----===== END REMD =====---- #ifdef MPI end do ! Loop over REMD exchanges (mdloop) ! Cleanup REMD files. if (rem /= 0) call remd_cleanup() #endif call timer_stop(TIME_RUNMD) if (master) call amflsh(6) if (erstop) then ! This error condition stems from subroutine shake; ! furthermore, it seems that erstop can never be true since shake ! can never return with its third last argument, niter, equal to 0. ! SRB, Sep 24, 2003 if (master) then write(6, *) 'FATAL ERROR' end if call mexit(6,1) end if case (1) !--- Minimization: ! Input flag ntmin determines the method of minimization select case (ntmin) case (0, 1, 2) call runmin(x,ix,ih,ipairs,x(lcrd),x(lforce),x(lvel), & ix(iibh),ix(ijbh),x(l50),x(lwinv),ix(ibellygp), & x(l95),ene, carrms, qsetup) ! If a conventional minimisation is being done, ! the restart file is written inside the runmin routine. case (LMOD_NTMIN_XMIN) write(6,'(a)') ' LMOD XMIN Minimization.' write(6,'(a)') '' write(6,'(a)') ' Note: Owing to the behaviour of the XMIN algorithm,' write(6,'(a)') ' coordinates in the trajectory and intermediate' write(6,'(a)') ' restart files will not match up with energies' write(6,'(a)') ' in the mdout and mdinfo files. The final energy' write(6,'(a)') ' and final coordinates do match.' write(6,'(a)') '' xmin_iter = 0 call run_xmin( x, ix, ih, ipairs, & x(lcrd), x(lforce), ene, qsetup, xmin_iter, ntpr ) if (master) call minrit(0,nrp,ntxo,x(lcrd)) ! Write the restart file case (LMOD_NTMIN_LMOD) write(6,'(a)') ' LMOD LMOD Minimization.' write(6,'(a)') '' write(6,'(a)') ' Note: Owing to the behaviour of the XMIN algorithm,' write(6,'(a)') ' coordinates in the trajectory and intermediate' write(6,'(a)') ' restart files will not match up with energies' write(6,'(a)') ' in the mdout and mdinfo files. The final energy' write(6,'(a)') ' and final coordinates do match.' write(6,'(a)') '' call run_lmod( x, ix, ih, ipairs, & x(lcrd), x(lforce), ene, qsetup ) if (master) call minrit(0,nrp,ntxo,x(lcrd)) ! Write the restart file case default ! invalid ntmin ! ntmin input validation occurs in mdread.f ASSERT( .false. ) end select case (5) ! ---carlos modified for reading trajectories (trajene option) write (6,*) "POST-PROCESSING OF TRAJECTORY ENERGIES" ! ---read trajectories and calculate energies for each frame call trajene(x,ix,ih,ipairs,ene,ok,qsetup) if (.not. ok) then write (6,*) 'error in trajene()' call mexit(6,1) end if case default ! invalid imin ! imin input validation should be transferred to mdread.f write(6,'(/2x,a,i3,a)') 'Error: Invalid IMIN (',imin,').' ASSERT( .false. ) end select #ifdef MPI /* SOFT CORE */ if (master) then if (icfe /=0 .and. ifsc == 1) call summarize_ti_changes(natom,resat) end if #endif ! finish up EMAP if (temap) call qemap() if (abfqmmm_param%abfqmmm /= 1) then ! lam81 exit ! lam81 else ! lam81 #ifdef MPI /* lam81 */ call mpi_barrier(commsander,ierr) ! lam81 #endif /* lam81 */ if (abfqmmm_param%qmstep /= abfqmmm_param%maxqmstep .or. abfqmmm_param%system /= 2) then ! lam81 if(qmmm_nml%ifqnt) call deallocate_qmmm(qmmm_nml, qmmm_struct, qmmm_vsolv, qm2_params) ! lam81 end if ! lam81 call deallocate_m1m2m3() ! lam81 if (abfqmmm_param%system == 2 .and. master) then ! lam81 call abfqmmm_write_idrst() ! lam81 call abfqmmm_write_pdb(x(lcrd),ix(i70)) ! lam81 end if ! lam81 call abfqmmm_next_step() ! lam81 end if ! lam81 end do ! lam81 if(abfqmmm_param%abfqmmm == 1) then ! lam81 deallocate(abfqmmm_param%id, stat=ier) ! lam81 REQUIRE( ier == 0 ) ! lam81 deallocate(abfqmmm_param%v, stat=ier) ! lam81 REQUIRE( ier == 0 ) ! lam81 deallocate(abfqmmm_param%f, stat=ier) ! lam81 REQUIRE( ier == 0 ) ! lam81 deallocate(abfqmmm_param%f1, stat=ier) ! lam81 REQUIRE( ier == 0 ) ! lam81 deallocate(abfqmmm_param%f2, stat=ier) ! lam81 REQUIRE( ier == 0 ) ! lam81 if(master) then ! lam81 deallocate(abfqmmm_param%isqm, stat=ier) ! lam81 REQUIRE( ier == 0 ) ! lam81 endif ! lam81 end if ! lam81 ! -- calc time spent running vs setup call timer_stop(TIME_TOTAL) call wallclock( time1 ) call date_and_time( final_date, final_time ) call profile_time( time1 - time0, num_calls_nblist, profile_mpi) #ifdef MPI ! =========================== AMBER/MPI =========================== ! Set and broadcast notdone in mpi_orig case to inform ! other nodes that we are finished calling force(). (tec3) if (mpi_orig) then notdone = 0 call mpi_bcast(notdone,1,mpi_integer,0, commsander,ierr) end if ! ========================= END AMBER/MPI ========================= #endif ! ========================= PUPIL INTERFACE ========================= #ifdef PUPIL_SUPPORT ! Finalize Corba Interface puperror = 0 call killcorbaintfc(puperror) if (puperror /= 0) then write(6,*) 'Error ending PUPIL CORBA interface.' end if write(6,'(a)') 'PUPIL CORBA interface finalized.' pupactive = .false. #endif ! ========================= PUPIL INTERFACE ========================= #ifdef _XRAY !(out_lun,residue_pointer,residue_label,atom_name,coor,num_bonds,ibond,jbond) call xray_fini() #endif call amflsh(6) if (master) then #ifdef MPI ! adaptive QM/MM (qmmm_nml%vsolv=2) via multisander ! all groups have identical coords and velocities ! only master of first group needs to dump results if ( (qmmm_nml%vsolv < 2) .or. (worldrank == 0) ) then #endif call close_dump_files #ifdef MPI end if #endif if (icnstph /= 0) & call cnstph_finalize ! --- write out final times, taking REMD into account #ifdef MPI if (rem .ne. 0) then nstlim_total = nstlim * numexchg else nstlim_total = nstlim end if #else nstlim_total = nstlim #endif if (imin == 0) & call print_ongoing_time_summary(nstlim_total,nstlim_total,dt,6) write(6,'(12(a))') '| Job began at ', initial_time(1:2), & ':', initial_time(3:4), ':', initial_time(5:10), ' on ',& initial_date(5:6), '/', initial_date(7:8), '/', initial_date(1:4) write(6,'(12(a))') '| Setup done at ', setup_end_time(1:2), & ':', setup_end_time(3:4), ':', setup_end_time(5:10), ' on ', & setup_end_date(5:6), '/',setup_end_date(7:8),'/',setup_end_date(1:4) write(6,'(12(a))') '| Run done at ', final_time(1:2), & ':', final_time(3:4), ':', final_time(5:10), ' on ', & final_date(5:6), '/', final_date(7:8), '/', final_date(1:4) call nwallclock( ncalls ) write(6, '(''|'',5x,''wallclock() was called'',I8,'' times'')') ncalls call amflsh(6) if (iesp > 0) then call esp(natom,x(lcrd),x(linddip)) end if end if call amflsh(6) ! --- dynamic memory deallocation: 999 continue if(qmmm_nml%ifqnt .and. qmmm_nml%qmtheory%EXTERN .and. master) then call qm2_extern_finalize() endif if (qmmm_nml%ifqnt .and. .not. qmmm_struct%qm_mm_first_call) then ! If first_call is still true, this thread never really ! called the QMMM routine. E.g. more threads than PIMD replicates call deallocate_qmmm(qmmm_nml, qmmm_struct, qmmm_vsolv, qm2_params) end if if (ipimd > 0) call pimd_finalize(ipimd) if (ineb > 0) call neb_finalize() if (idecomp > 0) then call deallocate_real_decomp() call deallocate_int_decomp() end if if (master .and. idecomp == 0) call deallocate_int_decomp() #ifdef MPI /* SOFT CORE */ if (ifsc /= 0) call cleanup_sc() if (ifmbar /= 0) call cleanup_mbar() #endif ! finalize LIE module if initiated above if ( ilrt /= 0 ) then call cleanup_linear_response(master) end if #ifdef RISMSANDER call rism_finalize() #endif if ( ifcr /= 0 ) call cr_cleanup() if (sebomd_obj%do_sebomd) then if (master) then call sebomd_close_files end if call cleanup_sebomd_arrays end if if (master .and. iwrap == 2) then deallocate(iwrap_mask_atoms, stat=ier) REQUIRE(ier == 0) end if call nblist_deallocate() call deallocate_stacks() if ((igb /= 0 .and. igb /= 10 .and. ipb == 0) .or. hybridgb > 0 .or. icnstph > 1) then call deallocate_gb() end if if (master) then if(igb == 10 .or. ipb /= 0) then call pb_free() end if end if deallocate(ih, stat = ier) REQUIRE(ier == 0) deallocate(ipairs, stat = ier) REQUIRE(ier == 0) deallocate(ix, stat = ier) REQUIRE(ier == 0) deallocate(x, stat = ier) REQUIRE(ier == 0) if(ntb > 0 .and. ifbox == 1 .and. ew_type == 0 .and. mpoltype == 0) & call deallocate_m1m2m3() call AMOEBA_deallocate ! GMS: -- Module molecule -- call deallocate_molecule() ! -------------------------- if (charmm_active) call charmm_deallocate_arrays() if (cmap_active) call deallocate_cmap_arrays() if (master.and.mdout /= 'stdout') close(6) return end subroutine sander !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+ Calculate the ElectroStatic Potential subroutine esp(natom,x,mom_ind) ! routine to calculate the ESP due to the induced moments (only) ! at the same spatial points as the reference QM. use constants, only : zero, BOHRS_TO_A, INV_AMBER_ELECTROSTATIC use file_io_dat implicit none integer natom _REAL_ x(3,*) _REAL_ mom_ind(3,*) # include "ew_mpole.h" integer dat_unit, new_unit, minus_new_unit parameter(dat_unit=30, new_unit=31, minus_new_unit=33) integer inat, nesp, idum _REAL_ xin, yin, zin integer jn, kn _REAL_ esp_qm, xb_esp, yb_esp, zb_esp _REAL_ x_esp, y_esp, z_esp _REAL_ e_x, e_y, e_z, e_q, esp_new _REAL_ dist, dist3 integer iptr call amopen(dat_unit,"esp.dat",'O','F','R') call amopen(new_unit,"esp.induced",owrite,'F','W') call amopen(minus_new_unit,"esp.qm-induced",owrite,'F','W') read (dat_unit,'(3i5)')inat,nesp,idum write(6,'(t2,''inat = '',i5)')inat write(6,'(t2,''nesp = '',i5)')nesp write(new_unit,'(2i5)')inat,nesp write(minus_new_unit,'(2i5)')inat,nesp if (inat /= natom) then write(6,'(t2,''natom mismatch with esp file'')') call mexit(6,1) end if do jn = 1,inat read (dat_unit,'(17x,3e16.0)')xin,yin,zin write(new_unit,'(17x,3e16.7)')xin,yin,zin write(minus_new_unit,'(17x,3e16.7)')xin,yin,zin end do do jn = 1,nesp e_x = zero e_y = zero e_z = zero e_q = zero read(dat_unit,'(1x,4e16.0)')esp_qm,xb_esp,yb_esp,zb_esp x_esp = xb_esp * BOHRS_TO_A y_esp = yb_esp * BOHRS_TO_A z_esp = zb_esp * BOHRS_TO_A do kn = 1,natom dist = (sqrt((x(1,kn)-x_esp)**2 + & (x(2,kn)-y_esp)**2 + & (x(3,kn)-z_esp)**2)) dist3 = dist**3 e_x = e_x - mom_ind(1,kn )*(x(1,kn)-x_esp)/dist3 e_y = e_y - mom_ind(2,kn )*(x(2,kn)-y_esp)/dist3 e_z = e_z - mom_ind(3,kn )*(x(3,kn)-z_esp)/dist3 end do e_x = e_x * BOHRS_TO_A * INV_AMBER_ELECTROSTATIC e_y = e_y * BOHRS_TO_A * INV_AMBER_ELECTROSTATIC e_z = e_z * BOHRS_TO_A * INV_AMBER_ELECTROSTATIC e_q = e_q * BOHRS_TO_A * INV_AMBER_ELECTROSTATIC esp_new = e_x + e_y + e_z write(new_unit, '(1x,4e16.7)')esp_new, & xb_esp,yb_esp,zb_esp write(minus_new_unit,'(1x,4e16.7)')esp_qm-esp_new, & xb_esp,yb_esp,zb_esp end do close(dat_unit) close(new_unit) close(minus_new_unit) return end subroutine esp