subroutine MD c------------------------------------------------ c The driver for molecular dynamics subroutines c------------------------------------------------ implicit real*8 (a-h,o-z) include 'DIMENSIONS' #ifdef MPI include "mpif.h" integer IERROR,ERRCODE #endif include 'COMMON.SETUP' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' #ifndef LANG0 include 'COMMON.LANGEVIN' #else include 'COMMON.LANGEVIN.lang0' #endif include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.NAMES' include 'COMMON.TIME1' double precision cm(3),L(3),vcm(3) #ifdef VOUT double precision v_work(maxres6),v_transf(maxres6) #endif integer ilen,rstcount external ilen character*50 tytul common /gucio/ cm integer itime c #ifdef MPI if (ilen(tmpdir).gt.0) & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_" & //liczba(:ilen(liczba))//'.rst') #else if (ilen(tmpdir).gt.0) & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"//'.rst') #endif t_MDsetup=0.0d0 t_langsetup=0.0d0 t_MD=0.0d0 t_enegrad=0.0d0 t_sdsetup=0.0d0 write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started" #ifdef MPI tt0=MPI_Wtime() #else tt0 = tcpu() #endif c Determine the inverse of the inertia matrix. call setup_MD_matrices c Initialize MD call init_MD #ifdef MPI t_MDsetup = MPI_Wtime()-tt0 #else t_MDsetup = tcpu()-tt0 #endif rstcount=0 c Entering the MD loop #ifdef MPI tt0 = MPI_Wtime() #else tt0 = tcpu() #endif if (lang.eq.2 .or. lang.eq.3) then #ifndef LANG0 call setup_fricmat if (lang.eq.2) then call sd_verlet_p_setup else call sd_verlet_ciccotti_setup endif do i=1,dimen do j=1,dimen pfric0_mat(i,j,0)=pfric_mat(i,j) afric0_mat(i,j,0)=afric_mat(i,j) vfric0_mat(i,j,0)=vfric_mat(i,j) prand0_mat(i,j,0)=prand_mat(i,j) vrand0_mat1(i,j,0)=vrand_mat1(i,j) vrand0_mat2(i,j,0)=vrand_mat2(i,j) enddo enddo flag_stoch(0)=.true. do i=1,maxflag_stoch flag_stoch(i)=.false. enddo #else write (iout,*) & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0" #ifdef MPI call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE) #endif stop #endif else if (lang.eq.1 .or. lang.eq.4) then call setup_fricmat endif #ifdef MPI t_langsetup=MPI_Wtime()-tt0 tt0=MPI_Wtime() #else t_langsetup=tcpu()-tt0 tt0=tcpu() #endif do itime=1,n_timestep rstcount=rstcount+1 if (lang.gt.0 .and. surfarea .and. & mod(itime,reset_fricmat).eq.0) then if (lang.eq.2 .or. lang.eq.3) then #ifndef LANG0 call setup_fricmat if (lang.eq.2) then call sd_verlet_p_setup else call sd_verlet_ciccotti_setup endif do i=1,dimen do j=1,dimen pfric0_mat(i,j,0)=pfric_mat(i,j) afric0_mat(i,j,0)=afric_mat(i,j) vfric0_mat(i,j,0)=vfric_mat(i,j) prand0_mat(i,j,0)=prand_mat(i,j) vrand0_mat1(i,j,0)=vrand_mat1(i,j) vrand0_mat2(i,j,0)=vrand_mat2(i,j) enddo enddo flag_stoch(0)=.true. do i=1,maxflag_stoch flag_stoch(i)=.false. enddo #endif else if (lang.eq.1 .or. lang.eq.4) then call setup_fricmat endif write (iout,'(a,i10)') & "Friction matrix reset based on surface area, itime",itime endif if (reset_vel .and. tbf .and. lang.eq.0 & .and. mod(itime,count_reset_vel).eq.0) then call random_vel write(iout,'(a,f20.2)') & "Velocities reset to random values, time",totT do i=0,2*nres do j=1,3 d_t_old(j,i)=d_t(j,i) enddo enddo endif if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then call inertia_tensor call vcm_vel(vcm) do j=1,3 d_t(j,0)=d_t(j,0)-vcm(j) enddo call kinetic(EK) kinetic_T=2.0d0/(dimen3*Rb)*EK scalfac=dsqrt(T_bath/kinetic_T) write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT do i=0,2*nres do j=1,3 d_t_old(j,i)=scalfac*d_t(j,i) enddo enddo endif if (lang.ne.4) then if (RESPA) then c Time-reversible RESPA algorithm c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992) call RESPA_step(itime) else c Variable time step algorithm. call velverlet_step(itime) endif else #ifdef BROWN call brown_step(itime) #else print *,"Brown dynamics not here!" #ifdef MPI call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE) #endif stop #endif endif if (ntwe.ne.0) then if (mod(itime,ntwe).eq.0) call statout(itime) #ifdef VOUT do j=1,3 v_work(j)=d_t(j,0) enddo ind=3 do i=nnt,nct-1 do j=1,3 ind=ind+1 v_work(ind)=d_t(j,i) enddo enddo do i=nnt,nct if (itype(i).ne.10) then do j=1,3 ind=ind+1 v_work(ind)=d_t(j,i+nres) enddo endif enddo write (66,'(80f10.5)') & ((d_t(j,i),j=1,3),i=0,nres-1),((d_t(j,i+nres),j=1,3),i=1,nres) do i=1,ind v_transf(i)=0.0d0 do j=1,ind v_transf(i)=v_transf(i)+gvec(j,i)*v_work(j) enddo v_transf(i)= v_transf(i)*dsqrt(geigen(i)) enddo write (67,'(80f10.5)') (v_transf(i),i=1,ind) #endif endif if (mod(itime,ntwx).eq.0) then write (tytul,'("time",f8.2)') totT if(mdpdb) then call pdbout(potE,tytul,ipdb) else call cartout(totT) endif endif if (rstcount.eq.1000.or.itime.eq.n_timestep) then open(irest2,file=rest2name,status='unknown') write(irest2,*) totT,EK,potE,totE,t_bath do i=1,2*nres write (irest2,'(3e15.5)') (d_t(j,i),j=1,3) enddo do i=1,2*nres write (irest2,'(3e15.5)') (dc(j,i),j=1,3) enddo close(irest2) rstcount=0 endif enddo #ifdef MPI t_MD=MPI_Wtime()-tt0 #else t_MD=tcpu()-tt0 #endif write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))') & ' Timing ', & 'MD calculations setup:',t_MDsetup, & 'Energy & gradient evaluation:',t_enegrad, & 'Stochastic MD setup:',t_langsetup, & 'Stochastic MD step setup:',t_sdsetup, & 'MD steps:',t_MD write (iout,'(/28(1h=),a25,27(1h=))') & ' End of MD calculation ' #ifdef TIMING_ENE write (iout,*) "time for etotal",t_etotal," elong",t_elong, & " eshort",t_eshort write (iout,*) "time_fric",time_fric," time_stoch",time_stoch, & " time_fricmatmult",time_fricmatmult," time_fsample ", & time_fsample #endif return end