c------------------------------------------------
c The driver for molecular dynamics subroutines
c------------------------------------------------
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
#ifdef MPI
include "mpif.h"
include 'COMMON.CONTROL'
include 'COMMON.VAR'
include 'COMMON.MD'
+#ifdef FIVEDIAG
+ include 'COMMON.LAGRANGE.5diag'
+#else
+ include 'COMMON.LAGRANGE'
+#endif
#ifndef LANG0
include 'COMMON.LANGEVIN'
#else
+#ifdef FIVEDIAG
+ include 'COMMON.LANGEVIN.lang0.5diag'
+#else
include 'COMMON.LANGEVIN.lang0'
#endif
+#endif
include 'COMMON.CHAIN'
include 'COMMON.DERIV'
include 'COMMON.GEO'
common /gucio/ cm
integer itime
logical ovrtim
+ integer i,j,icount_scale,itime_scal
+ integer nharp,iharp(4,maxres/3)
+ double precision scalfac
+ double precision tt0
c
#ifdef MPI
if (ilen(tmpdir).gt.0)
if (ilen(tmpdir).gt.0)
& call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"//'.rst')
#endif
+ write (iout,*) "MD lang",lang
t_MDsetup=0.0d0
t_langsetup=0.0d0
t_MD=0.0d0
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)
+ implicit none
include 'DIMENSIONS'
#ifdef MPI
include 'mpif.h'
- integer ierror,ierrcode
+ integer ierror,ierrcode,errcode
#endif
include 'COMMON.SETUP'
include 'COMMON.CONTROL'
include 'COMMON.VAR'
include 'COMMON.MD'
+#ifdef FIVEDIAG
+ include 'COMMON.LAGRANGE.5diag'
+#else
+ include 'COMMON.LAGRANGE'
+#endif
#ifndef LANG0
include 'COMMON.LANGEVIN'
#else
+#ifdef FIVEDIAG
+ include 'COMMON.LANGEVIN.lang0.5diag'
+#else
include 'COMMON.LANGEVIN.lang0'
#endif
+#endif
include 'COMMON.CHAIN'
include 'COMMON.DERIV'
include 'COMMON.GEO'
common /gucio/ cm
double precision stochforcvec(MAXRES6)
common /stochcalc/ stochforcvec
- integer itime
+ integer itime,icount_scale,itime_scal,ifac_time,i,j,itt
logical scale
+ double precision epdrift,fac_time
+ double precision tt0
c
scale=.true.
icount_scale=0
call zerograd
call etotal(potEcomp)
! AL 4/17/17: Reduce the steps if NaNs occurred.
- if (potEcomp(0).gt.0.99e20 .or. isnan(potEcomp(0))) then
+ if (potEcomp(0).gt.0.99e20 .or. isnan(potEcomp(0)).gt.0) then
d_time=d_time/2
cycle
endif
c-------------------------------------------------------------------------------
c Perform a single RESPA step.
c-------------------------------------------------------------------------------
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
#ifdef MPI
include 'mpif.h'
include 'COMMON.CONTROL'
include 'COMMON.VAR'
include 'COMMON.MD'
+#ifdef FIVEDIAG
+ include 'COMMON.LAGRANGE.5diag'
+#else
+ include 'COMMON.LAGRANGE'
+#endif
#ifndef LANG0
include 'COMMON.LANGEVIN'
#else
+#ifdef FIVEDIAG
+ include 'COMMON.LANGEVIN.lang0.5diag'
+#else
include 'COMMON.LANGEVIN.lang0'
#endif
+#endif
include 'COMMON.CHAIN'
include 'COMMON.DERIV'
include 'COMMON.GEO'
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 i,j
+ double precision fac_time
logical PRINT_AMTS_MSG /.false./
integer ilen,count,rstcount
external ilen
common /stochcalc/ stochforcvec
integer itime
logical scale
+ integer itt
common /cipiszcze/ itt
+ integer itsplit
+ double precision epdrift,epdriftmax
+ double precision tt0
itt=itime
if (ntwe.ne.0) then
if (large.and. mod(itime,ntwe).eq.0) then
write (iout,*) "Cartesian and internal coordinates: step 2"
c call cartprint
call pdbout(0.0d0,
- & cipiszcze ,iout)
+ & 'cipiszcze ',iout)
call intout
write (iout,*) "Accelerations from long-range forces"
do i=0,nres
if (ntwe.ne.0) then
if (large.and. mod(itime,ntwe).eq.0) then
call enerprint(potEcomp)
- write (iout,*) "potE",potD
+ write (iout,*) "potE",potE
endif
endif
c potE=energia_short(0)+energia_long(0)
subroutine RESPA_vel
c First and last RESPA step (incrementing velocities using long-range
c forces).
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
include 'COMMON.CONTROL'
include 'COMMON.VAR'
include 'COMMON.MD'
+#ifdef FIVEDIAG
+ include 'COMMON.LAGRANGE.5diag'
+#else
+ include 'COMMON.LAGRANGE'
+#endif
include 'COMMON.CHAIN'
include 'COMMON.DERIV'
include 'COMMON.GEO'
include 'COMMON.INTERACT'
include 'COMMON.IOUNITS'
include 'COMMON.NAMES'
+ integer i,j,inres
do j=1,3
d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
enddo
c-----------------------------------------------------------------
subroutine verlet1
c Applying velocity Verlet algorithm - step 1 to coordinates
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
include 'COMMON.CONTROL'
include 'COMMON.VAR'
include 'COMMON.MD'
+#ifdef FIVEDIAG
+ include 'COMMON.LAGRANGE.5diag'
+#else
+ include 'COMMON.LAGRANGE'
+#endif
include 'COMMON.CHAIN'
include 'COMMON.DERIV'
include 'COMMON.GEO'
include 'COMMON.IOUNITS'
include 'COMMON.NAMES'
double precision adt,adt2
-
+ integer i,j,inres
#ifdef DEBUG
write (iout,*) "VELVERLET1 START: DC"
do i=0,nres
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 i=nnt,nct-1
C SPYTAC ADAMA
C do i=0,nres
+#ifdef DEBUG
+ write (iout,*) "i",i," d_a_old",(d_a_old(j,i),j=1,3)
+#endif
do j=1,3
adt=d_a_old(j,i)*d_time
adt2=0.5d0*adt
c---------------------------------------------------------------------
subroutine verlet2
c Step 2 of the velocity Verlet algorithm: update velocities
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
include 'COMMON.CONTROL'
include 'COMMON.VAR'
include 'COMMON.MD'
+#ifdef FIVEDIAG
+ include 'COMMON.LAGRANGE.5diag'
+#else
+ include 'COMMON.LAGRANGE'
+#endif
include 'COMMON.CHAIN'
include 'COMMON.DERIV'
include 'COMMON.GEO'
include 'COMMON.INTERACT'
include 'COMMON.IOUNITS'
include 'COMMON.NAMES'
+ integer i,j,inres
do j=1,3
d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
enddo
c-----------------------------------------------------------------
subroutine sddir_precalc
c Applying velocity Verlet algorithm - step 1 to coordinates
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
#ifdef MPI
include 'mpif.h'
include 'COMMON.CONTROL'
include 'COMMON.VAR'
include 'COMMON.MD'
+#ifdef FIVEDIAG
+ include 'COMMON.LAGRANGE.5diag'
+#else
+ include 'COMMON.LAGRANGE'
+#endif
#ifndef LANG0
include 'COMMON.LANGEVIN'
#else
+#ifdef FIVEDIAG
+ include 'COMMON.LANGEVIN.lang0.5diag'
+#else
include 'COMMON.LANGEVIN.lang0'
#endif
+#endif
include 'COMMON.CHAIN'
include 'COMMON.DERIV'
include 'COMMON.GEO'
include 'COMMON.IOUNITS'
include 'COMMON.NAMES'
include 'COMMON.TIME1'
+ double precision time00
double precision stochforcvec(MAXRES6)
common /stochcalc/ stochforcvec
+ integer i
c
c Compute friction and stochastic forces
c
time00=tcpu()
#endif
call friction_force
+c write (iout,*) "After friction_force"
#ifdef MPI
time_fric=time_fric+MPI_Wtime()-time00
time00=MPI_Wtime()
time00=tcpu()
#endif
call stochastic_force(stochforcvec)
+c write (iout,*) "After stochastic_force"
#ifdef MPI
time_stoch=time_stoch+MPI_Wtime()-time00
#else
c Compute the acceleration due to friction forces (d_af_work) and stochastic
c forces (d_as_work)
c
+#ifdef FIVEDIAG
+c write (iout,*) "friction accelerations"
+ call fivediaginv_mult(dimen,fric_work, d_af_work)
+c write (iout,*) "stochastic acceleratios"
+ call fivediaginv_mult(dimen,stochforcvec, d_as_work)
+c write (iout,*) "Leaving sddir_precalc"
+#else
call ginv_mult(fric_work, d_af_work)
call ginv_mult(stochforcvec, d_as_work)
+#endif
+#ifdef DEBUG
+ write (iout,*) "d_af_work"
+ write (iout,'(3f10.5)') (d_af_work(i),i=1,dimen3)
+ write (iout,*) "d_as_work"
+ write (iout,'(3f10.5)') (d_as_work(i),i=1,dimen3)
+#endif
return
end
c---------------------------------------------------------------------
subroutine sddir_verlet1
c Applying velocity Verlet algorithm - step 1 to velocities
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
include 'COMMON.CONTROL'
include 'COMMON.VAR'
include 'COMMON.MD'
+#ifdef FIVEDIAG
+ include 'COMMON.LAGRANGE.5diag'
+#else
+ include 'COMMON.LAGRANGE'
+#endif
#ifndef LANG0
include 'COMMON.LANGEVIN'
#else
+#ifdef FIVEDIAG
+ include 'COMMON.LANGEVIN.lang0.5diag'
+#else
include 'COMMON.LANGEVIN.lang0'
#endif
+#endif
include 'COMMON.CHAIN'
include 'COMMON.DERIV'
include 'COMMON.GEO'
c position and velocity increments included.
double precision sqrt13 /0.57735026918962576451d0/ ! 1/sqrt(3)
double precision adt,adt2
+ integer i,j,ind,inres
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
d_t(j,0)=d_t_old(j,0)+adt
enddo
ind=3
- do i=nnt,nct-1
+ 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
c---------------------------------------------------------------------
subroutine sddir_verlet2
c Calculating the adjusted velocities for accelerations
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
include 'COMMON.CONTROL'
include 'COMMON.VAR'
include 'COMMON.MD'
+#ifdef FIVEDIAG
+ include 'COMMON.LAGRANGE.5diag'
+#else
+ include 'COMMON.LAGRANGE'
+#endif
#ifndef LANG0
include 'COMMON.LANGEVIN'
#else
+#ifdef FIVEDIAG
+ include 'COMMON.LANGEVIN.lang0.5diag'
+#else
include 'COMMON.LANGEVIN.lang0'
#endif
+#endif
include 'COMMON.CHAIN'
include 'COMMON.DERIV'
include 'COMMON.GEO'
include 'COMMON.NAMES'
double precision stochforcvec(MAXRES6),d_as_work1(MAXRES6)
double precision cos60 /0.5d0/, sin60 /0.86602540378443864676d0/
+ integer i,j,inres,ind
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 Compute the acceleration due to friction forces (d_af_work) and stochastic
c forces (d_as_work)
c
+#ifdef FIVEDIAG
+ call fivediaginv_mult(maxres6,stochforcvec, d_as_work1)
+#else
call ginv_mult(stochforcvec, d_as_work1)
-
+#endif
c
c Update velocities
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)
+ implicit none
include 'DIMENSIONS'
include 'COMMON.CONTROL'
include 'COMMON.VAR'
include 'COMMON.MD'
+#ifdef FIVEDIAG
+ include 'COMMON.LAGRANGE.5diag'
+#else
+ include 'COMMON.LAGRANGE'
+#endif
include 'COMMON.CHAIN'
include 'COMMON.DERIV'
include 'COMMON.GEO'
include 'COMMON.LOCAL'
include 'COMMON.INTERACT'
include 'COMMON.IOUNITS'
+ integer i,j
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)
c
c Predict the drift of the potential energy
c
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
include 'COMMON.CONTROL'
include 'COMMON.VAR'
include 'COMMON.MD'
+#ifdef FIVEDIAG
+ include 'COMMON.LAGRANGE.5diag'
+#else
+ include 'COMMON.LAGRANGE'
+#endif
include 'COMMON.CHAIN'
include 'COMMON.DERIV'
include 'COMMON.GEO'
include 'COMMON.IOUNITS'
include 'COMMON.MUCA'
double precision epdrift,epdriftij
+ integer i,j
c Drift of the potential energy
epdrift=0.0d0
do i=nnt,nct
c
c Coupling to the thermostat by using the Berendsen algorithm
c
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
include 'COMMON.CONTROL'
include 'COMMON.VAR'
include 'COMMON.MD'
+#ifdef FIVEDIAG
+ include 'COMMON.LAGRANGE.5diag'
+#else
+ include 'COMMON.LAGRANGE'
+#endif
+#ifdef LANG0
+#ifdef FIVEDIAG
+ include 'COMMON.LANGEVIN.lang0.5diag'
+#else
+ include 'COMMON.LANGEVIN.lang0'
+#endif
+#else
+ include 'COMMON.LANGEVIN'
+#endif
include 'COMMON.CHAIN'
include 'COMMON.DERIV'
include 'COMMON.GEO'
include 'COMMON.IOUNITS'
include 'COMMON.NAMES'
double precision T_half,fact
-c
+ integer i,j ,inres
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---------------------------------------------------------
subroutine init_MD
c Set up the initial conditions of a MD simulation
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
#ifdef MP
include 'mpif.h'
character*16 form
- integer IERROR,ERRCODE
+ integer IERROR,ERRCODE,error_msg,ierr,ierrcode
#endif
include 'COMMON.SETUP'
include 'COMMON.CONTROL'
include 'COMMON.VAR'
include 'COMMON.MD'
+#ifdef FIVEDIAG
+ include 'COMMON.LAGRANGE.5diag'
+#else
+ include 'COMMON.LAGRANGE'
+#endif
+ include 'COMMON.QRESTR'
#ifndef LANG0
include 'COMMON.LANGEVIN'
#else
+#ifdef FIVEDIAG
+ include 'COMMON.LANGEVIN.lang0.5diag'
+#else
+#ifdef FIVEDIAG
+ include 'COMMON.LANGEVIN.lang0.5diag'
+#else
include 'COMMON.LANGEVIN.lang0'
#endif
+#endif
+#endif
include 'COMMON.CHAIN'
include 'COMMON.DERIV'
include 'COMMON.GEO'
include 'COMMON.IOUNITS'
include 'COMMON.NAMES'
include 'COMMON.REMD'
+ include 'COMMON.TIME1'
+#ifdef LBFGS
+ character*9 status
+ integer niter
+ common /lbfgstat/ status,niter,nfun
+#endif
+ integer n_model_try,list_model_try(max_template),k
+ double precision tt0
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
character*50 tytul
logical file_exist
common /gucio/ cm
+ integer i,ipos,iq,iw,j,iranmin,nft_sc,iretcode,nfun,itrial,itmp,
+ & i_model,itime
+ integer iran_num
+ double precision etot
+ logical fail
write (iout,*) "init_MD INDPDB",indpdb
d_time0=d_time
c write(iout,*) "d_time", d_time
endif
call flush(iout)
write (iout,*) "init_MD before initial structure REST ",rest
- if (.not.rest) then
+ if (.not.rest) then
+ 122 continue
if (iranconf.ne.0) then
c 8/22/17 AL Loop to produce a low-energy random conformation
do iranmin=1,10
44 continue
else if (preminim) then
if (start_from_model) then
- i_model=iran_num(1,constr_homology)
- write (iout,*) 'starting from model ',i_model
- do i=1,2*nres
- do j=1,3
- c(j,i)=chomo(j,i,i_model)
+ n_model_try=0
+ do while (fail .and. n_model_try.lt.constr_homology)
+ do
+ i_model=iran_num(1,constr_homology)
+ do k=1,n_model_try
+ if (i_model.eq.list_model_try(k)) exit
+ enddo
+ if (k.gt.n_model_try) exit
enddo
- enddo
- call int_from_cart(.true.,.false.)
- call sc_loc_geom(.false.)
- do i=1,nres-1
- do j=1,3
- dc(j,i)=c(j,i+1)-c(j,i)
- dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
+ n_model_try=n_model_try+1
+ list_model_try(n_model_try)=i_model
+ write (iout,*) 'starting from model ',i_model
+ do i=1,2*nres
+ do j=1,3
+ c(j,i)=chomo(j,i,i_model)
+ enddo
enddo
- enddo
- do i=2,nres-1
- do j=1,3
- dc(j,i+nres)=c(j,i+nres)-c(j,i)
- dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
+ call int_from_cart(.true.,.false.)
+ call sc_loc_geom(.false.)
+ do i=1,nres-1
+ do j=1,3
+ dc(j,i)=c(j,i+1)-c(j,i)
+ dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
+ enddo
enddo
- enddo
- endif
+ do i=2,nres-1
+ do j=1,3
+ dc(j,i+nres)=c(j,i+nres)-c(j,i)
+ dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
+ enddo
+ enddo
+ if (me.eq.king.or..not.out1file) then
+ write (iout,*) "Energies before removing overlaps"
+ call etotal(energia(0))
+ call enerprint(energia(0))
+ endif
! Remove SC overlaps if requested
- if (overlapsc) then
- write (iout,*) 'Calling OVERLAP_SC'
- call overlap_sc(fail)
- endif
+ if (overlapsc) then
+ write (iout,*) 'Calling OVERLAP_SC'
+ call overlap_sc(fail)
+ if (fail) then
+ write (iout,*)
+ & "Failed to remove overlap from model",i_model
+ cycle
+ endif
+ endif
+#ifdef SEARCHSC
+ if (me.eq.king.or..not.out1file) then
+ write (iout,*) "Energies after removing overlaps"
+ call etotal(energia(0))
+ call enerprint(energia(0))
+ endif
! Search for better SC rotamers if requested
- 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
+ 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
+ call etotal(energia(0))
+#endif
+ enddo
+ if (n_model_try.gt.constr_homology) then
+ write (iout,*)
+ & "All models have irreparable overlaps. Trying randoms starts."
+ iranconf=1
+ goto 122
+ endif
endif
- call etotal(energia(0))
C 8/22/17 AL Minimize initial structure
if (dccart) then
if (me.eq.king.or..not.out1file) write(iout,*)
- & 'Minimizing initial PDB structure: Calling MINIM_DC'
+ & 'Minimizing initial PDB structure: Calling MINIM_DC'
call minim_dc(etot,iretcode,nfun)
else
call geom_to_var(nvar,varia)
if(me.eq.king.or..not.out1file) write (iout,*)
- & 'Minimizing initial PDB structure: Calling MINIMIZE.'
+ & 'Minimizing initial PDB structure: Calling MINIMIZE.'
call minimize(etot,varia,iretcode,nfun)
call var_to_geom(nvar,varia)
- endif
- if (me.eq.king.or..not.out1file)
+#ifdef LBFGS
+ if (me.eq.king.or..not.out1file)
+ & write(iout,*) 'LBFGS return code is ',status,' eval ',nfun
+ if(me.eq.king.or..not.out1file)
+ & write(iout,*) 'LBFGS return code is ',status,' eval ',nfun
+#else
+ if (me.eq.king.or..not.out1file)
& write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
- if(me.eq.king.or..not.out1file)
+ if(me.eq.king.or..not.out1file)
& write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
+#endif
+ endif
endif
- endif
+ endif ! .not. rest
call chainbuild_cart
call kinetic(EK)
if (tbf) then
call verlet_bath
- endif
+ endif
kinetic_T=2.0d0/(dimen3*Rb)*EK
if(me.eq.king.or..not.out1file)then
call cartprint
end
c-----------------------------------------------------------
subroutine random_vel
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
include 'COMMON.CONTROL'
include 'COMMON.VAR'
include 'COMMON.MD'
+#ifdef FIVEDIAG
+ include 'COMMON.LAGRANGE.5diag'
+#else
+ include 'COMMON.LAGRANGE'
+#endif
#ifndef LANG0
include 'COMMON.LANGEVIN'
#else
+#ifdef FIVEDIAG
+ include 'COMMON.LANGEVIN.lang0.5diag'
+#else
include 'COMMON.LANGEVIN.lang0'
#endif
+#endif
include 'COMMON.CHAIN'
include 'COMMON.DERIV'
include 'COMMON.GEO'
include 'COMMON.IOUNITS'
include 'COMMON.NAMES'
include 'COMMON.TIME1'
- double precision xv,sigv,lowb,highb,vec_afm(3)
+ double precision xv,sigv,lowb,highb,vec_afm(3),Ek1,Ek2,Ek3,aux
+ integer i,ii,j,k,l,ind
+ double precision anorm_distr
+ logical lprn /.true./
+#ifdef FIVEDIAG
+ integer ichain,n,innt,inct,ibeg,ierr
+ double precision work(8*maxres6)
+ integer iwork(maxres6)
+ double precision Ghalf(mmaxres2_chain),Geigen(maxres2_chain),
+ & Gvec(maxres2_chain,maxres2_chain)
+ common /przechowalnia/Ghalf,Geigen,Gvec
+#ifdef DEBUG
+ double precision inertia(maxres2_chain,maxres2_chain)
+#endif
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)
+#ifdef DEBUG
+ write (iout,*) "Random_vel, fivediag"
+#endif
+ d_t=0.0d0
+ Ek2=0.0d0
+ EK=0.0d0
+ Ek3=0.0d0
+ do ichain=1,nchain
+ ind=0
+ ghalf=0.0d0
+ n=dimen_chain(ichain)
+ innt=iposd_chain(ichain)
+ inct=innt+n-1
+#ifdef DEBUG
+ write (iout,*) "Chain",ichain," n",n," start",innt
+ do i=innt,inct
+ if (i.lt.inct-1) then
+ write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i),DU1orig(i),
+ & DU2orig(i)
+ else if (i.eq.inct-1) then
+ write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i),DU1orig(i)
+ else
+ write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i)
+ endif
+ enddo
+#endif
+ ghalf(ind+1)=dmorig(innt)
+ ghalf(ind+2)=du1orig(innt)
+ ghalf(ind+3)=dmorig(innt+1)
+ ind=ind+3
+ do i=3,n
+ ind=ind+i-3
+c write (iout,*) "i",i," ind",ind," indu2",innt+i-2,
+c & " indu1",innt+i-1," indm",innt+i
+ ghalf(ind+1)=du2orig(innt-1+i-2)
+ ghalf(ind+2)=du1orig(innt-1+i-1)
+ ghalf(ind+3)=dmorig(innt-1+i)
+c write (iout,'(3(a,i2,1x))') "DU2",innt-1+i-2,
+c & "DU1",innt-1+i-1,"DM ",innt-1+i
+ ind=ind+3
+ enddo
+#ifdef DEBUG
+ ind=0
+ do i=1,n
+ do j=1,i
+ ind=ind+1
+ inertia(i,j)=ghalf(ind)
+ inertia(j,i)=ghalf(ind)
+ enddo
+ enddo
+#endif
+#ifdef DEBUG
+ write (iout,*) "Chain ",ichain," ind",ind," dim",n*(n+1)/2
+ write (iout,*) "Five-diagonal inertia matrix, lower triangle"
+ call matoutr(n,ghalf)
+#endif
+ call gldiag(maxres2_chain,n,n,Ghalf,work,Geigen,Gvec,ierr,iwork)
+ if (large) then
+ write (iout,'(//a,i3)')
+ & "Eigenvectors and eigenvalues of the G matrix chain",ichain
+ call eigout(n,n,maxres2_chain,maxres2_chain,Gvec,Geigen)
+ endif
+#ifdef DIAGCHECK
+c check diagonalization
+ do i=1,n
+ do j=1,n
+ aux=0.0d0
+ do k=1,n
+ do l=1,n
+ aux=aux+gvec(k,i)*gvec(l,j)*inertia(k,l)
+ enddo
+ enddo
+ if (i.eq.j) then
+ write (iout,*) i,j,aux,geigen(i)
+ else
+ write (iout,*) i,j,aux
+ endif
+ enddo
+ enddo
+#endif
+ xv=0.0d0
+ ii=0
+ do i=1,n
+ 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)
+ EK=EK+0.5d0*geigen(i)*d_t_work_new(ii)**2
+c write (iout,*) "i",i," ii",ii," geigen",geigen(i),
+c & " d_t_work_new",d_t_work_new(ii)
+ enddo
+ enddo
+ do k=1,3
+ do i=1,n
+ ind=(i-1)*3+k
+ d_t_work(ind)=0.0d0
+ do j=1,n
+ d_t_work(ind)=d_t_work(ind)
+ & +Gvec(i,j)*d_t_work_new((j-1)*3+k)
+ enddo
+c write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
+c call flush(iout)
+ enddo
+ enddo
+#ifdef DEBUG
+ aux=0.0d0
+ do k=1,3
+ do i=1,n
+ do j=1,n
+ aux=aux+inertia(i,j)*d_t_work(3*(i-1)+k)*d_t_work(3*(j-1)+k)
+ enddo
+ enddo
+ enddo
+ Ek3=Ek3+aux/2
+#endif
+c Transfer to the d_t vector
+ innt=chain_border(1,ichain)
+ inct=chain_border(2,ichain)
+ ind=0
+c write (iout,*) "ichain",ichain," innt",innt," inct",inct
+ do i=innt,inct
+ do j=1,3
+ ind=ind+1
+ d_t(j,i)=d_t_work(ind)
+ enddo
+ if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ do j=1,3
+ ind=ind+1
+ d_t(j,i+nres)=d_t_work(ind)
+ enddo
+ endif
+ enddo
+ enddo
+ if (large) then
+ write (iout,*)
+ write (iout,*) "Random velocities in the Calpha,SC space"
+ do i=1,nres
+ write (iout,'(a3,1h(,i5,1h),3f10.5,3x,3f10.5)')
+ & restyp(itype(i)),i,(d_t(j,i),j=1,3),(d_t(j,i+nres),j=1,3)
+ enddo
+ endif
+ call kinetic_CASC(Ek1)
+!
+! Transform the velocities to virtual-bond space
+!
+#define WLOS
+#ifdef WLOS
+ do i=1,nres
+ if (itype(i).eq.10 .or. itype(i).eq.ntyp1) then
+ do j=1,3
+ d_t(j,i)=d_t(j,i+1)-d_t(j,i)
+ enddo
+ else
+ do j=1,3
+ d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
+ d_t(j,i)=d_t(j,i+1)-d_t(j,i)
+ enddo
+ end if
+ enddo
+ d_t(:,nct)=0.0d0
+c d_a(:,0)=d_a(:,1)
+c d_a(:,1)=0.0d0
+c write (iout,*) "Shifting accelerations"
+ do ichain=1,nchain
+c write (iout,*) "ichain",chain_border1(1,ichain)-1,
+c & chain_border1(1,ichain)
+ d_t(:,chain_border1(1,ichain)-1)=d_t(:,chain_border1(1,ichain))
+ d_t(:,chain_border1(1,ichain))=0.0d0
+ enddo
+c write (iout,*) "Adding accelerations"
+ do ichain=2,nchain
+c write (iout,*) "chain",ichain,chain_border1(1,ichain)-1,
+c & chain_border(2,ichain-1)
+ d_t(:,chain_border1(1,ichain)-1)=
+ & d_t(:,chain_border1(1,ichain)-1)+d_t(:,chain_border(2,ichain-1))
+ d_t(:,chain_border(2,ichain-1))=0.0d0
+ enddo
+ do ichain=2,nchain
+ write (iout,*) "chain",ichain,chain_border1(1,ichain)-1,
+ & chain_border(2,ichain-1)
+ d_t(:,chain_border1(1,ichain)-1)=
+ & d_t(:,chain_border1(1,ichain)-1)+d_t(:,chain_border(2,ichain-1))
+ d_t(:,chain_border(2,ichain-1))=0.0d0
+ enddo
+#else
+ ibeg=0
+c do j=1,3
+c d_t(j,0)=d_t(j,nnt)
+c enddo
+ do ichain=1,nchain
+ innt=chain_border(1,ichain)
+ inct=chain_border(2,ichain)
+c write (iout,*) "ichain",ichain," innt",innt," inct",inct
+c write (iout,*) "ibeg",ibeg
+ do j=1,3
+ d_t(j,ibeg)=d_t(j,innt)
+ enddo
+ ibeg=inct+1
+ do i=innt,inct
+ if (iabs(itype(i).eq.10)) then
+c write (iout,*) "i",i,(d_t(j,i),j=1,3),(d_t(j,i+1),j=1,3)
+ do j=1,3
+ d_t(j,i)=d_t(j,i+1)-d_t(j,i)
+ enddo
+ else
+ do j=1,3
+ d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
+ d_t(j,i)=d_t(j,i+1)-d_t(j,i)
+ enddo
+ end if
+ enddo
+ enddo
+#endif
+ if (large) then
+ write (iout,*)
+ write (iout,*)
+ & "Random velocities in the virtual-bond-vector space"
+ write (iout,'(3hORG,1h(,i5,1h),3f10.5)') 0,(d_t(j,0),j=1,3)
+ do i=1,nres
+ write (iout,'(a3,1h(,i5,1h),3f10.5,3x,3f10.5)')
+ & restyp(itype(i)),i,(d_t(j,i),j=1,3),(d_t(j,i+nres),j=1,3)
+ enddo
+ write (iout,*)
+ write (iout,*) "Kinetic energy from inertia matrix eigenvalues",
+ & Ek
+ write (iout,*)
+ & "Kinetic temperatures from inertia matrix eigenvalues",
+ & 2*Ek/(3*dimen*Rb)
+#ifdef DEBUG
+ write (iout,*) "Kinetic energy from inertia matrix",Ek3
+ write (iout,*) "Kinetic temperatures from inertia",
+ & 2*Ek3/(3*dimen*Rb)
+#endif
+ write (iout,*) "Kinetic energy from velocities in CA-SC space",
+ & Ek1
+ write (iout,*)
+ & "Kinetic temperatures from velovities in CA-SC space",
+ & 2*Ek1/(3*dimen*Rb)
+ call kinetic(Ek1)
+ write (iout,*)
+ & "Kinetic energy from virtual-bond-vector velocities",Ek1
+ write (iout,*)
+ & "Kinetic temperature from virtual-bond-vector velocities ",
+ & 2*Ek1/(dimen3*Rb)
+ endif
+#else
xv=0.0d0
ii=0
do i=1,dimen
c write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
c & 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
c call flush(iout)
+#endif
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)
+ implicit none
include 'DIMENSIONS'
#ifdef MPI
include 'mpif.h'
#ifndef LANG0
include 'COMMON.LANGEVIN'
#else
+#ifdef FIVEDIAG
+ include 'COMMON.LANGEVIN.lang0.5diag'
+#else
include 'COMMON.LANGEVIN.lang0'
#endif
+#endif
include 'COMMON.CHAIN'
include 'COMMON.DERIV'
include 'COMMON.GEO'
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
c-------------------------------------------------------------
subroutine sd_verlet1
c Applying stochastic velocity Verlet algorithm - step 1 to velocities
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
include 'COMMON.CONTROL'
include 'COMMON.VAR'
include 'COMMON.MD'
+#ifdef FIVEDIAG
+ include 'COMMON.LAGRANGE.5diag'
+#else
+ include 'COMMON.LAGRANGE'
+#endif
#ifndef LANG0
include 'COMMON.LANGEVIN'
#else
+#ifdef FIVEDIAG
+ include 'COMMON.LANGEVIN.lang0.5diag'
+#else
include 'COMMON.LANGEVIN.lang0'
#endif
+#endif
include 'COMMON.CHAIN'
include 'COMMON.DERIV'
include 'COMMON.GEO'
double precision stochforcvec(MAXRES6)
common /stochcalc/ stochforcvec
logical lprn /.false./
+ integer i,j,ind,inres
c write (iout,*) "dc_old"
c do i=0,nres
c--------------------------------------------------------------------------
subroutine sd_verlet2
c Calculating the adjusted velocities for accelerations
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
include 'COMMON.CONTROL'
include 'COMMON.VAR'
include 'COMMON.MD'
+#ifdef FIVEDIAG
+ include 'COMMON.LAGRANGE.5diag'
+#else
+ include 'COMMON.LAGRANGE'
+#endif
#ifndef LANG0
include 'COMMON.LANGEVIN'
#else
+#ifdef FIVEDIAG
+ include 'COMMON.LANGEVIN.lang0.5diag'
+#else
include 'COMMON.LANGEVIN.lang0'
#endif
+#endif
include 'COMMON.CHAIN'
include 'COMMON.DERIV'
include 'COMMON.GEO'
include 'COMMON.NAMES'
double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
common /stochcalc/ stochforcvec
+ integer i,j,ind,inres
c
c Compute the stochastic forces which contribute to velocity change
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)
+ implicit none
include 'DIMENSIONS'
#ifdef MPI
include 'mpif.h'
include 'COMMON.CONTROL'
include 'COMMON.VAR'
include 'COMMON.MD'
-#ifndef LANG0
- include 'COMMON.LANGEVIN'
+#ifdef FIVEDIAG
+ include 'COMMON.LAGRANGE.5diag'
#else
- include 'COMMON.LANGEVIN.lang0'
+ include 'COMMON.LAGRANGE'
#endif
+ include 'COMMON.LANGEVIN'
include 'COMMON.CHAIN'
include 'COMMON.DERIV'
include 'COMMON.GEO'
logical lprn /.false./
double precision zero /1.0d-8/, gdt_radius /0.05d0/
double precision ktm
+ integer i
#ifdef MPI
tt0 = MPI_Wtime()
#else
c-------------------------------------------------------------
subroutine sd_verlet1_ciccotti
c Applying stochastic velocity Verlet algorithm - step 1 to velocities
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
#ifdef MPI
include 'mpif.h'
include 'COMMON.CONTROL'
include 'COMMON.VAR'
include 'COMMON.MD'
-#ifndef LANG0
- include 'COMMON.LANGEVIN'
+#ifdef FIVEDIAG
+ include 'COMMON.LAGRANGE.5diag'
#else
- include 'COMMON.LANGEVIN.lang0'
+ include 'COMMON.LAGRANGE'
#endif
+ include 'COMMON.LANGEVIN'
include 'COMMON.CHAIN'
include 'COMMON.DERIV'
include 'COMMON.GEO'
double precision stochforcvec(MAXRES6)
common /stochcalc/ stochforcvec
logical lprn /.false./
+ integer i,j
c write (iout,*) "dc_old"
c do i=0,nres
c--------------------------------------------------------------------------
subroutine sd_verlet2_ciccotti
c Calculating the adjusted velocities for accelerations
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
include 'COMMON.CONTROL'
include 'COMMON.VAR'
include 'COMMON.MD'
-#ifndef LANG0
- include 'COMMON.LANGEVIN'
+#ifdef FIVEDIAG
+ include 'COMMON.LAGRANGE.5diag'
#else
- include 'COMMON.LANGEVIN.lang0'
+ include 'COMMON.LAGRANGE'
#endif
+ include 'COMMON.LANGEVIN'
include 'COMMON.CHAIN'
include 'COMMON.DERIV'
include 'COMMON.GEO'
include 'COMMON.NAMES'
double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
common /stochcalc/ stochforcvec
+ integer i,j
c
c Compute the stochastic forces which contribute to velocity change
c
projects / unres.git / blobdiff