X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M-newcorr%2FMD_A-MTS.F_safe;fp=source%2Funres%2Fsrc_MD-M-newcorr%2FMD_A-MTS.F_safe;h=db8058fd5a0a879c345dddf7a56257983a73ea67;hb=7308760ff07636ef6b1ee28d8c3a67a23c14b34b;hp=0000000000000000000000000000000000000000;hpb=9a54ab407f6d0d9d564d52763b3e2136450b9ffc;p=unres.git diff --git a/source/unres/src_MD-M-newcorr/MD_A-MTS.F_safe b/source/unres/src_MD-M-newcorr/MD_A-MTS.F_safe new file mode 100644 index 0000000..db8058f --- /dev/null +++ b/source/unres/src_MD-M-newcorr/MD_A-MTS.F_safe @@ -0,0 +1,2327 @@ + 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 ' + return + end +c------------------------------------------------------------------------------- + subroutine velverlet_step(itime) +c------------------------------------------------------------------------------- +c Perform a single velocity Verlet step; the time step can be rescaled if +c increments in accelerations exceed the threshold +c------------------------------------------------------------------------------- + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' +#ifdef MPI + include 'mpif.h' + integer ierror,ierrcode +#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' + include 'COMMON.MUCA' + double precision vcm(3),incr(3) + double precision cm(3),L(3) + integer ilen,count,rstcount + external ilen + character*50 tytul + integer maxcount_scale /20/ + common /gucio/ cm + double precision stochforcvec(MAXRES6) + common /stochcalc/ stochforcvec + integer itime + logical scale +c + scale=.true. + icount_scale=0 + if (lang.eq.1) then + call sddir_precalc + else if (lang.eq.2 .or. lang.eq.3) then +#ifndef LANG0 + call stochastic_force(stochforcvec) +#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 + endif + itime_scal=0 + do while (scale) + icount_scale=icount_scale+1 + if (icount_scale.gt.maxcount_scale) then + write (iout,*) + & "ERROR: too many attempts at scaling down the time step. ", + & "amax=",amax,"epdrift=",epdrift, + & "damax=",damax,"edriftmax=",edriftmax, + & "d_time=",d_time + call flush(iout) +#ifdef MPI + call MPI_Abort(MPI_COMM_WORLD,IERROR,IERRCODE) +#endif + stop + endif +c First step of the velocity Verlet algorithm + if (lang.eq.2) then +#ifndef LANG0 + call sd_verlet1 +#endif + else if (lang.eq.3) then +#ifndef LANG0 + call sd_verlet1_ciccotti +#endif + else if (lang.eq.1) then + call sddir_verlet1 + else + call verlet1 + endif +c Build the chain from the newly calculated coordinates + call chainbuild_cart + if (rattle) call rattle1 + if (ntwe.ne.0) then + if (large.and. mod(itime,ntwe).eq.0) then + write (iout,*) "Cartesian and internal coordinates: step 1" + call cartprint + call intout + write (iout,*) "dC" + do i=0,nres + write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3), + & (dc(j,i+nres),j=1,3) + enddo + write (iout,*) "Accelerations" + do i=0,nres + write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3), + & (d_a(j,i+nres),j=1,3) + enddo + write (iout,*) "Velocities, step 1" + do i=0,nres + write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3), + & (d_t(j,i+nres),j=1,3) + enddo + endif + endif +#ifdef MPI + tt0 = MPI_Wtime() +#else + tt0 = tcpu() +#endif +c Calculate energy and forces + call zerograd + call etotal(potEcomp) + potE=potEcomp(0)-potEcomp(20) + call cartgrad +c Get the new accelerations + call lagrangian +#ifdef MPI + t_enegrad=t_enegrad+MPI_Wtime()-tt0 +#else + t_enegrad=t_enegrad+tcpu()-tt0 +#endif +c Determine maximum acceleration and scale down the timestep if needed + call max_accel + amax=amax/(itime_scal+1)**2 + call predict_edrift(epdrift) + if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then +c Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step + scale=.true. + ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax)) + & /dlog(2.0d0)+1 + itime_scal=itime_scal+ifac_time +c fac_time=dmin1(damax/amax,0.5d0) + fac_time=0.5d0**ifac_time + d_time=d_time*fac_time + if (lang.eq.2 .or. lang.eq.3) then +#ifndef LANG0 +c write (iout,*) "Calling sd_verlet_setup: 1" +c Rescale the stochastic forces and recalculate or restore +c the matrices of tinker integrator + if (itime_scal.gt.maxflag_stoch) then + if (large) write (iout,'(a,i5,a)') + & "Calculate matrices for stochastic step;", + & " itime_scal ",itime_scal + if (lang.eq.2) then + call sd_verlet_p_setup + else + call sd_verlet_ciccotti_setup + endif + write (iout,'(2a,i3,a,i3,1h.)') + & "Warning: cannot store matrices for stochastic", + & " integration because the index",itime_scal, + & " is greater than",maxflag_stoch + write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct", + & " integration Langevin algorithm for better efficiency." + else if (flag_stoch(itime_scal)) then + if (large) write (iout,'(a,i5,a,l1)') + & "Restore matrices for stochastic step; itime_scal ", + & itime_scal," flag ",flag_stoch(itime_scal) + do i=1,dimen + do j=1,dimen + pfric_mat(i,j)=pfric0_mat(i,j,itime_scal) + afric_mat(i,j)=afric0_mat(i,j,itime_scal) + vfric_mat(i,j)=vfric0_mat(i,j,itime_scal) + prand_mat(i,j)=prand0_mat(i,j,itime_scal) + vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal) + vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal) + enddo + enddo + else + if (large) write (iout,'(2a,i5,a,l1)') + & "Calculate & store matrices for stochastic step;", + & " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal) + if (lang.eq.2) then + call sd_verlet_p_setup + else + call sd_verlet_ciccotti_setup + endif + flag_stoch(ifac_time)=.true. + do i=1,dimen + do j=1,dimen + pfric0_mat(i,j,itime_scal)=pfric_mat(i,j) + afric0_mat(i,j,itime_scal)=afric_mat(i,j) + vfric0_mat(i,j,itime_scal)=vfric_mat(i,j) + prand0_mat(i,j,itime_scal)=prand_mat(i,j) + vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j) + vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j) + enddo + enddo + endif + fac_time=1.0d0/dsqrt(fac_time) + do i=1,dimen + stochforcvec(i)=fac_time*stochforcvec(i) + enddo +#endif + else if (lang.eq.1) then +c Rescale the accelerations due to stochastic forces + fac_time=1.0d0/dsqrt(fac_time) + do i=1,dimen + d_as_work(i)=d_as_work(i)*fac_time + enddo + endif + if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)') + & "itime",itime," Timestep scaled down to ", + & d_time," ifac_time",ifac_time," itime_scal",itime_scal + else +c Second step of the velocity Verlet algorithm + if (lang.eq.2) then +#ifndef LANG0 + call sd_verlet2 +#endif + else if (lang.eq.3) then +#ifndef LANG0 + call sd_verlet2_ciccotti +#endif + else if (lang.eq.1) then + call sddir_verlet2 + else + call verlet2 + endif + if (rattle) call rattle2 + totT=totT+d_time + if (d_time.ne.d_time0) then + d_time=d_time0 +#ifndef LANG0 + if (lang.eq.2 .or. lang.eq.3) then + if (large) write (iout,'(a)') + & "Restore original matrices for stochastic step" +c write (iout,*) "Calling sd_verlet_setup: 2" +c Restore the matrices of tinker integrator if the time step has been restored + do i=1,dimen + do j=1,dimen + pfric_mat(i,j)=pfric0_mat(i,j,0) + afric_mat(i,j)=afric0_mat(i,j,0) + vfric_mat(i,j)=vfric0_mat(i,j,0) + prand_mat(i,j)=prand0_mat(i,j,0) + vrand_mat1(i,j)=vrand0_mat1(i,j,0) + vrand_mat2(i,j)=vrand0_mat2(i,j,0) + enddo + enddo + endif +#endif + endif + scale=.false. + endif + enddo +c Calculate the kinetic and the total energy and the kinetic temperature + call kinetic(EK) + totE=EK+potE +c diagnostics +c call kinetic1(EK1) +c write (iout,*) "step",itime," EK",EK," EK1",EK1 +c end diagnostics +c Couple the system to Berendsen bath if needed + if (tbf .and. lang.eq.0) then + call verlet_bath + endif + kinetic_T=2.0d0/(dimen3*Rb)*EK +c Backup the coordinates, velocities, and accelerations + do i=0,2*nres + do j=1,3 + dc_old(j,i)=dc(j,i) + d_t_old(j,i)=d_t(j,i) + d_a_old(j,i)=d_a(j,i) + enddo + enddo + if (ntwe.ne.0) then + if (mod(itime,ntwe).eq.0 .and. large) then + write (iout,*) "Velocities, step 2" + do i=0,nres + write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3), + & (d_t(j,i+nres),j=1,3) + enddo + endif + endif + return + end +c------------------------------------------------------------------------------- + subroutine RESPA_step(itime) +c------------------------------------------------------------------------------- +c Perform a single RESPA step. +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 energia_short(0:n_ene), + & energia_long(0:n_ene) + double precision cm(3),L(3),vcm(3),incr(3) + double precision dc_old0(3,0:maxres2),d_t_old0(3,0:maxres2), + & d_a_old0(3,0:maxres2) + integer ilen,count,rstcount + external ilen + character*50 tytul + integer maxcount_scale /10/ + common /gucio/ cm + double precision stochforcvec(MAXRES6) + common /stochcalc/ stochforcvec + integer itime + logical scale + common /cipiszcze/ itt + itt=itime + if (ntwe.ne.0) then + if (large.and. mod(itime,ntwe).eq.0) then + write (iout,*) "***************** RESPA itime",itime + write (iout,*) "Cartesian and internal coordinates: step 0" +c call cartprint + call pdbout(0.0d0,"cipiszcze",iout) + call intout + write (iout,*) "Accelerations from long-range forces" + do i=0,nres + write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3), + & (d_a(j,i+nres),j=1,3) + enddo + write (iout,*) "Velocities, step 0" + do i=0,nres + write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3), + & (d_t(j,i+nres),j=1,3) + enddo + endif + endif +c +c Perform the initial RESPA step (increment velocities) +c write (iout,*) "*********************** RESPA ini" + call RESPA_vel + if (ntwe.ne.0) then + if (mod(itime,ntwe).eq.0 .and. large) then + write (iout,*) "Velocities, end" + do i=0,nres + write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3), + & (d_t(j,i+nres),j=1,3) + enddo + endif + endif +c Compute the short-range forces +#ifdef MPI + tt0 =MPI_Wtime() +#else + tt0 = tcpu() +#endif + call zerograd + call etotal_short(energia_short) + call cartgrad + call lagrangian + if (ntwe.ne.0) then + if (large.and. mod(itime,ntwe).eq.0) then + write (iout,*) "energia_short",energia_short(0) + write (iout,*) "Accelerations from short-range forces" + do i=0,nres + write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3), + & (d_a(j,i+nres),j=1,3) + enddo + endif + endif +#ifdef MPI + t_enegrad=t_enegrad+MPI_Wtime()-tt0 +#else + t_enegrad=t_enegrad+tcpu()-tt0 +#endif + do i=0,2*nres + do j=1,3 + dc_old(j,i)=dc(j,i) + d_t_old(j,i)=d_t(j,i) + d_a_old(j,i)=d_a(j,i) + enddo + enddo +c 6/30/08 A-MTS: attempt at increasing the split number + do i=0,2*nres + do j=1,3 + dc_old0(j,i)=dc_old(j,i) + d_t_old0(j,i)=d_t_old(j,i) + d_a_old0(j,i)=d_a_old(j,i) + enddo + enddo + if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2 + if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0 +c + scale=.true. + d_time0=d_time + do while (scale) + + scale=.false. +c write (iout,*) "itime",itime," ntime_split",ntime_split +c Split the time step + d_time=d_time0/ntime_split +c Perform the short-range RESPA steps (velocity Verlet increments of +c positions and velocities using short-range forces) +c write (iout,*) "*********************** RESPA split" + do itsplit=1,ntime_split + if (lang.eq.1) then + call sddir_precalc + else if (lang.eq.2 .or. lang.eq.3) then +#ifndef LANG0 + call stochastic_force(stochforcvec) +#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 + endif +c First step of the velocity Verlet algorithm + if (lang.eq.2) then +#ifndef LANG0 + call sd_verlet1 +#endif + else if (lang.eq.3) then +#ifndef LANG0 + call sd_verlet1_ciccotti +#endif + else if (lang.eq.1) then + call sddir_verlet1 + else + call verlet1 + endif +c Build the chain from the newly calculated coordinates + call chainbuild_cart + if (rattle) call rattle1 + if (ntwe.ne.0) then + if (large.and. mod(itime,ntwe).eq.0) then + write (iout,*) "***** ITSPLIT",itsplit + write (iout,*) "Cartesian and internal coordinates: step 1" + call pdbout(0.0d0,"cipiszcze",iout) +c call cartprint + call intout + write (iout,*) "Velocities, step 1" + do i=0,nres + write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3), + & (d_t(j,i+nres),j=1,3) + enddo + endif + endif +#ifdef MPI + tt0 = MPI_Wtime() +#else + tt0 = tcpu() +#endif +c Calculate energy and forces + call zerograd + call etotal_short(energia_short) + call cartgrad +c Get the new accelerations + call lagrangian + if (ntwe.ne.0) then + if (large.and. mod(itime,ntwe).eq.0) then + write (iout,*)"energia_short",energia_short(0) + write (iout,*) "Accelerations from short-range forces" + do i=0,nres + write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3), + & (d_a(j,i+nres),j=1,3) + enddo + endif + endif +c 6/30/08 A-MTS +c Determine maximum acceleration and scale down the timestep if needed + call max_accel + amax=amax/ntime_split**2 + call predict_edrift(epdrift) + if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0) + & write (iout,*) "amax",amax," damax",damax, + & " epdrift",epdrift," epdriftmax",epdriftmax +c Exit loop and try with increased split number if the change of +c acceleration is too big + if (amax.gt.damax .or. epdrift.gt.edriftmax) then + if (ntime_split.lt.maxtime_split) then + scale=.true. + ntime_split=ntime_split*2 + do i=0,2*nres + do j=1,3 + dc_old(j,i)=dc_old0(j,i) + d_t_old(j,i)=d_t_old0(j,i) + d_a_old(j,i)=d_a_old0(j,i) + enddo + enddo + write (iout,*) "acceleration/energy drift too large",amax, + & epdrift," split increased to ",ntime_split," itime",itime, + & " itsplit",itsplit + exit + else + write (iout,*) + & "Uh-hu. Bumpy landscape. Maximum splitting number", + & maxtime_split, + & " already reached!!! Trying to carry on!" + endif + endif +#ifdef MPI + t_enegrad=t_enegrad+MPI_Wtime()-tt0 +#else + t_enegrad=t_enegrad+tcpu()-tt0 +#endif +c Second step of the velocity Verlet algorithm + if (lang.eq.2) then +#ifndef LANG0 + call sd_verlet2 +#endif + else if (lang.eq.3) then +#ifndef LANG0 + call sd_verlet2_ciccotti +#endif + else if (lang.eq.1) then + call sddir_verlet2 + else + call verlet2 + endif + if (rattle) call rattle2 +c Backup the coordinates, velocities, and accelerations + do i=0,2*nres + do j=1,3 + dc_old(j,i)=dc(j,i) + d_t_old(j,i)=d_t(j,i) + d_a_old(j,i)=d_a(j,i) + enddo + enddo + enddo + + enddo ! while scale + +c Restore the time step + d_time=d_time0 +c Compute long-range forces +#ifdef MPI + tt0 =MPI_Wtime() +#else + tt0 = tcpu() +#endif + call zerograd + call etotal_long(energia_long) + call cartgrad + call lagrangian +#ifdef MPI + t_enegrad=t_enegrad+MPI_Wtime()-tt0 +#else + t_enegrad=t_enegrad+tcpu()-tt0 +#endif +c Compute accelerations from long-range forces + if (ntwe.ne.0) then + if (large.and. mod(itime,ntwe).eq.0) then + write (iout,*) "energia_long",energia_long(0) + write (iout,*) "Cartesian and internal coordinates: step 2" +c call cartprint + call pdbout(0.0d0,"cipiszcze",iout) + call intout + write (iout,*) "Accelerations from long-range forces" + do i=0,nres + write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3), + & (d_a(j,i+nres),j=1,3) + enddo + write (iout,*) "Velocities, step 2" + do i=0,nres + write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3), + & (d_t(j,i+nres),j=1,3) + enddo + endif + endif +c Compute the final RESPA step (increment velocities) +c write (iout,*) "*********************** RESPA fin" + call RESPA_vel +c Compute the complete potential energy + do i=0,n_ene + potEcomp(i)=energia_short(i)+energia_long(i) + enddo + potE=potEcomp(0)-potEcomp(20) +c potE=energia_short(0)+energia_long(0) + totT=totT+d_time +c Calculate the kinetic and the total energy and the kinetic temperature + call kinetic(EK) + totE=EK+potE +c Couple the system to Berendsen bath if needed + if (tbf .and. lang.eq.0) then + call verlet_bath + endif + kinetic_T=2.0d0/(dimen3*Rb)*EK +c Backup the coordinates, velocities, and accelerations + if (ntwe.ne.0) then + if (mod(itime,ntwe).eq.0 .and. large) then + write (iout,*) "Velocities, end" + do i=0,nres + write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3), + & (d_t(j,i+nres),j=1,3) + enddo + endif + endif + return + end +c--------------------------------------------------------------------- + subroutine RESPA_vel +c First and last RESPA step (incrementing velocities using long-range +c forces). + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.CONTROL' + include 'COMMON.VAR' + include 'COMMON.MD' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.GEO' + include 'COMMON.LOCAL' + include 'COMMON.INTERACT' + include 'COMMON.IOUNITS' + include 'COMMON.NAMES' + do j=1,3 + d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time + enddo + do i=nnt,nct-1 + do j=1,3 + d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time + enddo + enddo + do i=nnt,nct + if (itype(i).ne.10) then + inres=i+nres + do j=1,3 + d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time + enddo + endif + enddo + return + end +c----------------------------------------------------------------- + subroutine verlet1 +c Applying velocity Verlet algorithm - step 1 to coordinates + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.CONTROL' + include 'COMMON.VAR' + include 'COMMON.MD' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.GEO' + include 'COMMON.LOCAL' + include 'COMMON.INTERACT' + include 'COMMON.IOUNITS' + include 'COMMON.NAMES' + double precision adt,adt2 + +#ifdef DEBUG + write (iout,*) "VELVERLET1 START: DC" + do i=0,nres + write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3), + & (dc(j,i+nres),j=1,3) + enddo +#endif + do j=1,3 + adt=d_a_old(j,0)*d_time + adt2=0.5d0*adt + dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time + d_t_new(j,0)=d_t_old(j,0)+adt2 + d_t(j,0)=d_t_old(j,0)+adt + enddo + do i=nnt,nct-1 + do j=1,3 + adt=d_a_old(j,i)*d_time + adt2=0.5d0*adt + dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time + d_t_new(j,i)=d_t_old(j,i)+adt2 + d_t(j,i)=d_t_old(j,i)+adt + enddo + enddo + do i=nnt,nct + if (itype(i).ne.10) then + inres=i+nres + do j=1,3 + adt=d_a_old(j,inres)*d_time + adt2=0.5d0*adt + dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time + d_t_new(j,inres)=d_t_old(j,inres)+adt2 + d_t(j,inres)=d_t_old(j,inres)+adt + enddo + endif + enddo +#ifdef DEBUG + write (iout,*) "VELVERLET1 END: DC" + do i=0,nres + write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3), + & (dc(j,i+nres),j=1,3) + enddo +#endif + return + end +c--------------------------------------------------------------------- + subroutine verlet2 +c Step 2 of the velocity Verlet algorithm: update velocities + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.CONTROL' + include 'COMMON.VAR' + include 'COMMON.MD' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.GEO' + include 'COMMON.LOCAL' + include 'COMMON.INTERACT' + include 'COMMON.IOUNITS' + include 'COMMON.NAMES' + do j=1,3 + d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time + enddo + do i=nnt,nct-1 + do j=1,3 + d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time + enddo + enddo + do i=nnt,nct + if (itype(i).ne.10) then + inres=i+nres + do j=1,3 + d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time + enddo + endif + enddo + return + end +c----------------------------------------------------------------- + subroutine sddir_precalc +c Applying velocity Verlet algorithm - step 1 to coordinates + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + 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' + double precision stochforcvec(MAXRES6) + common /stochcalc/ stochforcvec +c +c Compute friction and stochastic forces +c + call friction_force + call stochastic_force(stochforcvec) +c +c Compute the acceleration due to friction forces (d_af_work) and stochastic +c forces (d_as_work) +c + call ginv_mult(fric_work, d_af_work) + call ginv_mult(stochforcvec, d_as_work) + return + end +c--------------------------------------------------------------------- + subroutine sddir_verlet1 +c Applying velocity Verlet algorithm - step 1 to velocities + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + 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' +c Revised 3/31/05 AL: correlation between random contributions to +c position and velocity increments included. + double precision sqrt13 /0.57735026918962576451d0/ ! 1/sqrt(3) + double precision adt,adt2 +c +c Add the contribution from BOTH friction and stochastic force to the +c coordinates, but ONLY the contribution from the friction forces to velocities +c + do j=1,3 + adt=(d_a_old(j,0)+d_af_work(j))*d_time + adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time + dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time + d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt + d_t(j,0)=d_t_old(j,0)+adt + enddo + ind=3 + do i=nnt,nct-1 + do j=1,3 + adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time + adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time + dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time + d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt + d_t(j,i)=d_t_old(j,i)+adt + enddo + ind=ind+3 + enddo + do i=nnt,nct + if (itype(i).ne.10) then + inres=i+nres + do j=1,3 + adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time + adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time + dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time + d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt + d_t(j,inres)=d_t_old(j,inres)+adt + enddo + ind=ind+3 + endif + enddo + return + end +c--------------------------------------------------------------------- + subroutine sddir_verlet2 +c Calculating the adjusted velocities for accelerations + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + 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' + double precision stochforcvec(MAXRES6),d_as_work1(MAXRES6) + double precision cos60 /0.5d0/, sin60 /0.86602540378443864676d0/ +c Revised 3/31/05 AL: correlation between random contributions to +c position and velocity increments included. +c The correlation coefficients are calculated at low-friction limit. +c Also, friction forces are now not calculated with new velocities. + +c call friction_force + call stochastic_force(stochforcvec) +c +c Compute the acceleration due to friction forces (d_af_work) and stochastic +c forces (d_as_work) +c + call ginv_mult(stochforcvec, d_as_work1) + +c +c Update velocities +c + do j=1,3 + d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j)) + & +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time + enddo + ind=3 + do i=nnt,nct-1 + do j=1,3 + d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j)) + & +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time + enddo + ind=ind+3 + enddo + do i=nnt,nct + if (itype(i).ne.10) then + inres=i+nres + do j=1,3 + d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres) + & +d_af_work(ind+j))+sin60*d_as_work(ind+j) + & +cos60*d_as_work1(ind+j))*d_time + enddo + ind=ind+3 + endif + enddo + return + end +c--------------------------------------------------------------------- + subroutine max_accel +c +c Find the maximum difference in the accelerations of the the sites +c at the beginning and the end of the time step. +c + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.CONTROL' + include 'COMMON.VAR' + include 'COMMON.MD' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.GEO' + include 'COMMON.LOCAL' + include 'COMMON.INTERACT' + include 'COMMON.IOUNITS' + double precision aux(3),accel(3),accel_old(3),dacc + do j=1,3 +c aux(j)=d_a(j,0)-d_a_old(j,0) + accel_old(j)=d_a_old(j,0) + accel(j)=d_a(j,0) + enddo + amax=0.0d0 + do i=nnt,nct +c Backbone + if (i.lt.nct) then +c 7/3/08 changed to asymmetric difference + do j=1,3 +c accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i)) + accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i) + accel(j)=accel(j)+0.5d0*d_a(j,i) +c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j)) + if (dabs(accel(j)).gt.dabs(accel_old(j))) then + dacc=dabs(accel(j)-accel_old(j)) + if (dacc.gt.amax) amax=dacc + endif + enddo + endif + enddo +c Side chains + do j=1,3 +c accel(j)=aux(j) + accel_old(j)=d_a_old(j,0) + accel(j)=d_a(j,0) + enddo + if (nnt.eq.2) then + do j=1,3 + accel_old(j)=accel_old(j)+d_a_old(j,1) + accel(j)=accel(j)+d_a(j,1) + enddo + endif + do i=nnt,nct + if (itype(i).ne.10) then + do j=1,3 +c accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres) + accel_old(j)=accel_old(j)+d_a_old(j,i+nres) + accel(j)=accel(j)+d_a(j,i+nres) + enddo + endif + do j=1,3 +c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j)) + if (dabs(accel(j)).gt.dabs(accel_old(j))) then + dacc=dabs(accel(j)-accel_old(j)) + if (dacc.gt.amax) amax=dacc + endif + enddo + do j=1,3 + accel_old(j)=accel_old(j)+d_a_old(j,i) + accel(j)=accel(j)+d_a(j,i) +c aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i) + enddo + enddo + return + end +c--------------------------------------------------------------------- + subroutine predict_edrift(epdrift) +c +c Predict the drift of the potential energy +c + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.CONTROL' + include 'COMMON.VAR' + include 'COMMON.MD' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.GEO' + include 'COMMON.LOCAL' + include 'COMMON.INTERACT' + include 'COMMON.IOUNITS' + include 'COMMON.MUCA' + double precision epdrift,epdriftij +c Drift of the potential energy + epdrift=0.0d0 + do i=nnt,nct +c Backbone + if (i.lt.nct) then + do j=1,3 + epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i)) + if (lmuca) epdriftij=epdriftij*factor +c write (iout,*) "back",i,j,epdriftij + if (epdriftij.gt.epdrift) epdrift=epdriftij + enddo + endif +c Side chains + if (itype(i).ne.10) then + do j=1,3 + epdriftij= + & dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i)) + if (lmuca) epdriftij=epdriftij*factor +c write (iout,*) "side",i,j,epdriftij + if (epdriftij.gt.epdrift) epdrift=epdriftij + enddo + endif + enddo + epdrift=0.5d0*epdrift*d_time*d_time +c write (iout,*) "epdrift",epdrift + return + end +c----------------------------------------------------------------------- + subroutine verlet_bath +c +c Coupling to the thermostat by using the Berendsen algorithm +c + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.CONTROL' + include 'COMMON.VAR' + include 'COMMON.MD' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.GEO' + include 'COMMON.LOCAL' + include 'COMMON.INTERACT' + include 'COMMON.IOUNITS' + include 'COMMON.NAMES' + double precision T_half,fact +c + T_half=2.0d0/(dimen3*Rb)*EK + fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0)) +c write(iout,*) "T_half", T_half +c write(iout,*) "EK", EK +c write(iout,*) "fact", fact + do j=1,3 + d_t(j,0)=fact*d_t(j,0) + enddo + do i=nnt,nct-1 + do j=1,3 + d_t(j,i)=fact*d_t(j,i) + enddo + enddo + do i=nnt,nct + if (itype(i).ne.10) then + inres=i+nres + do j=1,3 + d_t(j,inres)=fact*d_t(j,inres) + enddo + endif + enddo + return + end +c--------------------------------------------------------- + subroutine init_MD +c Set up the initial conditions of a MD simulation + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' +#ifdef MP + include 'mpif.h' + character*16 form + 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.REMD' + real*8 energia_long(0:n_ene), + & energia_short(0:n_ene),vcm(3),incr(3) + double precision cm(3),L(3),xv,sigv,lowb,highb + double precision varia(maxvar) + character*256 qstr + integer ilen + external ilen + character*50 tytul + logical file_exist + common /gucio/ cm + d_time0=d_time +c write(iout,*) "d_time", d_time +c Compute the standard deviations of stochastic forces for Langevin dynamics +c if the friction coefficients do not depend on surface area + if (lang.gt.0 .and. .not.surfarea) then + do i=nnt,nct-1 + stdforcp(i)=stdfp*dsqrt(gamp) + enddo + do i=nnt,nct + stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamsc(itype(i))) + enddo + endif +c Open the pdb file for snapshotshots +#ifdef MPI + if(mdpdb) then + if (ilen(tmpdir).gt.0) + & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// + & liczba(:ilen(liczba))//".pdb") + open(ipdb, + & file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) + & //".pdb") + else +#ifdef NOXDR + if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) + & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// + & liczba(:ilen(liczba))//".x") + cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) + & //".x" +#else + if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) + & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// + & liczba(:ilen(liczba))//".cx") + cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) + & //".cx" +#endif + endif +#else + if(mdpdb) then + if (ilen(tmpdir).gt.0) + & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb") + open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb") + else + if (ilen(tmpdir).gt.0) + & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx") + cartname=prefix(:ilen(prefix))//"_MD.cx" + endif +#endif + if (usampl) then + write (qstr,'(256(1h ))') + ipos=1 + do i=1,nfrag + iq = qinfrag(i,iset)*10 + iw = wfrag(i,iset)/100 + if (iw.gt.0) then + if(me.eq.king.or..not.out1file) + & write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw + write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw + ipos=ipos+7 + endif + enddo + do i=1,npair + iq = qinpair(i,iset)*10 + iw = wpair(i,iset)/100 + if (iw.gt.0) then + if(me.eq.king.or..not.out1file) + & write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw + write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw + ipos=ipos+7 + endif + enddo +c pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb' +#ifdef NOXDR +c cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x' +#else +c cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx' +#endif +c statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat' + endif + icg=1 + if (rest) then + if (restart1file) then + if (me.eq.king) + & inquire(file=mremd_rst_name,exist=file_exist) + write (*,*) me," Before broadcast: file_exist",file_exist + call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM, + & IERR) + write (*,*) me," After broadcast: file_exist",file_exist +c inquire(file=mremd_rst_name,exist=file_exist) + if(me.eq.king.or..not.out1file) + & write(iout,*) "Initial state read by master and distributed" + else + if (ilen(tmpdir).gt.0) + & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_' + & //liczba(:ilen(liczba))//'.rst') + inquire(file=rest2name,exist=file_exist) + endif + if(file_exist) then + if(.not.restart1file) then + if(me.eq.king.or..not.out1file) + & write(iout,*) "Initial state will be read from file ", + & rest2name(:ilen(rest2name)) + call readrst + endif + call rescale_weights(t_bath) + else + if(me.eq.king.or..not.out1file)then + if (restart1file) then + write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)), + & " does not exist" + else + write(iout,*) "File ",rest2name(:ilen(rest2name)), + & " does not exist" + endif + write(iout,*) "Initial velocities randomly generated" + endif + call random_vel + totT=0.0d0 + endif + else +c Generate initial velocities + if(me.eq.king.or..not.out1file) + & write(iout,*) "Initial velocities randomly generated" + call random_vel + totT=0.0d0 + endif +c rest2name = prefix(:ilen(prefix))//'.rst' + if(me.eq.king.or..not.out1file)then + write(iout,*) "Initial backbone velocities" + do i=nnt,nct-1 + write(iout,*) (d_t(j,i),j=1,3) + enddo + write(iout,*) "Initial side-chain velocities" + do i=nnt,nct + write(iout,*) (d_t(j,i+nres),j=1,3) + enddo +c Zeroing the total angular momentum of the system + write(iout,*) "Calling the zero-angular + & momentum subroutine" + endif + call inertia_tensor +c Getting the potential energy and forces and velocities and accelerations + call vcm_vel(vcm) +c write (iout,*) "velocity of the center of the mass:" +c write (iout,*) (vcm(j),j=1,3) + do j=1,3 + d_t(j,0)=d_t(j,0)-vcm(j) + enddo +c Removing the velocity of the center of mass + call vcm_vel(vcm) + if(me.eq.king.or..not.out1file)then + write (iout,*) "vcm right after adjustment:" + write (iout,*) (vcm(j),j=1,3) + endif + if (.not.rest) then + call chainbuild + if(iranconf.ne.0) then + if (overlapsc) then + print *, 'Calling OVERLAP_SC' + call overlap_sc(fail) + endif + + if (searchsc) then + call sc_move(2,nres-1,10,1d10,nft_sc,etot) + print *,'SC_move',nft_sc,etot + if(me.eq.king.or..not.out1file) + & write(iout,*) 'SC_move',nft_sc,etot + endif + + if(dccart)then + print *, 'Calling MINIM_DC' + call minim_dc(etot,iretcode,nfun) + else + call geom_to_var(nvar,varia) + print *,'Calling MINIMIZE.' + call minimize(etot,varia,iretcode,nfun) + call var_to_geom(nvar,varia) + endif + if(me.eq.king.or..not.out1file) + & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun + endif + endif + call chainbuild_cart + call kinetic(EK) + if (tbf) then + call verlet_bath(EK) + endif + kinetic_T=2.0d0/(dimen3*Rb)*EK + if(me.eq.king.or..not.out1file)then + call cartprint + call intout + endif + call zerograd + call etotal(potEcomp) + potE=potEcomp(0) + call cartgrad + call lagrangian + call max_accel + if (amax*d_time .gt. dvmax) then + d_time=d_time*dvmax/amax + if(me.eq.king.or..not.out1file) write (iout,*) + & "Time step reduced to",d_time, + & " because of too large initial acceleration." + endif + if(me.eq.king.or..not.out1file)then + write(iout,*) "Potential energy and its components" + write(iout,*) (potEcomp(i),i=0,n_ene) + endif + potE=potEcomp(0)-potEcomp(20) + totE=EK+potE + itime=0 + if (ntwe.ne.0) call statout(itime) + if(me.eq.king.or..not.out1file) + & write (iout,*) "Initial:", + & " Kinetic energy",EK," potential energy",potE, + & " total energy",totE," maximum acceleration ", + & amax + if (large) then + write (iout,*) "Initial coordinates" + do i=1,nres + write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3), + & (c(j,i+nres),j=1,3) + enddo + write (iout,*) "Initial dC" + do i=0,nres + write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3), + & (dc(j,i+nres),j=1,3) + enddo + write (iout,*) "Initial velocities" + do i=0,nres + write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3), + & (d_t(j,i+nres),j=1,3) + enddo + write (iout,*) "Initial accelerations" + do i=0,nres +c write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3), + write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3), + & (d_a(j,i+nres),j=1,3) + enddo + endif + do i=0,2*nres + do j=1,3 + dc_old(j,i)=dc(j,i) + d_t_old(j,i)=d_t(j,i) + d_a_old(j,i)=d_a(j,i) + enddo +c write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3) + enddo + if (RESPA) then +#ifdef MPI + tt0 =MPI_Wtime() +#else + tt0 = tcpu() +#endif + call zerograd + call etotal_long(energia_long) + if(.not.out1file .and. large) then + write (iout,*) "energia_long",energia_long(0) + write (iout,*) "Initial slow-force accelerations" + do i=0,nres + write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3), + & (d_a(j,i+nres),j=1,3) + enddo + endif + call cartgrad + call lagrangian +#ifdef MPI + t_enegrad=t_enegrad+MPI_Wtime()-tt0 +#else + t_enegrad=t_enegrad+tcpu()-tt0 +#endif +c call etotal_short(energia_short) +c write (iout,*) "energia_long",energia_long(0), +c & " energia_short",energia_short(0), +c & " total",energia_long(0)+energia_short(0) + endif + return + end +c----------------------------------------------------------- + subroutine random_vel + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + 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 xv,sigv,lowb,highb +c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m +c First generate velocities in the eigenspace of the G matrix +c write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3 +c call flush(iout) + xv=0.0d0 + ii=0 + do i=1,dimen + do k=1,3 + ii=ii+1 + sigv=dsqrt((Rb*t_bath)/geigen(i)) + lowb=-5*sigv + highb=5*sigv + d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb) +c write (iout,*) "ii",ii," d_t_work_new",d_t_work_new(ii) + enddo + enddo +c diagnostics +c Ek1=0.0d0 +c ii=0 +c do i=1,dimen +c do k=1,3 +c ii=ii+1 +c Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2 +c enddo +c enddo +c write (iout,*) "Ek from eigenvectors",Ek1 +c end diagnostics +c Transform velocities to UNRES coordinate space + do k=0,2 + do i=1,dimen + ind=(i-1)*3+k+1 + d_t_work(ind)=0.0d0 + do j=1,dimen + d_t_work(ind)=d_t_work(ind) + & +Gvec(i,j)*d_t_work_new((j-1)*3+k+1) + enddo +c write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind) +c call flush(iout) + enddo + enddo +c Transfer to the d_t vector + do j=1,3 + d_t(j,0)=d_t_work(j) + enddo + ind=3 + do i=nnt,nct-1 + do j=1,3 + ind=ind+1 + d_t(j,i)=d_t_work(ind) + enddo + enddo + do i=nnt,nct + if (itype(i).ne.10) then + do j=1,3 + ind=ind+1 + d_t(j,i+nres)=d_t_work(ind) + enddo + endif + enddo +c call kinetic(EK) +c write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature", +c & 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1 +c call flush(iout) + return + end +#ifndef LANG0 +c----------------------------------------------------------- + subroutine sd_verlet_p_setup +c Sets up the parameters of stochastic Verlet algorithm + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' +#ifdef MPI + include 'mpif.h' +#endif + 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 emgdt(MAXRES6), + & pterm,vterm,rho,rhoc,vsig, + & pfric_vec(MAXRES6),vfric_vec(MAXRES6), + & afric_vec(MAXRES6),prand_vec(MAXRES6), + & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6) + logical lprn /.false./ + double precision zero /1.0d-8/, gdt_radius /0.05d0/ + double precision ktm +#ifdef MPI + tt0 = MPI_Wtime() +#else + tt0 = tcpu() +#endif +c +c AL 8/17/04 Code adapted from tinker +c +c Get the frictional and random terms for stochastic dynamics in the +c eigenspace of mass-scaled UNRES friction matrix +c + do i = 1, dimen + gdt = fricgam(i) * d_time +c +c Stochastic dynamics reduces to simple MD for zero friction +c + if (gdt .le. zero) then + pfric_vec(i) = 1.0d0 + vfric_vec(i) = d_time + afric_vec(i) = 0.5d0 * d_time * d_time + prand_vec(i) = 0.0d0 + vrand_vec1(i) = 0.0d0 + vrand_vec2(i) = 0.0d0 +c +c Analytical expressions when friction coefficient is large +c + else + if (gdt .ge. gdt_radius) then + egdt = dexp(-gdt) + pfric_vec(i) = egdt + vfric_vec(i) = (1.0d0-egdt) / fricgam(i) + afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i) + pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt + vterm = 1.0d0 - egdt**2 + rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm) +c +c Use series expansions when friction coefficient is small +c + else + gdt2 = gdt * gdt + gdt3 = gdt * gdt2 + gdt4 = gdt2 * gdt2 + gdt5 = gdt2 * gdt3 + gdt6 = gdt3 * gdt3 + gdt7 = gdt3 * gdt4 + gdt8 = gdt4 * gdt4 + gdt9 = gdt4 * gdt5 + afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0 + & - gdt5/120.0d0 + gdt6/720.0d0 + & - gdt7/5040.0d0 + gdt8/40320.0d0 + & - gdt9/362880.0d0) / fricgam(i)**2 + vfric_vec(i) = d_time - fricgam(i)*afric_vec(i) + pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i) + pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0 + & + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0 + & + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0 + & + 127.0d0*gdt9/90720.0d0 + vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0 + & - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0 + & - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0 + & - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0 + rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0 + & - 17.0d0*gdt2/1280.0d0 + & + 17.0d0*gdt3/6144.0d0 + & + 40967.0d0*gdt4/34406400.0d0 + & - 57203.0d0*gdt5/275251200.0d0 + & - 1429487.0d0*gdt6/13212057600.0d0) + end if +c +c Compute the scaling factors of random terms for the nonzero friction case +c + ktm = 0.5d0*d_time/fricgam(i) + psig = dsqrt(ktm*pterm) / fricgam(i) + vsig = dsqrt(ktm*vterm) + rhoc = dsqrt(1.0d0 - rho*rho) + prand_vec(i) = psig + vrand_vec1(i) = vsig * rho + vrand_vec2(i) = vsig * rhoc + end if + end do + if (lprn) then + write (iout,*) + & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,", + & " vrand_vec2" + do i=1,dimen + write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i), + & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i) + enddo + endif +c +c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables +c +#ifndef LANG0 + call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat) + call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat) + call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat) + call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat) + call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1) + call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2) +#endif +#ifdef MPI + t_sdsetup=t_sdsetup+MPI_Wtime() +#else + t_sdsetup=t_sdsetup+tcpu()-tt0 +#endif + return + end +c------------------------------------------------------------- + subroutine eigtransf1(n,ndim,ab,d,c) + implicit none + integer n,ndim + double precision ab(ndim,ndim,n),c(ndim,n),d(ndim) + integer i,j,k + do i=1,n + do j=1,n + c(i,j)=0.0d0 + do k=1,n + c(i,j)=c(i,j)+ab(k,j,i)*d(k) + enddo + enddo + enddo + return + end +c------------------------------------------------------------- + subroutine eigtransf(n,ndim,a,b,d,c) + implicit none + integer n,ndim + double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim) + integer i,j,k + do i=1,n + do j=1,n + c(i,j)=0.0d0 + do k=1,n + c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k) + enddo + enddo + enddo + return + end +c------------------------------------------------------------- + subroutine sd_verlet1 +c Applying stochastic velocity Verlet algorithm - step 1 to velocities + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + 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' + double precision stochforcvec(MAXRES6) + common /stochcalc/ stochforcvec + logical lprn /.false./ + +c write (iout,*) "dc_old" +c do i=0,nres +c write (iout,'(i5,3f10.5,5x,3f10.5)') +c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3) +c enddo + do j=1,3 + dc_work(j)=dc_old(j,0) + d_t_work(j)=d_t_old(j,0) + d_a_work(j)=d_a_old(j,0) + enddo + ind=3 + do i=nnt,nct-1 + do j=1,3 + dc_work(ind+j)=dc_old(j,i) + d_t_work(ind+j)=d_t_old(j,i) + d_a_work(ind+j)=d_a_old(j,i) + enddo + ind=ind+3 + enddo + do i=nnt,nct + if (itype(i).ne.10) then + do j=1,3 + dc_work(ind+j)=dc_old(j,i+nres) + d_t_work(ind+j)=d_t_old(j,i+nres) + d_a_work(ind+j)=d_a_old(j,i+nres) + enddo + ind=ind+3 + endif + enddo +#ifndef LANG0 + if (lprn) then + write (iout,*) + & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,", + & " vrand_mat2" + do i=1,dimen + do j=1,dimen + write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j), + & vfric_mat(i,j),afric_mat(i,j), + & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j) + enddo + enddo + endif + do i=1,dimen + ddt1=0.0d0 + ddt2=0.0d0 + do j=1,dimen + dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) + & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j) + ddt1=ddt1+pfric_mat(i,j)*d_t_work(j) + ddt2=ddt2+vfric_mat(i,j)*d_a_work(j) + enddo + d_t_work_new(i)=ddt1+0.5d0*ddt2 + d_t_work(i)=ddt1+ddt2 + enddo +#endif + do j=1,3 + dc(j,0)=dc_work(j) + d_t(j,0)=d_t_work(j) + enddo + ind=3 + do i=nnt,nct-1 + do j=1,3 + dc(j,i)=dc_work(ind+j) + d_t(j,i)=d_t_work(ind+j) + enddo + ind=ind+3 + enddo + do i=nnt,nct + if (itype(i).ne.10) then + inres=i+nres + do j=1,3 + dc(j,inres)=dc_work(ind+j) + d_t(j,inres)=d_t_work(ind+j) + enddo + ind=ind+3 + endif + enddo + return + end +c-------------------------------------------------------------------------- + subroutine sd_verlet2 +c Calculating the adjusted velocities for accelerations + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + 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' + double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6) + common /stochcalc/ stochforcvec +c +c Compute the stochastic forces which contribute to velocity change +c + call stochastic_force(stochforcvecV) + +#ifndef LANG0 + do i=1,dimen + ddt1=0.0d0 + ddt2=0.0d0 + do j=1,dimen + ddt1=ddt1+vfric_mat(i,j)*d_a_work(j) + ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+ + & vrand_mat2(i,j)*stochforcvecV(j) + enddo + d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2 + enddo +#endif + do j=1,3 + d_t(j,0)=d_t_work(j) + enddo + ind=3 + do i=nnt,nct-1 + do j=1,3 + d_t(j,i)=d_t_work(ind+j) + enddo + ind=ind+3 + enddo + do i=nnt,nct + if (itype(i).ne.10) then + inres=i+nres + do j=1,3 + d_t(j,inres)=d_t_work(ind+j) + enddo + ind=ind+3 + endif + enddo + return + end +c----------------------------------------------------------- + subroutine sd_verlet_ciccotti_setup +c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's +c version + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' +#ifdef MPI + include 'mpif.h' +#endif + 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 emgdt(MAXRES6), + & pterm,vterm,rho,rhoc,vsig, + & pfric_vec(MAXRES6),vfric_vec(MAXRES6), + & afric_vec(MAXRES6),prand_vec(MAXRES6), + & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6) + logical lprn /.false./ + double precision zero /1.0d-8/, gdt_radius /0.05d0/ + double precision ktm +#ifdef MPI + tt0 = MPI_Wtime() +#else + tt0 = tcpu() +#endif +c +c AL 8/17/04 Code adapted from tinker +c +c Get the frictional and random terms for stochastic dynamics in the +c eigenspace of mass-scaled UNRES friction matrix +c + do i = 1, dimen + write (iout,*) "i",i," fricgam",fricgam(i) + gdt = fricgam(i) * d_time +c +c Stochastic dynamics reduces to simple MD for zero friction +c + if (gdt .le. zero) then + pfric_vec(i) = 1.0d0 + vfric_vec(i) = d_time + afric_vec(i) = 0.5d0*d_time*d_time + prand_vec(i) = afric_vec(i) + vrand_vec2(i) = vfric_vec(i) +c +c Analytical expressions when friction coefficient is large +c + else + egdt = dexp(-gdt) + pfric_vec(i) = egdt + vfric_vec(i) = dexp(-0.5d0*gdt)*d_time + afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time + prand_vec(i) = afric_vec(i) + vrand_vec2(i) = vfric_vec(i) +c +c Compute the scaling factors of random terms for the nonzero friction case +c +c ktm = 0.5d0*d_time/fricgam(i) +c psig = dsqrt(ktm*pterm) / fricgam(i) +c vsig = dsqrt(ktm*vterm) +c prand_vec(i) = psig*afric_vec(i) +c vrand_vec2(i) = vsig*vfric_vec(i) + end if + end do + if (lprn) then + write (iout,*) + & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,", + & " vrand_vec2" + do i=1,dimen + write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i), + & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i) + enddo + endif +c +c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables +c + call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat) + call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat) + call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat) + call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat) + call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2) +#ifdef MPI + t_sdsetup=t_sdsetup+MPI_Wtime() +#else + t_sdsetup=t_sdsetup+tcpu()-tt0 +#endif + return + end +c------------------------------------------------------------- + subroutine sd_verlet1_ciccotti +c Applying stochastic velocity Verlet algorithm - step 1 to velocities + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' +#ifdef MPI + include 'mpif.h' +#endif + 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' + double precision stochforcvec(MAXRES6) + common /stochcalc/ stochforcvec + logical lprn /.false./ + +c write (iout,*) "dc_old" +c do i=0,nres +c write (iout,'(i5,3f10.5,5x,3f10.5)') +c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3) +c enddo + do j=1,3 + dc_work(j)=dc_old(j,0) + d_t_work(j)=d_t_old(j,0) + d_a_work(j)=d_a_old(j,0) + enddo + ind=3 + do i=nnt,nct-1 + do j=1,3 + dc_work(ind+j)=dc_old(j,i) + d_t_work(ind+j)=d_t_old(j,i) + d_a_work(ind+j)=d_a_old(j,i) + enddo + ind=ind+3 + enddo + do i=nnt,nct + if (itype(i).ne.10) then + do j=1,3 + dc_work(ind+j)=dc_old(j,i+nres) + d_t_work(ind+j)=d_t_old(j,i+nres) + d_a_work(ind+j)=d_a_old(j,i+nres) + enddo + ind=ind+3 + endif + enddo + +#ifndef LANG0 + if (lprn) then + write (iout,*) + & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,", + & " vrand_mat2" + do i=1,dimen + do j=1,dimen + write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j), + & vfric_mat(i,j),afric_mat(i,j), + & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j) + enddo + enddo + endif + do i=1,dimen + ddt1=0.0d0 + ddt2=0.0d0 + do j=1,dimen + dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) + & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j) + ddt1=ddt1+pfric_mat(i,j)*d_t_work(j) + ddt2=ddt2+vfric_mat(i,j)*d_a_work(j) + enddo + d_t_work_new(i)=ddt1+0.5d0*ddt2 + d_t_work(i)=ddt1+ddt2 + enddo +#endif + do j=1,3 + dc(j,0)=dc_work(j) + d_t(j,0)=d_t_work(j) + enddo + ind=3 + do i=nnt,nct-1 + do j=1,3 + dc(j,i)=dc_work(ind+j) + d_t(j,i)=d_t_work(ind+j) + enddo + ind=ind+3 + enddo + do i=nnt,nct + if (itype(i).ne.10) then + inres=i+nres + do j=1,3 + dc(j,inres)=dc_work(ind+j) + d_t(j,inres)=d_t_work(ind+j) + enddo + ind=ind+3 + endif + enddo + return + end +c-------------------------------------------------------------------------- + subroutine sd_verlet2_ciccotti +c Calculating the adjusted velocities for accelerations + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + 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' + double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6) + common /stochcalc/ stochforcvec +c +c Compute the stochastic forces which contribute to velocity change +c + call stochastic_force(stochforcvecV) +#ifndef LANG0 + do i=1,dimen + ddt1=0.0d0 + ddt2=0.0d0 + do j=1,dimen + + ddt1=ddt1+vfric_mat(i,j)*d_a_work(j) +c ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j) + ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j) + enddo + d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2 + enddo +#endif + do j=1,3 + d_t(j,0)=d_t_work(j) + enddo + ind=3 + do i=nnt,nct-1 + do j=1,3 + d_t(j,i)=d_t_work(ind+j) + enddo + ind=ind+3 + enddo + do i=nnt,nct + if (itype(i).ne.10) then + inres=i+nres + do j=1,3 + d_t(j,inres)=d_t_work(ind+j) + enddo + ind=ind+3 + endif + enddo + return + end +#endif