+++ /dev/null
- 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
projects / unres.git / blobdiff