--- /dev/null
+ module MDyn
+!-----------------------------------------------------------------------------
+ use io_units
+ use names
+ use math
+ use md_calc
+ use geometry_data
+ use io_base
+ use geometry
+ use energy
+ use MD_data
+ use REMD
+
+ implicit none
+!-----------------------------------------------------------------------------
+! common.MD
+! common /mdgrad/ in module.energy
+! common /back_constr/ in module.energy
+! common /qmeas/ in module.energy
+! common /mdpar/
+! common /MDcalc/
+! common /lagrange/
+ real(kind=8),dimension(:),allocatable :: d_t_work,&
+ d_t_work_new,d_af_work,d_as_work,kinetic_force !(MAXRES6)
+ real(kind=8),dimension(:,:),allocatable :: d_t_new,&
+ d_a_old,d_a_short!,d_a !(3,0:MAXRES2)
+! real(kind=8),dimension(:),allocatable :: d_a_work !(6*MAXRES)
+! real(kind=8),dimension(:,:),allocatable :: Gmat,Ginv,A,&
+! Gsqrp,Gsqrm,Gvec !(maxres2,maxres2)
+! real(kind=8),dimension(:),allocatable :: Geigen !(maxres2)
+! integer :: dimen,dimen1,dimen3
+! integer :: lang,count_reset_moment,count_reset_vel
+! logical :: reset_moment,reset_vel,rattle,RESPA
+! common /inertia/
+! common /langevin/
+! real(kind=8) :: rwat,etawat,stdfp,cPoise
+! real(kind=8),dimension(:),allocatable :: gamsc !(ntyp1)
+! real(kind=8),dimension(:),allocatable :: stdfsc !(ntyp)
+ real(kind=8),dimension(:),allocatable :: stdforcp,stdforcsc !(MAXRES)
+!-----------------------------------------------------------------------------
+! 'sizes.i'
+!
+!
+! ###################################################
+! ## COPYRIGHT (C) 1992 by Jay William Ponder ##
+! ## All Rights Reserved ##
+! ###################################################
+!
+! #############################################################
+! ## ##
+! ## sizes.i -- parameter values to set array dimensions ##
+! ## ##
+! #############################################################
+!
+!
+! "sizes.i" sets values for critical array dimensions used
+! throughout the software; these parameters will fix the size
+! of the largest systems that can be handled; values too large
+! for the computer's memory and/or swap space to accomodate
+! will result in poor performance or outright failure
+!
+! parameter: maximum allowed number of:
+!
+! maxatm atoms in the molecular system
+! maxval atoms directly bonded to an atom
+! maxgrp ! user-defined groups of atoms
+! maxtyp force field atom type definitions
+! maxclass force field atom class definitions
+! maxkey lines in the keyword file
+! maxrot bonds for torsional rotation
+! maxvar optimization variables (vector storage)
+! maxopt optimization variables (matrix storage)
+! maxhess off-diagonal Hessian elements
+! maxlight sites for method of lights neighbors
+! maxvib vibrational frequencies
+! maxgeo distance geometry points
+! maxcell unit cells in replicated crystal
+! maxring 3-, 4-, or 5-membered rings
+! maxfix geometric restraints
+! maxbio biopolymer atom definitions
+! maxres residues in the macromolecule
+! maxamino amino acid residue types
+! maxnuc nucleic acid residue types
+! maxbnd covalent bonds in molecular system
+! maxang bond angles in molecular system
+! maxtors torsional angles in molecular system
+! maxpi atoms in conjugated pisystem
+! maxpib covalent bonds involving pisystem
+! maxpit torsional angles involving pisystem
+!
+!
+!el integer maxatm,maxval,maxgrp
+!el integer maxtyp,maxclass,maxkey
+!el integer maxrot,maxopt
+!el integer maxhess,maxlight,maxvib
+!el integer maxgeo,maxcell,maxring
+!el integer maxfix,maxbio
+!el integer maxamino,maxnuc,maxbnd
+!el integer maxang,maxtors,maxpi
+!el integer maxpib,maxpit
+ integer :: maxatm !=2*nres !maxres2 maxres2=2*maxres
+ integer,parameter :: maxval=8
+ integer,parameter :: maxgrp=1000
+ integer,parameter :: maxtyp=3000
+ integer,parameter :: maxclass=500
+ integer,parameter :: maxkey=10000
+ integer,parameter :: maxrot=1000
+ integer,parameter :: maxopt=1000
+ integer,parameter :: maxhess=1000000
+ integer :: maxlight !=8*maxatm
+ integer,parameter :: maxvib=1000
+ integer,parameter :: maxgeo=1000
+ integer,parameter :: maxcell=10000
+ integer,parameter :: maxring=10000
+ integer,parameter :: maxfix=10000
+ integer,parameter :: maxbio=10000
+ integer,parameter :: maxamino=31
+ integer,parameter :: maxnuc=12
+ integer :: maxbnd !=2*maxatm
+ integer :: maxang !=3*maxatm
+ integer :: maxtors !=4*maxatm
+ integer,parameter :: maxpi=100
+ integer,parameter :: maxpib=2*maxpi
+ integer,parameter :: maxpit=4*maxpi
+!-----------------------------------------------------------------------------
+! Maximum number of seed
+ integer,parameter :: max_seed=1
+!-----------------------------------------------------------------------------
+ real(kind=8),dimension(:),allocatable :: stochforcvec !(MAXRES6) maxres6=6*maxres
+! common /stochcalc/ stochforcvec
+!-----------------------------------------------------------------------------
+! common /przechowalnia/ subroutines: rattle1,rattle2,rattle_brown
+ real(kind=8),dimension(:,:),allocatable :: GGinv !(2*nres,2*nres) maxres2=2*maxres
+ real(kind=8),dimension(:,:,:),allocatable :: gdc !(3,2*nres,2*nres) maxres2=2*maxres
+ real(kind=8),dimension(:,:),allocatable :: Cmat !(2*nres,2*nres) maxres2=2*maxres
+!-----------------------------------------------------------------------------
+! common /syfek/ subroutines: friction_force,setup_fricmat
+!el real(kind=8),dimension(:),allocatable :: gamvec !(MAXRES6) or (MAXRES2)
+!-----------------------------------------------------------------------------
+! common /przechowalnia/ subroutines: friction_force,setup_fricmat
+ real(kind=8),dimension(:,:),allocatable :: ginvfric !(2*nres,2*nres) !maxres2=2*maxres
+!-----------------------------------------------------------------------------
+! common /przechowalnia/ subroutine: setup_fricmat
+!el real(kind=8),dimension(:,:),allocatable :: fcopy !(2*nres,2*nres)
+!-----------------------------------------------------------------------------
+!
+!
+!-----------------------------------------------------------------------------
+ contains
+!-----------------------------------------------------------------------------
+! brown_step.f
+!-----------------------------------------------------------------------------
+ subroutine brown_step(itime)
+!------------------------------------------------
+! Perform a single Euler integration step of Brownian dynamics
+!------------------------------------------------
+! implicit real*8 (a-h,o-z)
+ use comm_gucio
+ use control, only: tcpu
+ use control_data
+ use energy_data
+! use io_conf, only:cartprint
+! 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'
+ real(kind=8),dimension(6*nres) :: zapas !(MAXRES6) maxres6=6*maxres
+ integer :: rstcount !ilen,
+!el external ilen
+!el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
+ real(kind=8),dimension(6*nres,2*nres) :: Bmat,GBmat,Tmat !(MAXRES6,MAXRES2) (maxres2=2*maxres,maxres6=6*maxres)
+ real(kind=8),dimension(2*nres,2*nres) :: Cmat_,Cinv !(maxres2,maxres2) maxres2=2*maxres
+ real(kind=8),dimension(6*nres,6*nres) :: Pmat !(maxres6,maxres6) maxres6=6*maxres
+ real(kind=8),dimension(6*nres) :: Td !(maxres6) maxres6=6*maxres
+ real(kind=8),dimension(2*nres) :: ppvec !(maxres2) maxres2=2*maxres
+!el common /stochcalc/ stochforcvec
+!el real(kind=8),dimension(3) :: cm !el
+!el common /gucio/ cm
+ integer :: itime
+ logical :: lprn = .false.,lprn1 = .false.
+ integer :: maxiter = 5
+ real(kind=8) :: difftol = 1.0d-5
+ real(kind=8) :: xx,diffmax,blen2,diffbond,tt0
+ integer :: i,j,nbond,k,ind,ind1,iter
+ integer :: nres2,nres6
+ logical :: osob
+ nres2=2*nres
+ nres6=6*nres
+
+ if (.not.allocated(stochforcvec)) allocate(stochforcvec(nres6)) !(MAXRES6) maxres6=6*maxres
+
+ nbond=nct-nnt
+ do i=nnt,nct
+ if (itype(i).ne.10) nbond=nbond+1
+ enddo
+!
+ if (lprn1) then
+ write (iout,*) "Generalized inverse of fricmat"
+ call matout(dimen,dimen,nres6,nres6,fricmat)
+ endif
+ do i=1,dimen
+ do j=1,nbond
+ Bmat(i,j)=0.0d0
+ enddo
+ enddo
+ ind=3
+ ind1=0
+ do i=nnt,nct-1
+ ind1=ind1+1
+ do j=1,3
+ Bmat(ind+j,ind1)=dC_norm(j,i)
+ enddo
+ ind=ind+3
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ ind1=ind1+1
+ do j=1,3
+ Bmat(ind+j,ind1)=dC_norm(j,i+nres)
+ enddo
+ ind=ind+3
+ endif
+ enddo
+ if (lprn1) then
+ write (iout,*) "Matrix Bmat"
+ call MATOUT(nbond,dimen,nres6,nres6,Bmat)
+ endif
+ do i=1,dimen
+ do j=1,nbond
+ GBmat(i,j)=0.0d0
+ do k=1,dimen
+ GBmat(i,j)=GBmat(i,j)+fricmat(i,k)*Bmat(k,j)
+ enddo
+ enddo
+ enddo
+ if (lprn1) then
+ write (iout,*) "Matrix GBmat"
+ call MATOUT(nbond,dimen,nres6,nres2,Gbmat)
+ endif
+ do i=1,nbond
+ do j=1,nbond
+ Cmat_(i,j)=0.0d0
+ do k=1,dimen
+ Cmat_(i,j)=Cmat_(i,j)+Bmat(k,i)*GBmat(k,j)
+ enddo
+ enddo
+ enddo
+ if (lprn1) then
+ write (iout,*) "Matrix Cmat"
+ call MATOUT(nbond,nbond,nres2,nres2,Cmat_)
+ endif
+ call matinvert(nbond,nres2,Cmat_,Cinv,osob)
+ if (lprn1) then
+ write (iout,*) "Matrix Cinv"
+ call MATOUT(nbond,nbond,nres2,nres2,Cinv)
+ endif
+ do i=1,dimen
+ do j=1,nbond
+ Tmat(i,j)=0.0d0
+ do k=1,nbond
+ Tmat(i,j)=Tmat(i,j)+GBmat(i,k)*Cinv(k,j)
+ enddo
+ enddo
+ enddo
+ if (lprn1) then
+ write (iout,*) "Matrix Tmat"
+ call MATOUT(nbond,dimen,nres6,nres2,Tmat)
+ endif
+ do i=1,dimen
+ do j=1,dimen
+ if (i.eq.j) then
+ Pmat(i,j)=1.0d0
+ else
+ Pmat(i,j)=0.0d0
+ endif
+ do k=1,nbond
+ Pmat(i,j)=Pmat(i,j)-Tmat(i,k)*Bmat(j,k)
+ enddo
+ enddo
+ enddo
+ if (lprn1) then
+ write (iout,*) "Matrix Pmat"
+ call MATOUT(dimen,dimen,nres6,nres6,Pmat)
+ endif
+ do i=1,dimen
+ Td(i)=0.0d0
+ ind=0
+ do k=nnt,nct-1
+ ind=ind+1
+ Td(i)=Td(i)+vbl*Tmat(i,ind)
+ enddo
+ do k=nnt,nct
+ if (itype(k).ne.10) then
+ ind=ind+1
+ Td(i)=Td(i)+vbldsc0(1,itype(k))*Tmat(i,ind)
+ endif
+ enddo
+ enddo
+ if (lprn1) then
+ write (iout,*) "Vector Td"
+ do i=1,dimen
+ write (iout,'(i5,f10.5)') i,Td(i)
+ enddo
+ endif
+ call stochastic_force(stochforcvec)
+ if (lprn) then
+ write (iout,*) "stochforcvec"
+ do i=1,dimen
+ write (iout,*) i,stochforcvec(i)
+ enddo
+ endif
+ do j=1,3
+ zapas(j)=-gcart(j,0)+stochforcvec(j)
+ d_t_work(j)=d_t(j,0)
+ dC_work(j)=dC_old(j,0)
+ enddo
+ ind=3
+ do i=nnt,nct-1
+ do j=1,3
+ ind=ind+1
+ zapas(ind)=-gcart(j,i)+stochforcvec(ind)
+ dC_work(ind)=dC_old(j,i)
+ enddo
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ do j=1,3
+ ind=ind+1
+ zapas(ind)=-gxcart(j,i)+stochforcvec(ind)
+ dC_work(ind)=dC_old(j,i+nres)
+ enddo
+ endif
+ enddo
+
+ if (lprn) then
+ write (iout,*) "Initial d_t_work"
+ do i=1,dimen
+ write (iout,*) i,d_t_work(i)
+ enddo
+ endif
+
+ do i=1,dimen
+ d_t_work(i)=0.0d0
+ do j=1,dimen
+ d_t_work(i)=d_t_work(i)+fricmat(i,j)*zapas(j)
+ enddo
+ enddo
+
+ do i=1,dimen
+ zapas(i)=Td(i)
+ do j=1,dimen
+ zapas(i)=zapas(i)+Pmat(i,j)*(dC_work(j)+d_t_work(j)*d_time)
+ enddo
+ enddo
+ if (lprn1) then
+ write (iout,*) "Final d_t_work and zapas"
+ do i=1,dimen
+ write (iout,*) i,d_t_work(i),zapas(i)
+ enddo
+ endif
+
+ do j=1,3
+ d_t(j,0)=d_t_work(j)
+ dc(j,0)=zapas(j)
+ dc_work(j)=dc(j,0)
+ enddo
+ ind=3
+ do i=nnt,nct-1
+ do j=1,3
+ d_t(j,i)=d_t_work(i)
+ dc(j,i)=zapas(ind+j)
+ dc_work(ind+j)=dc(j,i)
+ enddo
+ ind=ind+3
+ enddo
+ do i=nnt,nct
+ do j=1,3
+ d_t(j,i+nres)=d_t_work(ind+j)
+ dc(j,i+nres)=zapas(ind+j)
+ dc_work(ind+j)=dc(j,i+nres)
+ enddo
+ ind=ind+3
+ enddo
+ if (lprn) then
+ call chainbuild_cart
+ write (iout,*) "Before correction for rotational lengthening"
+ write (iout,*) "New coordinates",&
+ " and differences between actual and standard bond lengths"
+ ind=0
+ do i=nnt,nct-1
+ ind=ind+1
+ xx=vbld(i+1)-vbl
+ write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
+ i,(dC(j,i),j=1,3),xx
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ ind=ind+1
+ xx=vbld(i+nres)-vbldsc0(1,itype(i))
+ write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
+ i,(dC(j,i+nres),j=1,3),xx
+ endif
+ enddo
+ endif
+! Second correction (rotational lengthening)
+! do iter=1,maxiter
+ diffmax=0.0d0
+ ind=0
+ do i=nnt,nct-1
+ ind=ind+1
+ blen2 = scalar(dc(1,i),dc(1,i))
+ ppvec(ind)=2*vbl**2-blen2
+ diffbond=dabs(vbl-dsqrt(blen2))
+ if (diffbond.gt.diffmax) diffmax=diffbond
+ if (ppvec(ind).gt.0.0d0) then
+ ppvec(ind)=dsqrt(ppvec(ind))
+ else
+ ppvec(ind)=0.0d0
+ endif
+ if (lprn) then
+ write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind)
+ endif
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ ind=ind+1
+ blen2 = scalar(dc(1,i+nres),dc(1,i+nres))
+ ppvec(ind)=2*vbldsc0(1,itype(i))**2-blen2
+ diffbond=dabs(vbldsc0(1,itype(i))-dsqrt(blen2))
+ if (diffbond.gt.diffmax) diffmax=diffbond
+ if (ppvec(ind).gt.0.0d0) then
+ ppvec(ind)=dsqrt(ppvec(ind))
+ else
+ ppvec(ind)=0.0d0
+ endif
+ if (lprn) then
+ write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind)
+ endif
+ endif
+ enddo
+ if (lprn) write (iout,*) "iter",iter," diffmax",diffmax
+ if (diffmax.lt.difftol) goto 10
+ do i=1,dimen
+ Td(i)=0.0d0
+ do j=1,nbond
+ Td(i)=Td(i)+ppvec(j)*Tmat(i,j)
+ enddo
+ enddo
+ do i=1,dimen
+ zapas(i)=Td(i)
+ do j=1,dimen
+ zapas(i)=zapas(i)+Pmat(i,j)*dc_work(j)
+ enddo
+ enddo
+ do j=1,3
+ dc(j,0)=zapas(j)
+ dc_work(j)=zapas(j)
+ enddo
+ ind=3
+ do i=nnt,nct-1
+ do j=1,3
+ dc(j,i)=zapas(ind+j)
+ dc_work(ind+j)=zapas(ind+j)
+ enddo
+ ind=ind+3
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ do j=1,3
+ dc(j,i+nres)=zapas(ind+j)
+ dc_work(ind+j)=zapas(ind+j)
+ enddo
+ ind=ind+3
+ endif
+ enddo
+! Building the chain from the newly calculated coordinates
+ call chainbuild_cart
+ 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,'(a)') "Potential forces"
+ do i=0,nres
+ write (iout,'(i3,3f10.5,3x,3f10.5)') i,(-gcart(j,i),j=1,3),&
+ (-gxcart(j,i),j=1,3)
+ enddo
+ write (iout,'(a)') "Stochastic forces"
+ do i=0,nres
+ write (iout,'(i3,3f10.5,3x,3f10.5)') i,(stochforc(j,i),j=1,3),&
+ (stochforc(j,i+nres),j=1,3)
+ enddo
+ write (iout,'(a)') "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
+ endif
+ endif
+ if (lprn) then
+ write (iout,*) "After correction for rotational lengthening"
+ write (iout,*) "New coordinates",&
+ " and differences between actual and standard bond lengths"
+ ind=0
+ do i=nnt,nct-1
+ ind=ind+1
+ xx=vbld(i+1)-vbl
+ write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
+ i,(dC(j,i),j=1,3),xx
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ ind=ind+1
+ xx=vbld(i+nres)-vbldsc0(1,itype(i))
+ write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
+ i,(dC(j,i+nres),j=1,3),xx
+ endif
+ enddo
+ endif
+! ENDDO
+! write (iout,*) "Too many attempts at correcting the bonds"
+! stop
+ 10 continue
+#ifdef MPI
+ tt0 =MPI_Wtime()
+#else
+ tt0 = tcpu()
+#endif
+! Calculate energy and forces
+ call zerograd
+ call etotal(potEcomp)
+ potE=potEcomp(0)-potEcomp(20)
+ call cartgrad
+ totT=totT+d_time
+! Calculate the kinetic and total energy and the kinetic temperature
+ call kinetic(EK)
+#ifdef MPI
+ t_enegrad=t_enegrad+MPI_Wtime()-tt0
+#else
+ t_enegrad=t_enegrad+tcpu()-tt0
+#endif
+ totE=EK+potE
+ kinetic_T=2.0d0/(dimen*Rb)*EK
+ return
+ end subroutine brown_step
+!-----------------------------------------------------------------------------
+! gauss.f
+!-----------------------------------------------------------------------------
+ subroutine gauss(RO,AP,MT,M,N,*)
+!
+! CALCULATES (RO**(-1))*AP BY GAUSS ELIMINATION
+! RO IS A SQUARE MATRIX
+! THE CALCULATED PRODUCT IS STORED IN AP
+! ABNORMAL EXIT IF RO IS SINGULAR
+!
+ integer :: MT, M, N, M1,I,J,IM,&
+ I1,MI,MI1
+ real(kind=8) :: RO(MT,M),AP(MT,N),X,RM,PR,Y
+ integer :: k
+! real(kind=8) ::
+
+ if(M.ne.1)goto 10
+ X=RO(1,1)
+ if(dabs(X).le.1.0D-13) return 1
+ X=1.0/X
+ do 16 I=1,N
+16 AP(1,I)=AP(1,I)*X
+ return
+10 continue
+ M1=M-1
+ DO 1 I=1,M1
+ IM=I
+ RM=DABS(RO(I,I))
+ I1=I+1
+ do 2 J=I1,M
+ if(DABS(RO(J,I)).LE.RM) goto 2
+ RM=DABS(RO(J,I))
+ IM=J
+2 continue
+ If(IM.eq.I)goto 17
+ do 3 J=1,N
+ PR=AP(I,J)
+ AP(I,J)=AP(IM,J)
+3 AP(IM,J)=PR
+ do 4 J=I,M
+ PR=RO(I,J)
+ RO(I,J)=RO(IM,J)
+4 RO(IM,J)=PR
+17 X=RO(I,I)
+ if(dabs(X).le.1.0E-13) return 1
+ X=1.0/X
+ do 5 J=1,N
+5 AP(I,J)=X*AP(I,J)
+ do 6 J=I1,M
+6 RO(I,J)=X*RO(I,J)
+ do 7 J=I1,M
+ Y=RO(J,I)
+ do 8 K=1,N
+8 AP(J,K)=AP(J,K)-Y*AP(I,K)
+ do 9 K=I1,M
+9 RO(J,K)=RO(J,K)-Y*RO(I,K)
+7 continue
+1 continue
+ X=RO(M,M)
+ if(dabs(X).le.1.0E-13) return 1
+ X=1.0/X
+ do 11 J=1,N
+11 AP(M,J)=X*AP(M,J)
+ do 12 I=1,M1
+ MI=M-I
+ MI1=MI+1
+ do 14 J=1,N
+ X=AP(MI,J)
+ do 15 K=MI1,M
+15 X=X-AP(K,J)*RO(MI,K)
+14 AP(MI,J)=X
+12 continue
+ return
+ end subroutine gauss
+!-----------------------------------------------------------------------------
+! kinetic_lesyng.f
+!-----------------------------------------------------------------------------
+ subroutine kinetic(KE_total)
+!----------------------------------------------------------------
+! This subroutine calculates the total kinetic energy of the chain
+!-----------------------------------------------------------------
+ use energy_data
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.VAR'
+! include 'COMMON.CHAIN'
+! include 'COMMON.DERIV'
+! include 'COMMON.GEO'
+! include 'COMMON.LOCAL'
+! include 'COMMON.INTERACT'
+! include 'COMMON.MD'
+! include 'COMMON.IOUNITS'
+ real(kind=8) :: KE_total
+
+ integer :: i,j,k,iti
+ real(kind=8) :: KEt_p,KEt_sc,KEr_p,KEr_sc,incr(3),&
+ mag1,mag2,v(3)
+
+ KEt_p=0.0d0
+ KEt_sc=0.0d0
+! write (iout,*) "ISC",(isc(itype(i)),i=1,nres)
+! The translational part for peptide virtual bonds
+ do j=1,3
+ incr(j)=d_t(j,0)
+ enddo
+ do i=nnt,nct-1
+! write (iout,*) "Kinetic trp:",i,(incr(j),j=1,3)
+ do j=1,3
+ v(j)=incr(j)+0.5d0*d_t(j,i)
+ enddo
+ vtot(i)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
+ KEt_p=KEt_p+(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
+ do j=1,3
+ incr(j)=incr(j)+d_t(j,i)
+ enddo
+ enddo
+! write(iout,*) 'KEt_p', KEt_p
+! The translational part for the side chain virtual bond
+! Only now we can initialize incr with zeros. It must be equal
+! to the velocities of the first Calpha.
+ do j=1,3
+ incr(j)=d_t(j,0)
+ enddo
+ do i=nnt,nct
+ iti=iabs(itype(i))
+ if (itype(i).eq.10) then
+ do j=1,3
+ v(j)=incr(j)
+ enddo
+ else
+ do j=1,3
+ v(j)=incr(j)+d_t(j,nres+i)
+ enddo
+ endif
+! write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3)
+! write (iout,*) "i",i," msc",msc(iti)," v",(v(j),j=1,3)
+ KEt_sc=KEt_sc+msc(iti)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
+ vtot(i+nres)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
+ do j=1,3
+ incr(j)=incr(j)+d_t(j,i)
+ enddo
+ enddo
+! goto 111
+! write(iout,*) 'KEt_sc', KEt_sc
+! The part due to stretching and rotation of the peptide groups
+ KEr_p=0.0D0
+ do i=nnt,nct-1
+! write (iout,*) "i",i
+! write (iout,*) "i",i," mag1",mag1," mag2",mag2
+ do j=1,3
+ incr(j)=d_t(j,i)
+ enddo
+! write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3)
+ KEr_p=KEr_p+(incr(1)*incr(1)+incr(2)*incr(2) &
+ +incr(3)*incr(3))
+ enddo
+! goto 111
+! write(iout,*) 'KEr_p', KEr_p
+! The rotational part of the side chain virtual bond
+ KEr_sc=0.0D0
+ do i=nnt,nct
+ iti=iabs(itype(i))
+ if (itype(i).ne.10) then
+ do j=1,3
+ incr(j)=d_t(j,nres+i)
+ enddo
+! write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3)
+ KEr_sc=KEr_sc+Isc(iti)*(incr(1)*incr(1)+incr(2)*incr(2)+ &
+ incr(3)*incr(3))
+ endif
+ enddo
+! The total kinetic energy
+ 111 continue
+! write(iout,*) 'KEr_sc', KEr_sc
+ KE_total=0.5d0*(mp*KEt_p+KEt_sc+0.25d0*Ip*KEr_p+KEr_sc)
+! write (iout,*) "KE_total",KE_total
+ return
+ end subroutine kinetic
+!-----------------------------------------------------------------------------
+! MD_A-MTS.F
+!-----------------------------------------------------------------------------
+ subroutine MD
+!------------------------------------------------
+! The driver for molecular dynamics subroutines
+!------------------------------------------------
+ use comm_gucio
+! use MPI
+ use control, only:tcpu,ovrtim
+! use io_comm, only:ilen
+ use control_data
+ use compare, only:secondary2,hairpin
+ use io, only:cartout,statout
+! 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'
+! include 'COMMON.HAIRPIN'
+ real(kind=8),dimension(3) :: L,vcm
+#ifdef VOUT
+ real(kind=8),dimension(6*nres) :: v_work,v_transf !(maxres6) maxres6=6*maxres
+#endif
+ integer :: rstcount !ilen,
+!el external ilen
+ character(len=50) :: tytul
+!el common /gucio/ cm
+ integer :: itime,i,j,nharp
+ integer,dimension(4,nres/3) :: iharp !(4,nres/3)(4,maxres/3)
+! logical :: ovrtim
+ real(kind=8) :: tt0,scalfac
+ integer :: nres2
+ nres2=2*nres
+!
+#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
+! Determine the inverse of the inertia matrix.
+ call setup_MD_matrices
+! Initialize MD
+ call init_MD
+#ifdef MPI
+ t_MDsetup = MPI_Wtime()-tt0
+#else
+ t_MDsetup = tcpu()-tt0
+#endif
+ rstcount=0
+! 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
+ if (ovrtim()) exit
+ if (large.and. mod(itime,ntwe).eq.0) &
+ write (iout,*) "itime",itime
+ 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
+! Time-reversible RESPA algorithm
+! (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
+ call RESPA_step(itime)
+ else
+! 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 .and. itype(i).ne.ntyp1) 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 hairpin(.true.,nharp,iharp)
+ call secondary2(.true.)
+ 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 subroutine MD
+!-----------------------------------------------------------------------------
+ subroutine velverlet_step(itime)
+!-------------------------------------------------------------------------------
+! Perform a single velocity Verlet step; the time step can be rescaled if
+! increments in accelerations exceed the threshold
+!-------------------------------------------------------------------------------
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+ use comm_gucio
+ use control, only:tcpu
+ use control_data
+#ifdef MPI
+ include 'mpif.h'
+ integer :: ierror,ierrcode
+ real(kind=8) :: errcode
+#endif
+! include 'COMMON.SETUP'
+! 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'
+ real(kind=8),dimension(3) :: vcm,incr
+ real(kind=8),dimension(3) :: L
+ integer :: count,rstcount !ilen,
+!el external ilen
+ character(len=50) :: tytul
+ integer :: maxcount_scale = 20
+!el common /gucio/ cm
+!el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
+!el common /stochcalc/ stochforcvec
+ integer :: itime,icount_scale,itime_scal,i,j,ifac_time
+ logical :: scale
+ real(kind=8) :: epdrift,tt0,fac_time
+!
+ if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres)) !(MAXRES6) maxres6=6*maxres
+
+ 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
+! 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
+! 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
+! Calculate energy and forces
+ call zerograd
+ call etotal(potEcomp)
+ if (large.and. mod(itime,ntwe).eq.0) &
+ call enerprint(potEcomp)
+#ifdef TIMING_ENE
+#ifdef MPI
+ t_etotal=t_etotal+MPI_Wtime()-tt0
+#else
+ t_etotal=t_etotal+tcpu()-tt0
+#endif
+#endif
+ potE=potEcomp(0)-potEcomp(20)
+ call cartgrad
+! Get the new accelerations
+ call lagrangian
+#ifdef MPI
+ t_enegrad=t_enegrad+MPI_Wtime()-tt0
+#else
+ t_enegrad=t_enegrad+tcpu()-tt0
+#endif
+! 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
+! 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
+! 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
+! write (iout,*) "Calling sd_verlet_setup: 1"
+! Rescale the stochastic forces and recalculate or restore
+! 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
+! 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
+! 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"
+! write (iout,*) "Calling sd_verlet_setup: 2"
+! 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
+! Calculate the kinetic and the total energy and the kinetic temperature
+ call kinetic(EK)
+ totE=EK+potE
+! diagnostics
+! call kinetic1(EK1)
+! write (iout,*) "step",itime," EK",EK," EK1",EK1
+! end diagnostics
+! 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
+! 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 subroutine velverlet_step
+!-----------------------------------------------------------------------------
+ subroutine RESPA_step(itime)
+!-------------------------------------------------------------------------------
+! Perform a single RESPA step.
+!-------------------------------------------------------------------------------
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+ use comm_gucio
+ use comm_cipiszcze
+! use MPI
+ use control, only:tcpu
+ use control_data
+! use io_conf, only:cartprint
+#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'
+ real(kind=8),dimension(0:n_ene) :: energia_short,energia_long
+ real(kind=8),dimension(3) :: L,vcm,incr
+ real(kind=8),dimension(3,0:2*nres) :: dc_old0,d_t_old0,d_a_old0 !(3,0:maxres2) maxres2=2*maxres
+ logical :: PRINT_AMTS_MSG = .false.
+ integer :: count,rstcount !ilen,
+!el external ilen
+ character(len=50) :: tytul
+ integer :: maxcount_scale = 10
+!el common /gucio/ cm
+!el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
+!el common /stochcalc/ stochforcvec
+ integer :: itime,itt,i,j,itsplit
+ logical :: scale
+!el common /cipiszcze/ itt
+
+ real(kind=8) :: epdrift,tt0,epdriftmax
+ itt = itt_comm
+
+ if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres)) !(MAXRES6) maxres6=6*maxres
+
+ 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"
+! 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
+!
+! Perform the initial RESPA step (increment velocities)
+! 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
+! Compute the short-range forces
+#ifdef MPI
+ tt0 =MPI_Wtime()
+#else
+ tt0 = tcpu()
+#endif
+! 7/2/2009 commented out
+! call zerograd
+! call etotal_short(energia_short)
+! call cartgrad
+! call lagrangian
+! 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
+ do i=0,2*nres
+ do j=1,3
+ d_a(j,i)=d_a_short(j,i)
+ enddo
+ enddo
+ 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
+! 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
+!
+ scale=.true.
+ d_time0=d_time
+ do while (scale)
+
+ scale=.false.
+! write (iout,*) "itime",itime," ntime_split",ntime_split
+! Split the time step
+ d_time=d_time0/ntime_split
+! Perform the short-range RESPA steps (velocity Verlet increments of
+! positions and velocities using short-range forces)
+! 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
+! 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
+! 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)
+! 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
+! Calculate energy and forces
+ call zerograd
+ call etotal_short(energia_short)
+ if (large.and. mod(itime,ntwe).eq.0) &
+ call enerprint(energia_short)
+#ifdef TIMING_ENE
+#ifdef MPI
+ t_eshort=t_eshort+MPI_Wtime()-tt0
+#else
+ t_eshort=t_eshort+tcpu()-tt0
+#endif
+#endif
+ call cartgrad
+! Get the new accelerations
+ call lagrangian
+! 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
+ do i=0,2*nres
+ do j=1,3
+ d_a_short(j,i)=d_a(j,i)
+ enddo
+ enddo
+ 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
+! 6/30/08 A-MTS
+! 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
+! Exit loop and try with increased split number if the change of
+! 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
+ if (PRINT_AMTS_MSG) then
+ write (iout,*) "acceleration/energy drift too large",amax,&
+ epdrift," split increased to ",ntime_split," itime",itime,&
+ " itsplit",itsplit
+ endif
+ 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
+! 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
+! 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
+
+! Restore the time step
+ d_time=d_time0
+! Compute long-range forces
+#ifdef MPI
+ tt0 =MPI_Wtime()
+#else
+ tt0 = tcpu()
+#endif
+ call zerograd
+ call etotal_long(energia_long)
+ if (large.and. mod(itime,ntwe).eq.0) &
+ call enerprint(energia_long)
+#ifdef TIMING_ENE
+#ifdef MPI
+ t_elong=t_elong+MPI_Wtime()-tt0
+#else
+ t_elong=t_elong+tcpu()-tt0
+#endif
+#endif
+ call cartgrad
+ call lagrangian
+#ifdef MPI
+ t_enegrad=t_enegrad+MPI_Wtime()-tt0
+#else
+ t_enegrad=t_enegrad+tcpu()-tt0
+#endif
+! 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"
+! 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
+! Compute the final RESPA step (increment velocities)
+! write (iout,*) "*********************** RESPA fin"
+ call RESPA_vel
+! Compute the complete potential energy
+ do i=0,n_ene
+ potEcomp(i)=energia_short(i)+energia_long(i)
+ enddo
+ potE=potEcomp(0)-potEcomp(20)
+! potE=energia_short(0)+energia_long(0)
+ totT=totT+d_time
+! Calculate the kinetic and the total energy and the kinetic temperature
+ call kinetic(EK)
+ totE=EK+potE
+! 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
+! 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 subroutine RESPA_step
+!-----------------------------------------------------------------------------
+ subroutine RESPA_vel
+! First and last RESPA step (incrementing velocities using long-range
+! forces).
+ use energy_data
+! 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'
+ 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
+ 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 .and. itype(i).ne.ntyp1) 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 subroutine RESPA_vel
+!-----------------------------------------------------------------------------
+ subroutine verlet1
+! Applying velocity Verlet algorithm - step 1 to coordinates
+ use energy_data
+! 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'
+ real(kind=8) :: adt,adt2
+ integer :: i,j,inres
+
+#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 .and. itype(i).ne.ntyp1) 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 subroutine verlet1
+!-----------------------------------------------------------------------------
+ subroutine verlet2
+! Step 2 of the velocity Verlet algorithm: update velocities
+ use energy_data
+! 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'
+ 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
+ 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 .and. itype(i).ne.ntyp1) 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 subroutine verlet2
+!-----------------------------------------------------------------------------
+ subroutine sddir_precalc
+! Applying velocity Verlet algorithm - step 1 to coordinates
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+ use MPI_data
+ use control_data
+#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'
+!el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
+!el common /stochcalc/ stochforcvec
+ real(kind=8) :: time00
+!
+! Compute friction and stochastic forces
+!
+#ifdef MPI
+ time00=MPI_Wtime()
+ call friction_force
+ time_fric=time_fric+MPI_Wtime()-time00
+ time00=MPI_Wtime()
+ call stochastic_force(stochforcvec)
+ time_stoch=time_stoch+MPI_Wtime()-time00
+#endif
+!
+! Compute the acceleration due to friction forces (d_af_work) and stochastic
+! forces (d_as_work)
+!
+ call ginv_mult(fric_work, d_af_work)
+ call ginv_mult(stochforcvec, d_as_work)
+ return
+ end subroutine sddir_precalc
+!-----------------------------------------------------------------------------
+ subroutine sddir_verlet1
+! Applying velocity Verlet algorithm - step 1 to velocities
+!
+ use energy_data
+! 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'
+! Revised 3/31/05 AL: correlation between random contributions to
+! position and velocity increments included.
+ real(kind=8) :: sqrt13 = 0.57735026918962576451d0 ! 1/sqrt(3)
+ real(kind=8) :: adt,adt2
+ integer :: i,j,ind,inres
+!
+! Add the contribution from BOTH friction and stochastic force to the
+! coordinates, but ONLY the contribution from the friction forces to velocities
+!
+ 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 .and. itype(i).ne.ntyp1) 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 subroutine sddir_verlet1
+!-----------------------------------------------------------------------------
+ subroutine sddir_verlet2
+! Calculating the adjusted velocities for accelerations
+!
+ use energy_data
+! 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'
+ real(kind=8),dimension(6*nres) :: stochforcvec,d_as_work1 !(MAXRES6) maxres6=6*maxres
+ real(kind=8) :: cos60 = 0.5d0, sin60 = 0.86602540378443864676d0
+ integer :: i,j,ind,inres
+! Revised 3/31/05 AL: correlation between random contributions to
+! position and velocity increments included.
+! The correlation coefficients are calculated at low-friction limit.
+! Also, friction forces are now not calculated with new velocities.
+
+! call friction_force
+ call stochastic_force(stochforcvec)
+!
+! Compute the acceleration due to friction forces (d_af_work) and stochastic
+! forces (d_as_work)
+!
+ call ginv_mult(stochforcvec, d_as_work1)
+
+!
+! Update velocities
+!
+ 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 .and. itype(i).ne.ntyp1) 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 subroutine sddir_verlet2
+!-----------------------------------------------------------------------------
+ subroutine max_accel
+!
+! Find the maximum difference in the accelerations of the the sites
+! at the beginning and the end of the time step.
+!
+ use energy_data
+! 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'
+ real(kind=8),dimension(3) :: aux,accel,accel_old
+ real(kind=8) :: dacc
+ integer :: i,j
+
+ do j=1,3
+! 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
+! Backbone
+ if (i.lt.nct) then
+! 7/3/08 changed to asymmetric difference
+ do j=1,3
+! 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)
+! 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))
+! write (iout,*) i,dacc
+ if (dacc.gt.amax) amax=dacc
+ endif
+ enddo
+ endif
+ enddo
+! Side chains
+ do j=1,3
+! 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 .and. itype(i).ne.ntyp1) then
+ do j=1,3
+! 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
+! 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))
+! write (iout,*) "side-chain",i,dacc
+ 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)
+! aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
+ enddo
+ enddo
+ return
+ end subroutine max_accel
+!-----------------------------------------------------------------------------
+ subroutine predict_edrift(epdrift)
+!
+! Predict the drift of the potential energy
+!
+ use energy_data
+ use control_data, only: lmuca
+! 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'
+ real(kind=8) :: epdrift,epdriftij
+ integer :: i,j
+! Drift of the potential energy
+ epdrift=0.0d0
+ do i=nnt,nct
+! 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
+! write (iout,*) "back",i,j,epdriftij
+ if (epdriftij.gt.epdrift) epdrift=epdriftij
+ enddo
+ endif
+! Side chains
+ if (itype(i).ne.10 .and. itype(i).ne.ntyp1) 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
+! write (iout,*) "side",i,j,epdriftij
+ if (epdriftij.gt.epdrift) epdrift=epdriftij
+ enddo
+ endif
+ enddo
+ epdrift=0.5d0*epdrift*d_time*d_time
+! write (iout,*) "epdrift",epdrift
+ return
+ end subroutine predict_edrift
+!-----------------------------------------------------------------------------
+ subroutine verlet_bath
+!
+! Coupling to the thermostat by using the Berendsen algorithm
+!
+ use energy_data
+! 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'
+ real(kind=8) :: T_half,fact
+ 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))
+! write(iout,*) "T_half", T_half
+! write(iout,*) "EK", EK
+! 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 .and. itype(i).ne.ntyp1) then
+ inres=i+nres
+ do j=1,3
+ d_t(j,inres)=fact*d_t(j,inres)
+ enddo
+ endif
+ enddo
+ return
+ end subroutine verlet_bath
+!-----------------------------------------------------------------------------
+ subroutine init_MD
+! Set up the initial conditions of a MD simulation
+ use comm_gucio
+ use energy_data
+ use control, only:tcpu
+!el use io_basic, only:ilen
+ use control_data
+ use MPI_data
+ use minimm, only:minim_dc,minimize,sc_move
+ use io_config, only:readrst
+ use io, only:statout
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+#ifdef MP
+ include 'mpif.h'
+ character(len=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(kind=8),dimension(0:n_ene) :: energia_long,energia_short
+ real(kind=8),dimension(3) :: vcm,incr,L
+ real(kind=8) :: xv,sigv,lowb,highb
+ real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
+ character(len=256) :: qstr
+!el integer ilen
+!el external ilen
+ character(len=50) :: tytul
+ logical :: file_exist
+!el common /gucio/ cm
+ integer :: i,j,ipos,iq,iw,nft_sc,iretcode,nfun,itime,ierr
+ real(kind=8) :: etot,tt0
+ logical :: fail
+
+ d_time0=d_time
+! write(iout,*) "d_time", d_time
+! Compute the standard deviations of stochastic forces for Langevin dynamics
+! 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(iabs(itype(i))) &
+ *dsqrt(gamsc(iabs(itype(i))))
+ enddo
+ endif
+! 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
+! pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
+#ifdef NOXDR
+! cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
+#else
+! cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
+#endif
+! 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
+#ifdef MPI !el
+ call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,&
+ IERR)
+#endif !el
+ write (*,*) me," After broadcast: file_exist",file_exist
+! 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
+! Generate initial velocities
+ if(me.eq.king.or..not.out1file) &
+ write(iout,*) "Initial velocities randomly generated"
+ call random_vel
+ totT=0.0d0
+ endif
+! rest2name = prefix(:ilen(prefix))//'.rst'
+ if(me.eq.king.or..not.out1file)then
+ 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
+! Zeroing the total angular momentum of the system
+ write(iout,*) "Calling the zero-angular momentum subroutine"
+ endif
+ call inertia_tensor
+! Getting the potential energy and forces and velocities and accelerations
+ call vcm_vel(vcm)
+! write (iout,*) "velocity of the center of the mass:"
+! write (iout,*) (vcm(j),j=1,3)
+ do j=1,3
+ d_t(j,0)=d_t(j,0)-vcm(j)
+ enddo
+! 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
+ endif
+ kinetic_T=2.0d0/(dimen3*Rb)*EK
+ if(me.eq.king.or..not.out1file)then
+ call cartprint
+ call intout
+ endif
+#ifdef MPI
+ tt0=MPI_Wtime()
+#else
+ tt0=tcpu()
+#endif
+ call zerograd
+ call etotal(potEcomp)
+ if (large) call enerprint(potEcomp)
+#ifdef TIMING_ENE
+#ifdef MPI
+ t_etotal=t_etotal+MPI_Wtime()-tt0
+#else
+ t_etotal=t_etotal+tcpu()-tt0
+#endif
+#endif
+ 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"
+ call enerprint(potEcomp)
+! 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,'(/a/3(a25,1pe14.5/))') "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"
+ write (iout,"(13x,' backbone ',23x,' side chain')")
+ 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
+! 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
+! 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_short(energia_short)
+ if (large) call enerprint(potEcomp)
+#ifdef TIMING_ENE
+#ifdef MPI
+ t_eshort=t_eshort+MPI_Wtime()-tt0
+#else
+ t_eshort=t_eshort+tcpu()-tt0
+#endif
+#endif
+ call cartgrad
+ call lagrangian
+ if(.not.out1file .and. large) then
+ write (iout,*) "energia_long",energia_long(0),&
+ " energia_short",energia_short(0),&
+ " total",energia_long(0)+energia_short(0)
+ write (iout,*) "Initial fast-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
+! 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
+ do i=0,2*nres
+ do j=1,3
+ d_a_short(j,i)=d_a(j,i)
+ enddo
+ enddo
+#ifdef MPI
+ tt0=MPI_Wtime()
+#else
+ tt0=tcpu()
+#endif
+ call zerograd
+ call etotal_long(energia_long)
+ if (large) call enerprint(potEcomp)
+#ifdef TIMING_ENE
+#ifdef MPI
+ t_elong=t_elong+MPI_Wtime()-tt0
+#else
+ t_elong=t_elong+tcpu()-tt0
+#endif
+#endif
+ call cartgrad
+ call lagrangian
+ 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
+#ifdef MPI
+ t_enegrad=t_enegrad+MPI_Wtime()-tt0
+#else
+ t_enegrad=t_enegrad+tcpu()-tt0
+#endif
+ endif
+ return
+ end subroutine init_MD
+!-----------------------------------------------------------------------------
+ subroutine random_vel
+
+! implicit real*8 (a-h,o-z)
+ use energy_data
+ use random, only:anorm_distr
+! 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'
+ real(kind=8) :: xv,sigv,lowb,highb ,Ek1
+ integer :: i,j,ii,k,ind
+! Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
+! First generate velocities in the eigenspace of the G matrix
+! write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
+! 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)
+! write (iout,*) "i",i," ii",ii," geigen",geigen(i),&
+! " d_t_work_new",d_t_work_new(ii)
+ enddo
+ enddo
+! diagnostics
+! Ek1=0.0d0
+! ii=0
+! do i=1,dimen
+! do k=1,3
+! ii=ii+1
+! Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
+! enddo
+! enddo
+! write (iout,*) "Ek from eigenvectors",Ek1
+! end diagnostics
+! 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
+! write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
+! call flush(iout)
+ enddo
+ enddo
+! 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 .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
+! call kinetic(EK)
+! write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",&
+! 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
+! call flush(iout)
+ return
+ end subroutine random_vel
+!-----------------------------------------------------------------------------
+#ifndef LANG0
+ subroutine sd_verlet_p_setup
+! Sets up the parameters of stochastic Verlet algorithm
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+ use control, only: tcpu
+ use control_data
+#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'
+ real(kind=8),dimension(6*nres) :: emgdt !(MAXRES6) maxres6=6*maxres
+ real(kind=8) :: pterm,vterm,rho,rhoc,vsig
+ real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
+ prand_vec,vrand_vec1,vrand_vec2 !(MAXRES6) maxres6=6*maxres
+ logical :: lprn = .false.
+ real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
+ real(kind=8) :: ktm,gdt,egdt,gdt2,gdt3,gdt4,gdt5,gdt6,gdt7,gdt8,&
+ gdt9,psig,tt0
+ integer :: i,maxres2
+#ifdef MPI
+ tt0 = MPI_Wtime()
+#else
+ tt0 = tcpu()
+#endif
+!
+! AL 8/17/04 Code adapted from tinker
+!
+! Get the frictional and random terms for stochastic dynamics in the
+! eigenspace of mass-scaled UNRES friction matrix
+!
+ maxres2=2*nres
+ do i = 1, dimen
+ gdt = fricgam(i) * d_time
+!
+! Stochastic dynamics reduces to simple MD for zero friction
+!
+ 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
+!
+! Analytical expressions when friction coefficient is large
+!
+ 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)
+!
+! Use series expansions when friction coefficient is small
+!
+ 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
+!
+! Compute the scaling factors of random terms for the nonzero friction case
+!
+ 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
+!
+! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
+!
+#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 subroutine sd_verlet_p_setup
+!-----------------------------------------------------------------------------
+ subroutine eigtransf1(n,ndim,ab,d,c)
+
+!el implicit none
+ integer :: n,ndim
+ real(kind=8) :: 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 subroutine eigtransf1
+!-----------------------------------------------------------------------------
+ subroutine eigtransf(n,ndim,a,b,d,c)
+
+!el implicit none
+ integer :: n,ndim
+ real(kind=8) :: 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 subroutine eigtransf
+!-----------------------------------------------------------------------------
+ subroutine sd_verlet1
+
+! Applying stochastic velocity Verlet algorithm - step 1 to velocities
+ use energy_data
+! 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'
+!el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
+!el common /stochcalc/ stochforcvec
+ logical :: lprn = .false.
+ real(kind=8) :: ddt1,ddt2
+ integer :: i,j,ind,inres
+
+! write (iout,*) "dc_old"
+! do i=0,nres
+! write (iout,'(i5,3f10.5,5x,3f10.5)')
+! & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
+! 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 .and. itype(i).ne.ntyp1) 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 subroutine sd_verlet1
+!-----------------------------------------------------------------------------
+ subroutine sd_verlet2
+
+! Calculating the adjusted velocities for accelerations
+ use energy_data
+! 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'
+!el real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV !(MAXRES6) maxres6=6*maxres
+ real(kind=8),dimension(6*nres) :: stochforcvecV !(MAXRES6) maxres6=6*maxres
+!el common /stochcalc/ stochforcvec
+!
+ real(kind=8) :: ddt1,ddt2
+ integer :: i,j,ind,inres
+! Compute the stochastic forces which contribute to velocity change
+!
+ 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 .and. itype(i).ne.ntyp1) 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 subroutine sd_verlet2
+!-----------------------------------------------------------------------------
+ subroutine sd_verlet_ciccotti_setup
+
+! Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
+! version
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+ use control, only: tcpu
+ use control_data
+#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'
+ real(kind=8),dimension(6*nres) :: emgdt !(MAXRES6) maxres6=6*maxres
+ real(kind=8) :: pterm,vterm,rho,rhoc,vsig
+ real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
+ prand_vec,vrand_vec1,vrand_vec2 !(MAXRES6) maxres6=6*maxres
+ logical :: lprn = .false.
+ real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
+ real(kind=8) :: ktm,gdt,egdt,tt0
+ integer :: i,maxres2
+#ifdef MPI
+ tt0 = MPI_Wtime()
+#else
+ tt0 = tcpu()
+#endif
+!
+! AL 8/17/04 Code adapted from tinker
+!
+! Get the frictional and random terms for stochastic dynamics in the
+! eigenspace of mass-scaled UNRES friction matrix
+!
+ maxres2=2*nres
+ do i = 1, dimen
+ write (iout,*) "i",i," fricgam",fricgam(i)
+ gdt = fricgam(i) * d_time
+!
+! Stochastic dynamics reduces to simple MD for zero friction
+!
+ 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)
+!
+! Analytical expressions when friction coefficient is large
+!
+ 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)
+!
+! Compute the scaling factors of random terms for the nonzero friction case
+!
+! ktm = 0.5d0*d_time/fricgam(i)
+! psig = dsqrt(ktm*pterm) / fricgam(i)
+! vsig = dsqrt(ktm*vterm)
+! prand_vec(i) = psig*afric_vec(i)
+! 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
+!
+! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
+!
+ 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 subroutine sd_verlet_ciccotti_setup
+!-----------------------------------------------------------------------------
+ subroutine sd_verlet1_ciccotti
+
+! Applying stochastic velocity Verlet algorithm - step 1 to velocities
+! implicit real*8 (a-h,o-z)
+ use energy_data
+! 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'
+!el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
+!el common /stochcalc/ stochforcvec
+ logical :: lprn = .false.
+ real(kind=8) :: ddt1,ddt2
+ integer :: i,j,ind,inres
+! write (iout,*) "dc_old"
+! do i=0,nres
+! write (iout,'(i5,3f10.5,5x,3f10.5)')
+! & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
+! 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 .and. itype(i).ne.ntyp1) 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 .and. itype(i).ne.ntyp1) 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 subroutine sd_verlet1_ciccotti
+!-----------------------------------------------------------------------------
+ subroutine sd_verlet2_ciccotti
+
+! Calculating the adjusted velocities for accelerations
+ use energy_data
+! 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'
+!el real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV !(MAXRES6) maxres6=6*maxres
+ real(kind=8),dimension(6*nres) :: stochforcvecV !(MAXRES6) maxres6=6*maxres
+!el common /stochcalc/ stochforcvec
+ real(kind=8) :: ddt1,ddt2
+ integer :: i,j,ind,inres
+!
+! Compute the stochastic forces which contribute to velocity change
+!
+ 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_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 .and. itype(i).ne.ntyp1) 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 subroutine sd_verlet2_ciccotti
+#endif
+!-----------------------------------------------------------------------------
+! moments.f
+!-----------------------------------------------------------------------------
+ subroutine inertia_tensor
+
+! Calculating the intertia tensor for the entire protein in order to
+! remove the perpendicular components of velocity matrix which cause
+! the molecule to rotate.
+ use comm_gucio
+ use energy_data
+! 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'
+
+ real(kind=8),dimension(3,3) :: Im,Imcp,eigvec,Id
+ real(kind=8),dimension(3) :: pr,eigval,L,vp,vrot
+ real(kind=8) :: M_SC,mag,mag2
+ real(kind=8),dimension(3,0:nres) :: vpp !(3,0:MAXRES)
+ real(kind=8),dimension(3) :: vs_p,pp,incr,v
+ real(kind=8),dimension(3,3) :: pr1,pr2
+
+!el common /gucio/ cm
+ integer :: iti,inres,i,j,k
+ do i=1,3
+ do j=1,3
+ Im(i,j)=0.0d0
+ pr1(i,j)=0.0d0
+ pr2(i,j)=0.0d0
+ enddo
+ L(i)=0.0d0
+ cm(i)=0.0d0
+ vrot(i)=0.0d0
+ enddo
+! calculating the center of the mass of the protein
+ do i=nnt,nct-1
+ do j=1,3
+ cm(j)=cm(j)+c(j,i)+0.5d0*dc(j,i)
+ enddo
+ enddo
+ do j=1,3
+ cm(j)=mp*cm(j)
+ enddo
+ M_SC=0.0d0
+ do i=nnt,nct
+ iti=iabs(itype(i))
+ M_SC=M_SC+msc(iabs(iti))
+ inres=i+nres
+ do j=1,3
+ cm(j)=cm(j)+msc(iabs(iti))*c(j,inres)
+ enddo
+ enddo
+ do j=1,3
+ cm(j)=cm(j)/(M_SC+(nct-nnt)*mp)
+ enddo
+
+ do i=nnt,nct-1
+ do j=1,3
+ pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
+ enddo
+ Im(1,1)=Im(1,1)+mp*(pr(2)*pr(2)+pr(3)*pr(3))
+ Im(1,2)=Im(1,2)-mp*pr(1)*pr(2)
+ Im(1,3)=Im(1,3)-mp*pr(1)*pr(3)
+ Im(2,3)=Im(2,3)-mp*pr(2)*pr(3)
+ Im(2,2)=Im(2,2)+mp*(pr(3)*pr(3)+pr(1)*pr(1))
+ Im(3,3)=Im(3,3)+mp*(pr(1)*pr(1)+pr(2)*pr(2))
+ enddo
+
+ do i=nnt,nct
+ iti=iabs(itype(i))
+ inres=i+nres
+ do j=1,3
+ pr(j)=c(j,inres)-cm(j)
+ enddo
+ Im(1,1)=Im(1,1)+msc(iabs(iti))*(pr(2)*pr(2)+pr(3)*pr(3))
+ Im(1,2)=Im(1,2)-msc(iabs(iti))*pr(1)*pr(2)
+ Im(1,3)=Im(1,3)-msc(iabs(iti))*pr(1)*pr(3)
+ Im(2,3)=Im(2,3)-msc(iabs(iti))*pr(2)*pr(3)
+ Im(2,2)=Im(2,2)+msc(iabs(iti))*(pr(3)*pr(3)+pr(1)*pr(1))
+ Im(3,3)=Im(3,3)+msc(iabs(iti))*(pr(1)*pr(1)+pr(2)*pr(2))
+ enddo
+
+ do i=nnt,nct-1
+ Im(1,1)=Im(1,1)+Ip*(1-dc_norm(1,i)*dc_norm(1,i))* &
+ vbld(i+1)*vbld(i+1)*0.25d0
+ Im(1,2)=Im(1,2)+Ip*(-dc_norm(1,i)*dc_norm(2,i))* &
+ vbld(i+1)*vbld(i+1)*0.25d0
+ Im(1,3)=Im(1,3)+Ip*(-dc_norm(1,i)*dc_norm(3,i))* &
+ vbld(i+1)*vbld(i+1)*0.25d0
+ Im(2,3)=Im(2,3)+Ip*(-dc_norm(2,i)*dc_norm(3,i))* &
+ vbld(i+1)*vbld(i+1)*0.25d0
+ Im(2,2)=Im(2,2)+Ip*(1-dc_norm(2,i)*dc_norm(2,i))* &
+ vbld(i+1)*vbld(i+1)*0.25d0
+ Im(3,3)=Im(3,3)+Ip*(1-dc_norm(3,i)*dc_norm(3,i))* &
+ vbld(i+1)*vbld(i+1)*0.25d0
+ enddo
+
+
+ do i=nnt,nct
+ if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ iti=iabs(itype(i))
+ inres=i+nres
+ Im(1,1)=Im(1,1)+Isc(iti)*(1-dc_norm(1,inres)* &
+ dc_norm(1,inres))*vbld(inres)*vbld(inres)
+ Im(1,2)=Im(1,2)-Isc(iti)*(dc_norm(1,inres)* &
+ dc_norm(2,inres))*vbld(inres)*vbld(inres)
+ Im(1,3)=Im(1,3)-Isc(iti)*(dc_norm(1,inres)* &
+ dc_norm(3,inres))*vbld(inres)*vbld(inres)
+ Im(2,3)=Im(2,3)-Isc(iti)*(dc_norm(2,inres)* &
+ dc_norm(3,inres))*vbld(inres)*vbld(inres)
+ Im(2,2)=Im(2,2)+Isc(iti)*(1-dc_norm(2,inres)* &
+ dc_norm(2,inres))*vbld(inres)*vbld(inres)
+ Im(3,3)=Im(3,3)+Isc(iti)*(1-dc_norm(3,inres)* &
+ dc_norm(3,inres))*vbld(inres)*vbld(inres)
+ endif
+ enddo
+
+ call angmom(cm,L)
+! write(iout,*) "The angular momentum before adjustment:"
+! write(iout,*) (L(j),j=1,3)
+
+ Im(2,1)=Im(1,2)
+ Im(3,1)=Im(1,3)
+ Im(3,2)=Im(2,3)
+
+! Copying the Im matrix for the djacob subroutine
+ do i=1,3
+ do j=1,3
+ Imcp(i,j)=Im(i,j)
+ Id(i,j)=0.0d0
+ enddo
+ enddo
+
+! Finding the eigenvectors and eignvalues of the inertia tensor
+ call djacob(3,3,10000,1.0d-10,Imcp,eigvec,eigval)
+! write (iout,*) "Eigenvalues & Eigenvectors"
+! write (iout,'(5x,3f10.5)') (eigval(i),i=1,3)
+! write (iout,*)
+! do i=1,3
+! write (iout,'(i5,3f10.5)') i,(eigvec(i,j),j=1,3)
+! enddo
+! Constructing the diagonalized matrix
+ do i=1,3
+ if (dabs(eigval(i)).gt.1.0d-15) then
+ Id(i,i)=1.0d0/eigval(i)
+ else
+ Id(i,i)=0.0d0
+ endif
+ enddo
+ do i=1,3
+ do j=1,3
+ Imcp(i,j)=eigvec(j,i)
+ enddo
+ enddo
+ do i=1,3
+ do j=1,3
+ do k=1,3
+ pr1(i,j)=pr1(i,j)+Id(i,k)*Imcp(k,j)
+ enddo
+ enddo
+ enddo
+ do i=1,3
+ do j=1,3
+ do k=1,3
+ pr2(i,j)=pr2(i,j)+eigvec(i,k)*pr1(k,j)
+ enddo
+ enddo
+ enddo
+! Calculating the total rotational velocity of the molecule
+ do i=1,3
+ do j=1,3
+ vrot(i)=vrot(i)+pr2(i,j)*L(j)
+ enddo
+ enddo
+! Resetting the velocities
+ do i=nnt,nct-1
+ call vecpr(vrot(1),dc(1,i),vp)
+ do j=1,3
+ d_t(j,i)=d_t(j,i)-vp(j)
+ enddo
+ enddo
+ do i=nnt,nct
+ if(itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ inres=i+nres
+ call vecpr(vrot(1),dc(1,inres),vp)
+ do j=1,3
+ d_t(j,inres)=d_t(j,inres)-vp(j)
+ enddo
+ endif
+ enddo
+ call angmom(cm,L)
+! write(iout,*) "The angular momentum after adjustment:"
+! write(iout,*) (L(j),j=1,3)
+
+ return
+ end subroutine inertia_tensor
+!-----------------------------------------------------------------------------
+ subroutine angmom(cm,L)
+
+ use energy_data
+! 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'
+
+ real(kind=8),dimension(3) :: L,cm,pr,vp,vrot,incr,v,pp
+ integer :: iti,inres,i,j
+! Calculate the angular momentum
+ do j=1,3
+ L(j)=0.0d0
+ enddo
+ do j=1,3
+ incr(j)=d_t(j,0)
+ enddo
+ do i=nnt,nct-1
+ do j=1,3
+ pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
+ enddo
+ do j=1,3
+ v(j)=incr(j)+0.5d0*d_t(j,i)
+ enddo
+ do j=1,3
+ incr(j)=incr(j)+d_t(j,i)
+ enddo
+ call vecpr(pr(1),v(1),vp)
+ do j=1,3
+ L(j)=L(j)+mp*vp(j)
+ enddo
+ do j=1,3
+ pr(j)=0.5d0*dc(j,i)
+ pp(j)=0.5d0*d_t(j,i)
+ enddo
+ call vecpr(pr(1),pp(1),vp)
+ do j=1,3
+ L(j)=L(j)+Ip*vp(j)
+ enddo
+ enddo
+ do j=1,3
+ incr(j)=d_t(j,0)
+ enddo
+ do i=nnt,nct
+ iti=iabs(itype(i))
+ inres=i+nres
+ do j=1,3
+ pr(j)=c(j,inres)-cm(j)
+ enddo
+ if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ do j=1,3
+ v(j)=incr(j)+d_t(j,inres)
+ enddo
+ else
+ do j=1,3
+ v(j)=incr(j)
+ enddo
+ endif
+ call vecpr(pr(1),v(1),vp)
+! write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),
+! & " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
+ do j=1,3
+ L(j)=L(j)+msc(iabs(iti))*vp(j)
+ enddo
+! write (iout,*) "L",(l(j),j=1,3)
+ if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ do j=1,3
+ v(j)=incr(j)+d_t(j,inres)
+ enddo
+ call vecpr(dc(1,inres),d_t(1,inres),vp)
+ do j=1,3
+ L(j)=L(j)+Isc(iti)*vp(j)
+ enddo
+ endif
+ do j=1,3
+ incr(j)=incr(j)+d_t(j,i)
+ enddo
+ enddo
+ return
+ end subroutine angmom
+!-----------------------------------------------------------------------------
+ subroutine vcm_vel(vcm)
+
+ use energy_data
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.VAR'
+! include 'COMMON.MD'
+! include 'COMMON.CHAIN'
+! include 'COMMON.DERIV'
+! include 'COMMON.GEO'
+! include 'COMMON.LOCAL'
+! include 'COMMON.INTERACT'
+! include 'COMMON.IOUNITS'
+ real(kind=8),dimension(3) :: vcm,vv
+ real(kind=8) :: summas,amas
+ integer :: i,j
+
+ do j=1,3
+ vcm(j)=0.0d0
+ vv(j)=d_t(j,0)
+ enddo
+ summas=0.0d0
+ do i=nnt,nct
+ if (i.lt.nct) then
+ summas=summas+mp
+ do j=1,3
+ vcm(j)=vcm(j)+mp*(vv(j)+0.5d0*d_t(j,i))
+ enddo
+ endif
+ amas=msc(iabs(itype(i)))
+ summas=summas+amas
+ if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ do j=1,3
+ vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres))
+ enddo
+ else
+ do j=1,3
+ vcm(j)=vcm(j)+amas*vv(j)
+ enddo
+ endif
+ do j=1,3
+ vv(j)=vv(j)+d_t(j,i)
+ enddo
+ enddo
+! write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas
+ do j=1,3
+ vcm(j)=vcm(j)/summas
+ enddo
+ return
+ end subroutine vcm_vel
+!-----------------------------------------------------------------------------
+! rattle.F
+!-----------------------------------------------------------------------------
+ subroutine rattle1
+! RATTLE algorithm for velocity Verlet - step 1, UNRES
+! AL 9/24/04
+ use comm_przech
+ use energy_data
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+#ifdef RATTLE
+! 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'
+!el real(kind=8) :: gginv(2*nres,2*nres),&
+!el gdc(3,2*nres,2*nres)
+ real(kind=8) :: dC_uncor(3,2*nres) !,&
+!el real(kind=8) :: Cmat(2*nres,2*nres)
+ real(kind=8) :: x(2*nres),xcorr(3,2*nres) !maxres2=2*maxres
+!el common /przechowalnia/ GGinv,gdc,Cmat,nbond
+!el common /przechowalnia/ nbond
+ integer :: max_rattle = 5
+ logical :: lprn = .false., lprn1 = .false., not_done
+ real(kind=8) :: tol_rattle = 1.0d-5
+
+ integer :: ii,i,j,jj,l,ind,ind1,nres2
+ nres2=2*nres
+
+!el /common/ przechowalnia
+
+ if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
+ if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
+ if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
+!el--------
+ if (lprn) write (iout,*) "RATTLE1"
+ nbond=nct-nnt
+ do i=nnt,nct
+ if (itype(i).ne.10) nbond=nbond+1
+ enddo
+! Make a folded form of the Ginv-matrix
+ ind=0
+ ii=0
+ do i=nnt,nct-1
+ ii=ii+1
+ do j=1,3
+ ind=ind+1
+ ind1=0
+ jj=0
+ do k=nnt,nct-1
+ jj=jj+1
+ do l=1,3
+ ind1=ind1+1
+ if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
+ enddo
+ enddo
+ do k=nnt,nct
+ if (itype(k).ne.10) then
+ jj=jj+1
+ do l=1,3
+ ind1=ind1+1
+ if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
+ enddo
+ endif
+ enddo
+ enddo
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ ii=ii+1
+ do j=1,3
+ ind=ind+1
+ ind1=0
+ jj=0
+ do k=nnt,nct-1
+ jj=jj+1
+ do l=1,3
+ ind1=ind1+1
+ if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
+ enddo
+ enddo
+ do k=nnt,nct
+ if (itype(k).ne.10) then
+ jj=jj+1
+ do l=1,3
+ ind1=ind1+1
+ if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
+ enddo
+ endif
+ enddo
+ enddo
+ endif
+ enddo
+ if (lprn1) then
+ write (iout,*) "Matrix GGinv"
+ call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
+ endif
+ not_done=.true.
+ iter=0
+ do while (not_done)
+ iter=iter+1
+ if (iter.gt.max_rattle) then
+ write (iout,*) "Error - too many iterations in RATTLE."
+ stop
+ endif
+! Calculate the matrix C = GG**(-1) dC_old o dC
+ ind1=0
+ do i=nnt,nct-1
+ ind1=ind1+1
+ do j=1,3
+ dC_uncor(j,ind1)=dC(j,i)
+ enddo
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ ind1=ind1+1
+ do j=1,3
+ dC_uncor(j,ind1)=dC(j,i+nres)
+ enddo
+ endif
+ enddo
+ do i=1,nbond
+ ind=0
+ do k=nnt,nct-1
+ ind=ind+1
+ do j=1,3
+ gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
+ enddo
+ enddo
+ do k=nnt,nct
+ if (itype(k).ne.10) then
+ ind=ind+1
+ do j=1,3
+ gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
+ enddo
+ endif
+ enddo
+ enddo
+! Calculate deviations from standard virtual-bond lengths
+ ind=0
+ do i=nnt,nct-1
+ ind=ind+1
+ x(ind)=vbld(i+1)**2-vbl**2
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ ind=ind+1
+ x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i))**2
+ endif
+ enddo
+ if (lprn) then
+ write (iout,*) "Coordinates and violations"
+ do i=1,nbond
+ write(iout,'(i5,3f10.5,5x,e15.5)') &
+ i,(dC_uncor(j,i),j=1,3),x(i)
+ enddo
+ write (iout,*) "Velocities and violations"
+ ind=0
+ do i=nnt,nct-1
+ ind=ind+1
+ write (iout,'(2i5,3f10.5,5x,e15.5)') &
+ i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ ind=ind+1
+ write (iout,'(2i5,3f10.5,5x,e15.5)') &
+ i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
+ scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
+ endif
+ enddo
+! write (iout,*) "gdc"
+! do i=1,nbond
+! write (iout,*) "i",i
+! do j=1,nbond
+! write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
+! enddo
+! enddo
+ endif
+ xmax=dabs(x(1))
+ do i=2,nbond
+ if (dabs(x(i)).gt.xmax) then
+ xmax=dabs(x(i))
+ endif
+ enddo
+ if (xmax.lt.tol_rattle) then
+ not_done=.false.
+ goto 100
+ endif
+! Calculate the matrix of the system of equations
+ do i=1,nbond
+ do j=1,nbond
+ Cmat(i,j)=0.0d0
+ do k=1,3
+ Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
+ enddo
+ enddo
+ enddo
+ if (lprn1) then
+ write (iout,*) "Matrix Cmat"
+ call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
+ endif
+ call gauss(Cmat,X,MAXRES2,nbond,1,*10)
+! Add constraint term to positions
+ ind=0
+ do i=nnt,nct-1
+ ind=ind+1
+ do j=1,3
+ xx=0.0d0
+ do ii=1,nbond
+ xx = xx+x(ii)*gdc(j,ind,ii)
+ enddo
+ xx=0.5d0*xx
+ dC(j,i)=dC(j,i)-xx
+ d_t_new(j,i)=d_t_new(j,i)-xx/d_time
+ enddo
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ ind=ind+1
+ do j=1,3
+ xx=0.0d0
+ do ii=1,nbond
+ xx = xx+x(ii)*gdc(j,ind,ii)
+ enddo
+ xx=0.5d0*xx
+ dC(j,i+nres)=dC(j,i+nres)-xx
+ d_t_new(j,i+nres)=d_t_new(j,i+nres)-xx/d_time
+ enddo
+ endif
+ enddo
+! Rebuild the chain using the new coordinates
+ call chainbuild_cart
+ if (lprn) then
+ write (iout,*) "New coordinates, Lagrange multipliers,",&
+ " and differences between actual and standard bond lengths"
+ ind=0
+ do i=nnt,nct-1
+ ind=ind+1
+ xx=vbld(i+1)**2-vbl**2
+ write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
+ i,(dC(j,i),j=1,3),x(ind),xx
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ ind=ind+1
+ xx=vbld(i+nres)**2-vbldsc0(1,itype(i))**2
+ write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
+ i,(dC(j,i+nres),j=1,3),x(ind),xx
+ endif
+ enddo
+ write (iout,*) "Velocities and violations"
+ ind=0
+ do i=nnt,nct-1
+ ind=ind+1
+ write (iout,'(2i5,3f10.5,5x,e15.5)') &
+ i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ ind=ind+1
+ write (iout,'(2i5,3f10.5,5x,e15.5)') &
+ i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
+ scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
+ endif
+ enddo
+ endif
+ enddo
+ 100 continue
+ return
+ 10 write (iout,*) "Error - singularity in solving the system",&
+ " of equations for Lagrange multipliers."
+ stop
+#else
+ write (iout,*) &
+ "RATTLE inactive; use -DRATTLE switch at compile time."
+ stop
+#endif
+ end subroutine rattle1
+!-----------------------------------------------------------------------------
+ subroutine rattle2
+! RATTLE algorithm for velocity Verlet - step 2, UNRES
+! AL 9/24/04
+ use comm_przech
+ use energy_data
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+#ifdef RATTLE
+! 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'
+!el real(kind=8) :: gginv(2*nres,2*nres),&
+!el gdc(3,2*nres,2*nres)
+ real(kind=8) :: dC_uncor(3,2*nres) !,&
+!el Cmat(2*nres,2*nres)
+ real(kind=8) :: x(2*nres) !maxres2=2*maxres
+!el common /przechowalnia/ GGinv,gdc,Cmat,nbond
+!el common /przechowalnia/ nbond
+ integer :: max_rattle = 5
+ logical :: lprn = .false., lprn1 = .false., not_done
+ real(kind=8) :: tol_rattle = 1.0d-5
+ integer :: nres2
+ nres2=2*nres
+
+!el /common/ przechowalnia
+ if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
+ if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
+ if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
+!el--------
+ if (lprn) write (iout,*) "RATTLE2"
+ if (lprn) write (iout,*) "Velocity correction"
+! Calculate the matrix G dC
+ do i=1,nbond
+ ind=0
+ do k=nnt,nct-1
+ ind=ind+1
+ do j=1,3
+ gdc(j,i,ind)=GGinv(i,ind)*dC(j,k)
+ enddo
+ enddo
+ do k=nnt,nct
+ if (itype(k).ne.10) then
+ ind=ind+1
+ do j=1,3
+ gdc(j,i,ind)=GGinv(i,ind)*dC(j,k+nres)
+ enddo
+ endif
+ enddo
+ enddo
+! if (lprn) then
+! write (iout,*) "gdc"
+! do i=1,nbond
+! write (iout,*) "i",i
+! do j=1,nbond
+! write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
+! enddo
+! enddo
+! endif
+! Calculate the matrix of the system of equations
+ ind=0
+ do i=nnt,nct-1
+ ind=ind+1
+ do j=1,nbond
+ Cmat(ind,j)=0.0d0
+ do k=1,3
+ Cmat(ind,j)=Cmat(ind,j)+dC(k,i)*gdc(k,ind,j)
+ enddo
+ enddo
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ ind=ind+1
+ do j=1,nbond
+ Cmat(ind,j)=0.0d0
+ do k=1,3
+ Cmat(ind,j)=Cmat(ind,j)+dC(k,i+nres)*gdc(k,ind,j)
+ enddo
+ enddo
+ endif
+ enddo
+! Calculate the scalar product dC o d_t_new
+ ind=0
+ do i=nnt,nct-1
+ ind=ind+1
+ x(ind)=scalar(d_t(1,i),dC(1,i))
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ ind=ind+1
+ x(ind)=scalar(d_t(1,i+nres),dC(1,i+nres))
+ endif
+ enddo
+ if (lprn) then
+ write (iout,*) "Velocities and violations"
+ ind=0
+ do i=nnt,nct-1
+ ind=ind+1
+ write (iout,'(2i5,3f10.5,5x,e15.5)') &
+ i,ind,(d_t(j,i),j=1,3),x(ind)
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ ind=ind+1
+ write (iout,'(2i5,3f10.5,5x,e15.5)') &
+ i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind)
+ endif
+ enddo
+ endif
+ xmax=dabs(x(1))
+ do i=2,nbond
+ if (dabs(x(i)).gt.xmax) then
+ xmax=dabs(x(i))
+ endif
+ enddo
+ if (xmax.lt.tol_rattle) then
+ not_done=.false.
+ goto 100
+ endif
+ if (lprn1) then
+ write (iout,*) "Matrix Cmat"
+ call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
+ endif
+ call gauss(Cmat,X,MAXRES2,nbond,1,*10)
+! Add constraint term to velocities
+ ind=0
+ do i=nnt,nct-1
+ ind=ind+1
+ do j=1,3
+ xx=0.0d0
+ do ii=1,nbond
+ xx = xx+x(ii)*gdc(j,ind,ii)
+ enddo
+ d_t(j,i)=d_t(j,i)-xx
+ enddo
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ ind=ind+1
+ do j=1,3
+ xx=0.0d0
+ do ii=1,nbond
+ xx = xx+x(ii)*gdc(j,ind,ii)
+ enddo
+ d_t(j,i+nres)=d_t(j,i+nres)-xx
+ enddo
+ endif
+ enddo
+ if (lprn) then
+ write (iout,*) &
+ "New velocities, Lagrange multipliers violations"
+ ind=0
+ do i=nnt,nct-1
+ ind=ind+1
+ if (lprn) write (iout,'(2i5,3f10.5,5x,2e15.5)') &
+ i,ind,(d_t(j,i),j=1,3),x(ind),scalar(d_t(1,i),dC(1,i))
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ ind=ind+1
+ write (iout,'(2i5,3f10.5,5x,2e15.5)') &
+ i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind),&
+ scalar(d_t(1,i+nres),dC(1,i+nres))
+ endif
+ enddo
+ endif
+ 100 continue
+ return
+ 10 write (iout,*) "Error - singularity in solving the system",&
+ " of equations for Lagrange multipliers."
+ stop
+#else
+ write (iout,*) &
+ "RATTLE inactive; use -DRATTLE option at compile time."
+ stop
+#endif
+ end subroutine rattle2
+!-----------------------------------------------------------------------------
+ subroutine rattle_brown
+! RATTLE/LINCS algorithm for Brownian dynamics, UNRES
+! AL 9/24/04
+ use comm_przech
+ use energy_data
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+#ifdef RATTLE
+! 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'
+!el real(kind=8) :: gginv(2*nres,2*nres),&
+!el gdc(3,2*nres,2*nres)
+ real(kind=8) :: dC_uncor(3,2*nres) !,&
+!el real(kind=8) :: Cmat(2*nres,2*nres)
+ real(kind=8) :: x(2*nres) !maxres2=2*maxres
+!el common /przechowalnia/ GGinv,gdc,Cmat,nbond
+!el common /przechowalnia/ nbond
+ integer :: max_rattle = 5
+ logical :: lprn = .true., lprn1 = .true., not_done
+ real(kind=8) :: tol_rattle = 1.0d-5
+ integer :: nres2
+ nres2=2*nres
+
+!el /common/ przechowalnia
+ if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
+ if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
+ if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
+!el--------
+
+ if (lprn) write (iout,*) "RATTLE_BROWN"
+ nbond=nct-nnt
+ do i=nnt,nct
+ if (itype(i).ne.10) nbond=nbond+1
+ enddo
+! Make a folded form of the Ginv-matrix
+ ind=0
+ ii=0
+ do i=nnt,nct-1
+ ii=ii+1
+ do j=1,3
+ ind=ind+1
+ ind1=0
+ jj=0
+ do k=nnt,nct-1
+ jj=jj+1
+ do l=1,3
+ ind1=ind1+1
+ if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
+ enddo
+ enddo
+ do k=nnt,nct
+ if (itype(k).ne.10) then
+ jj=jj+1
+ do l=1,3
+ ind1=ind1+1
+ if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
+ enddo
+ endif
+ enddo
+ enddo
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ ii=ii+1
+ do j=1,3
+ ind=ind+1
+ ind1=0
+ jj=0
+ do k=nnt,nct-1
+ jj=jj+1
+ do l=1,3
+ ind1=ind1+1
+ if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
+ enddo
+ enddo
+ do k=nnt,nct
+ if (itype(k).ne.10) then
+ jj=jj+1
+ do l=1,3
+ ind1=ind1+1
+ if (j.eq.1 .and. l.eq.1)GGinv(ii,jj)=fricmat(ind,ind1)
+ enddo
+ endif
+ enddo
+ enddo
+ endif
+ enddo
+ if (lprn1) then
+ write (iout,*) "Matrix GGinv"
+ call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
+ endif
+ not_done=.true.
+ iter=0
+ do while (not_done)
+ iter=iter+1
+ if (iter.gt.max_rattle) then
+ write (iout,*) "Error - too many iterations in RATTLE."
+ stop
+ endif
+! Calculate the matrix C = GG**(-1) dC_old o dC
+ ind1=0
+ do i=nnt,nct-1
+ ind1=ind1+1
+ do j=1,3
+ dC_uncor(j,ind1)=dC(j,i)
+ enddo
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ ind1=ind1+1
+ do j=1,3
+ dC_uncor(j,ind1)=dC(j,i+nres)
+ enddo
+ endif
+ enddo
+ do i=1,nbond
+ ind=0
+ do k=nnt,nct-1
+ ind=ind+1
+ do j=1,3
+ gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
+ enddo
+ enddo
+ do k=nnt,nct
+ if (itype(k).ne.10) then
+ ind=ind+1
+ do j=1,3
+ gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
+ enddo
+ endif
+ enddo
+ enddo
+! Calculate deviations from standard virtual-bond lengths
+ ind=0
+ do i=nnt,nct-1
+ ind=ind+1
+ x(ind)=vbld(i+1)**2-vbl**2
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ ind=ind+1
+ x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i))**2
+ endif
+ enddo
+ if (lprn) then
+ write (iout,*) "Coordinates and violations"
+ do i=1,nbond
+ write(iout,'(i5,3f10.5,5x,e15.5)') &
+ i,(dC_uncor(j,i),j=1,3),x(i)
+ enddo
+ write (iout,*) "Velocities and violations"
+ ind=0
+ do i=nnt,nct-1
+ ind=ind+1
+ write (iout,'(2i5,3f10.5,5x,e15.5)') &
+ i,ind,(d_t(j,i),j=1,3),scalar(d_t(1,i),dC_old(1,i))
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ ind=ind+1
+ write (iout,'(2i5,3f10.5,5x,e15.5)') &
+ i+nres,ind,(d_t(j,i+nres),j=1,3),&
+ scalar(d_t(1,i+nres),dC_old(1,i+nres))
+ endif
+ enddo
+ write (iout,*) "gdc"
+ do i=1,nbond
+ write (iout,*) "i",i
+ do j=1,nbond
+ write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
+ enddo
+ enddo
+ endif
+ xmax=dabs(x(1))
+ do i=2,nbond
+ if (dabs(x(i)).gt.xmax) then
+ xmax=dabs(x(i))
+ endif
+ enddo
+ if (xmax.lt.tol_rattle) then
+ not_done=.false.
+ goto 100
+ endif
+! Calculate the matrix of the system of equations
+ do i=1,nbond
+ do j=1,nbond
+ Cmat(i,j)=0.0d0
+ do k=1,3
+ Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
+ enddo
+ enddo
+ enddo
+ if (lprn1) then
+ write (iout,*) "Matrix Cmat"
+ call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
+ endif
+ call gauss(Cmat,X,MAXRES2,nbond,1,*10)
+! Add constraint term to positions
+ ind=0
+ do i=nnt,nct-1
+ ind=ind+1
+ do j=1,3
+ xx=0.0d0
+ do ii=1,nbond
+ xx = xx+x(ii)*gdc(j,ind,ii)
+ enddo
+ xx=-0.5d0*xx
+ d_t(j,i)=d_t(j,i)+xx/d_time
+ dC(j,i)=dC(j,i)+xx
+ enddo
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ ind=ind+1
+ do j=1,3
+ xx=0.0d0
+ do ii=1,nbond
+ xx = xx+x(ii)*gdc(j,ind,ii)
+ enddo
+ xx=-0.5d0*xx
+ d_t(j,i+nres)=d_t(j,i+nres)+xx/d_time
+ dC(j,i+nres)=dC(j,i+nres)+xx
+ enddo
+ endif
+ enddo
+! Rebuild the chain using the new coordinates
+ call chainbuild_cart
+ if (lprn) then
+ write (iout,*) "New coordinates, Lagrange multipliers,",&
+ " and differences between actual and standard bond lengths"
+ ind=0
+ do i=nnt,nct-1
+ ind=ind+1
+ xx=vbld(i+1)**2-vbl**2
+ write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
+ i,(dC(j,i),j=1,3),x(ind),xx
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ ind=ind+1
+ xx=vbld(i+nres)**2-vbldsc0(1,itype(i))**2
+ write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
+ i,(dC(j,i+nres),j=1,3),x(ind),xx
+ endif
+ enddo
+ write (iout,*) "Velocities and violations"
+ ind=0
+ do i=nnt,nct-1
+ ind=ind+1
+ write (iout,'(2i5,3f10.5,5x,e15.5)') &
+ i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ ind=ind+1
+ write (iout,'(2i5,3f10.5,5x,e15.5)') &
+ i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
+ scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
+ endif
+ enddo
+ endif
+ enddo
+ 100 continue
+ return
+ 10 write (iout,*) "Error - singularity in solving the system",&
+ " of equations for Lagrange multipliers."
+ stop
+#else
+ write (iout,*) &
+ "RATTLE inactive; use -DRATTLE option at compile time"
+ stop
+#endif
+ end subroutine rattle_brown
+!-----------------------------------------------------------------------------
+! stochfric.F
+!-----------------------------------------------------------------------------
+ subroutine friction_force
+
+ use energy_data
+ use REMD_data
+ use comm_syfek
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.VAR'
+! include 'COMMON.CHAIN'
+! include 'COMMON.DERIV'
+! include 'COMMON.GEO'
+! include 'COMMON.LOCAL'
+! include 'COMMON.INTERACT'
+! include 'COMMON.MD'
+!#ifndef LANG0
+! include 'COMMON.LANGEVIN'
+!#else
+! include 'COMMON.LANGEVIN.lang0'
+!#endif
+! include 'COMMON.IOUNITS'
+!el real(kind=8),dimension(6*nres) :: gamvec !(MAXRES6) maxres6=6*maxres
+!el common /syfek/ gamvec
+ real(kind=8) :: vv(3),vvtot(3,nres),v_work(6*nres) !,&
+!el ginvfric(2*nres,2*nres) !maxres2=2*maxres
+!el common /przechowalnia/ ginvfric
+
+ logical :: lprn = .false., checkmode = .false.
+ integer :: i,j,ind,k,nres2,nres6
+ nres2=2*nres
+ nres6=6*nres
+
+ if(.not.allocated(gamvec)) allocate(gamvec(nres6)) !(MAXRES6)
+ if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
+ do i=0,nres2
+ do j=1,3
+ friction(j,i)=0.0d0
+ enddo
+ enddo
+
+ do j=1,3
+ d_t_work(j)=d_t(j,0)
+ enddo
+ ind=3
+ do i=nnt,nct-1
+ do j=1,3
+ d_t_work(ind+j)=d_t(j,i)
+ enddo
+ ind=ind+3
+ enddo
+ do i=nnt,nct
+ if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
+ do j=1,3
+ d_t_work(ind+j)=d_t(j,i+nres)
+ enddo
+ ind=ind+3
+ endif
+ enddo
+
+ call fricmat_mult(d_t_work,fric_work)
+
+ if (.not.checkmode) return
+
+ if (lprn) then
+ write (iout,*) "d_t_work and fric_work"
+ do i=1,3*dimen
+ write (iout,'(i3,2e15.5)') i,d_t_work(i),fric_work(i)
+ enddo
+ endif
+ do j=1,3
+ friction(j,0)=fric_work(j)
+ enddo
+ ind=3
+ do i=nnt,nct-1
+ do j=1,3
+ friction(j,i)=fric_work(ind+j)
+ enddo
+ ind=ind+3
+ enddo
+ do i=nnt,nct
+ if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
+ do j=1,3
+ friction(j,i+nres)=fric_work(ind+j)
+ enddo
+ ind=ind+3
+ endif
+ enddo
+ if (lprn) then
+ write(iout,*) "Friction backbone"
+ do i=0,nct-1
+ write(iout,'(i5,3e15.5,5x,3e15.5)') &
+ i,(friction(j,i),j=1,3),(d_t(j,i),j=1,3)
+ enddo
+ write(iout,*) "Friction side chain"
+ do i=nnt,nct
+ write(iout,'(i5,3e15.5,5x,3e15.5)') &
+ i,(friction(j,i+nres),j=1,3),(d_t(j,i+nres),j=1,3)
+ enddo
+ endif
+ if (lprn) then
+ do j=1,3
+ vv(j)=d_t(j,0)
+ enddo
+ do i=nnt,nct
+ do j=1,3
+ vvtot(j,i)=vv(j)+0.5d0*d_t(j,i)
+ vvtot(j,i+nres)=vv(j)+d_t(j,i+nres)
+ vv(j)=vv(j)+d_t(j,i)
+ enddo
+ enddo
+ write (iout,*) "vvtot backbone and sidechain"
+ do i=nnt,nct
+ write (iout,'(i5,3e15.5,5x,3e15.5)') i,(vvtot(j,i),j=1,3),&
+ (vvtot(j,i+nres),j=1,3)
+ enddo
+ ind=0
+ do i=nnt,nct-1
+ do j=1,3
+ v_work(ind+j)=vvtot(j,i)
+ enddo
+ ind=ind+3
+ enddo
+ do i=nnt,nct
+ do j=1,3
+ v_work(ind+j)=vvtot(j,i+nres)
+ enddo
+ ind=ind+3
+ enddo
+ write (iout,*) "v_work gamvec and site-based friction forces"
+ do i=1,dimen1
+ write (iout,'(i5,3e15.5)') i,v_work(i),gamvec(i),&
+ gamvec(i)*v_work(i)
+ enddo
+! do i=1,dimen
+! fric_work1(i)=0.0d0
+! do j=1,dimen1
+! fric_work1(i)=fric_work1(i)-A(j,i)*gamvec(j)*v_work(j)
+! enddo
+! enddo
+! write (iout,*) "fric_work and fric_work1"
+! do i=1,dimen
+! write (iout,'(i5,2e15.5)') i,fric_work(i),fric_work1(i)
+! enddo
+ do i=1,dimen
+ do j=1,dimen
+ ginvfric(i,j)=0.0d0
+ do k=1,dimen
+ ginvfric(i,j)=ginvfric(i,j)+ginv(i,k)*fricmat(k,j)
+ enddo
+ enddo
+ enddo
+ write (iout,*) "ginvfric"
+ do i=1,dimen
+ write (iout,'(i5,100f8.3)') i,(ginvfric(i,j),j=1,dimen)
+ enddo
+ write (iout,*) "symmetry check"
+ do i=1,dimen
+ do j=1,i-1
+ write (iout,*) i,j,ginvfric(i,j)-ginvfric(j,i)
+ enddo
+ enddo
+ endif
+ return
+ end subroutine friction_force
+!-----------------------------------------------------------------------------
+ subroutine setup_fricmat
+
+! use MPI
+ use energy_data
+ use control_data, only:time_Bcast
+ use control, only:tcpu
+ use comm_syfek
+! implicit real*8 (a-h,o-z)
+#ifdef MPI
+ use MPI_data
+ include 'mpif.h'
+ real(kind=8) :: time00
+#endif
+! include 'DIMENSIONS'
+! include 'COMMON.VAR'
+! include 'COMMON.CHAIN'
+! include 'COMMON.DERIV'
+! include 'COMMON.GEO'
+! include 'COMMON.LOCAL'
+! include 'COMMON.INTERACT'
+! include 'COMMON.MD'
+! include 'COMMON.SETUP'
+! include 'COMMON.TIME1'
+! integer licznik /0/
+! save licznik
+!#ifndef LANG0
+! include 'COMMON.LANGEVIN'
+!#else
+! include 'COMMON.LANGEVIN.lang0'
+!#endif
+! include 'COMMON.IOUNITS'
+ integer :: IERROR
+ integer :: i,j,ind,ind1,m
+ logical :: lprn = .false.
+ real(kind=8) :: dtdi !el ,gamvec(2*nres)
+!el real(kind=8),dimension(2*nres,2*nres) :: ginvfric,fcopy
+ real(kind=8),dimension(2*nres,2*nres) :: fcopy
+!el real(kind=8),dimension(2*nres*(2*nres+1)/2) :: Ghalf !(mmaxres2) (mmaxres2=(maxres2*(maxres2+1)/2))
+!el common /syfek/ gamvec
+ real(kind=8) :: work(8*2*nres)
+ integer :: iwork(2*nres)
+!el common /przechowalnia/ ginvfric,Ghalf,fcopy
+ integer :: ii,iti,k,l,nzero,nres2,nres6,ierr
+#ifdef MPI
+ if (fg_rank.ne.king) goto 10
+#endif
+ nres2=2*nres
+ nres6=6*nres
+
+ if(.not.allocated(gamvec)) allocate(gamvec(nres2)) !(MAXRES2)
+ if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
+!el if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
+!el allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
+ if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !maxres2=2*maxres
+
+!el if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
+! Zeroing out fricmat
+ do i=1,dimen
+ do j=1,dimen
+ fricmat(i,j)=0.0d0
+ enddo
+ enddo
+! Load the friction coefficients corresponding to peptide groups
+ ind1=0
+ do i=nnt,nct-1
+ ind1=ind1+1
+ gamvec(ind1)=gamp
+ enddo
+! Load the friction coefficients corresponding to side chains
+ m=nct-nnt
+ ind=0
+ gamsc(ntyp1)=1.0d0
+ do i=nnt,nct
+ ind=ind+1
+ ii = ind+m
+ iti=itype(i)
+ gamvec(ii)=gamsc(iabs(iti))
+ enddo
+ if (surfarea) call sdarea(gamvec)
+! if (lprn) then
+! write (iout,*) "Matrix A and vector gamma"
+! do i=1,dimen1
+! write (iout,'(i2,$)') i
+! do j=1,dimen
+! write (iout,'(f4.1,$)') A(i,j)
+! enddo
+! write (iout,'(f8.3)') gamvec(i)
+! enddo
+! endif
+ if (lprn) then
+ write (iout,*) "Vector gamvec"
+ do i=1,dimen1
+ write (iout,'(i5,f10.5)') i, gamvec(i)
+ enddo
+ endif
+
+! The friction matrix
+ do k=1,dimen
+ do i=1,dimen
+ dtdi=0.0d0
+ do j=1,dimen1
+ dtdi=dtdi+A(j,k)*A(j,i)*gamvec(j)
+ enddo
+ fricmat(k,i)=dtdi
+ enddo
+ enddo
+
+ if (lprn) then
+ write (iout,'(//a)') "Matrix fricmat"
+ call matout2(dimen,dimen,nres2,nres2,fricmat)
+ endif
+ if (lang.eq.2 .or. lang.eq.3) then
+! Mass-scale the friction matrix if non-direct integration will be performed
+ do i=1,dimen
+ do j=1,dimen
+ Ginvfric(i,j)=0.0d0
+ do k=1,dimen
+ do l=1,dimen
+ Ginvfric(i,j)=Ginvfric(i,j)+ &
+ Gsqrm(i,k)*Gsqrm(l,j)*fricmat(k,l)
+ enddo
+ enddo
+ enddo
+ enddo
+! Diagonalize the friction matrix
+ ind=0
+ do i=1,dimen
+ do j=1,i
+ ind=ind+1
+ Ghalf(ind)=Ginvfric(i,j)
+ enddo
+ enddo
+ call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
+ ierr,iwork)
+ if (lprn) then
+ write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
+ " mass-scaled friction matrix"
+ call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
+ endif
+! Precompute matrices for tinker stochastic integrator
+#ifndef LANG0
+ do i=1,dimen
+ do j=1,dimen
+ mt1(i,j)=0.0d0
+ mt2(i,j)=0.0d0
+ do k=1,dimen
+ mt1(i,j)=mt1(i,j)+fricvec(k,i)*gsqrm(k,j)
+ mt2(i,j)=mt2(i,j)+fricvec(k,i)*gsqrp(k,j)
+ enddo
+ mt3(j,i)=mt1(i,j)
+ enddo
+ enddo
+#endif
+ else if (lang.eq.4) then
+! Diagonalize the friction matrix
+ ind=0
+ do i=1,dimen
+ do j=1,i
+ ind=ind+1
+ Ghalf(ind)=fricmat(i,j)
+ enddo
+ enddo
+ call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
+ ierr,iwork)
+ if (lprn) then
+ write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
+ " friction matrix"
+ call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
+ endif
+! Determine the number of zero eigenvalues of the friction matrix
+ nzero=max0(dimen-dimen1,0)
+! do while (fricgam(nzero+1).le.1.0d-5 .and. nzero.lt.dimen)
+! nzero=nzero+1
+! enddo
+ write (iout,*) "Number of zero eigenvalues:",nzero
+ do i=1,dimen
+ do j=1,dimen
+ fricmat(i,j)=0.0d0
+ do k=nzero+1,dimen
+ fricmat(i,j)=fricmat(i,j) &
+ +fricvec(i,k)*fricvec(j,k)/fricgam(k)
+ enddo
+ enddo
+ enddo
+ if (lprn) then
+ write (iout,'(//a)') "Generalized inverse of fricmat"
+ call matout(dimen,dimen,nres6,nres6,fricmat)
+ endif
+ endif
+#ifdef MPI
+ 10 continue
+ if (nfgtasks.gt.1) then
+ if (fg_rank.eq.0) then
+! The matching BROADCAST for fg processors is called in ERGASTULUM
+#ifdef MPI
+ time00=MPI_Wtime()
+#else
+ time00=tcpu()
+#endif
+ call MPI_Bcast(10,1,MPI_INTEGER,king,FG_COMM,IERROR)
+#ifdef MPI
+ time_Bcast=time_Bcast+MPI_Wtime()-time00
+#else
+ time_Bcast=time_Bcast+tcpu()-time00
+#endif
+! print *,"Processor",myrank,
+! & " BROADCAST iorder in SETUP_FRICMAT"
+ endif
+! licznik=licznik+1
+ write (iout,*) "setup_fricmat licznik"!,licznik !sp
+#ifdef MPI
+ time00=MPI_Wtime()
+#else
+ time00=tcpu()
+#endif
+! Scatter the friction matrix
+ call MPI_Scatterv(fricmat(1,1),nginv_counts(0),&
+ nginv_start(0),MPI_DOUBLE_PRECISION,fcopy(1,1),&
+ myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
+#ifdef TIMING
+#ifdef MPI
+ time_scatter=time_scatter+MPI_Wtime()-time00
+ time_scatter_fmat=time_scatter_fmat+MPI_Wtime()-time00
+#else
+ time_scatter=time_scatter+tcpu()-time00
+ time_scatter_fmat=time_scatter_fmat+tcpu()-time00
+#endif
+#endif
+ do i=1,dimen
+ do j=1,2*my_ng_count
+ fricmat(j,i)=fcopy(i,j)
+ enddo
+ enddo
+! write (iout,*) "My chunk of fricmat"
+! call MATOUT2(my_ng_count,dimen,maxres2,maxres2,fcopy)
+ endif
+#endif
+ return
+ end subroutine setup_fricmat
+!-----------------------------------------------------------------------------
+ subroutine stochastic_force(stochforcvec)
+
+ use energy_data
+ use random, only:anorm_distr
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+ use control, only: tcpu
+ use control_data
+#ifdef MPI
+ include 'mpif.h'
+#endif
+! include 'COMMON.VAR'
+! include 'COMMON.CHAIN'
+! include 'COMMON.DERIV'
+! include 'COMMON.GEO'
+! include 'COMMON.LOCAL'
+! include 'COMMON.INTERACT'
+! include 'COMMON.MD'
+! include 'COMMON.TIME1'
+!#ifndef LANG0
+! include 'COMMON.LANGEVIN'
+!#else
+! include 'COMMON.LANGEVIN.lang0'
+!#endif
+! include 'COMMON.IOUNITS'
+
+ real(kind=8) :: x,sig,lowb,highb
+ real(kind=8) :: ff(3),force(3,0:2*nres),zeta2,lowb2
+ real(kind=8) :: highb2,sig2,forcvec(6*nres),stochforcvec(6*nres)
+ real(kind=8) :: time00
+ logical :: lprn = .false.
+ integer :: i,j,ind
+
+ do i=0,2*nres
+ do j=1,3
+ stochforc(j,i)=0.0d0
+ enddo
+ enddo
+ x=0.0d0
+
+#ifdef MPI
+ time00=MPI_Wtime()
+#else
+ time00=tcpu()
+#endif
+! Compute the stochastic forces acting on bodies. Store in force.
+ do i=nnt,nct-1
+ sig=stdforcp(i)
+ lowb=-5*sig
+ highb=5*sig
+ do j=1,3
+ force(j,i)=anorm_distr(x,sig,lowb,highb)
+ enddo
+ enddo
+ do i=nnt,nct
+ sig2=stdforcsc(i)
+ lowb2=-5*sig2
+ highb2=5*sig2
+ do j=1,3
+ force(j,i+nres)=anorm_distr(x,sig2,lowb2,highb2)
+ enddo
+ enddo
+#ifdef MPI
+ time_fsample=time_fsample+MPI_Wtime()-time00
+#else
+ time_fsample=time_fsample+tcpu()-time00
+#endif
+! Compute the stochastic forces acting on virtual-bond vectors.
+ do j=1,3
+ ff(j)=0.0d0
+ enddo
+ do i=nct-1,nnt,-1
+ do j=1,3
+ stochforc(j,i)=ff(j)+0.5d0*force(j,i)
+ enddo
+ do j=1,3
+ ff(j)=ff(j)+force(j,i)
+ enddo
+ if (itype(i+1).ne.ntyp1) then
+ do j=1,3
+ stochforc(j,i)=stochforc(j,i)+force(j,i+nres+1)
+ ff(j)=ff(j)+force(j,i+nres+1)
+ enddo
+ endif
+ enddo
+ do j=1,3
+ stochforc(j,0)=ff(j)+force(j,nnt+nres)
+ enddo
+ do i=nnt,nct
+ if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
+ do j=1,3
+ stochforc(j,i+nres)=force(j,i+nres)
+ enddo
+ endif
+ enddo
+
+ do j=1,3
+ stochforcvec(j)=stochforc(j,0)
+ enddo
+ ind=3
+ do i=nnt,nct-1
+ do j=1,3
+ stochforcvec(ind+j)=stochforc(j,i)
+ enddo
+ ind=ind+3
+ enddo
+ do i=nnt,nct
+ if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
+ do j=1,3
+ stochforcvec(ind+j)=stochforc(j,i+nres)
+ enddo
+ ind=ind+3
+ endif
+ enddo
+ if (lprn) then
+ write (iout,*) "stochforcvec"
+ do i=1,3*dimen
+ write(iout,'(i5,e15.5)') i,stochforcvec(i)
+ enddo
+ write(iout,*) "Stochastic forces backbone"
+ do i=0,nct-1
+ write(iout,'(i5,3e15.5)') i,(stochforc(j,i),j=1,3)
+ enddo
+ write(iout,*) "Stochastic forces side chain"
+ do i=nnt,nct
+ write(iout,'(i5,3e15.5)') &
+ i,(stochforc(j,i+nres),j=1,3)
+ enddo
+ endif
+
+ if (lprn) then
+
+ ind=0
+ do i=nnt,nct-1
+ write (iout,*) i,ind
+ do j=1,3
+ forcvec(ind+j)=force(j,i)
+ enddo
+ ind=ind+3
+ enddo
+ do i=nnt,nct
+ write (iout,*) i,ind
+ do j=1,3
+ forcvec(j+ind)=force(j,i+nres)
+ enddo
+ ind=ind+3
+ enddo
+
+ write (iout,*) "forcvec"
+ ind=0
+ do i=nnt,nct-1
+ do j=1,3
+ write (iout,'(2i3,2f10.5)') i,j,force(j,i),&
+ forcvec(ind+j)
+ enddo
+ ind=ind+3
+ enddo
+ do i=nnt,nct
+ do j=1,3
+ write (iout,'(2i3,2f10.5)') i,j,force(j,i+nres),&
+ forcvec(ind+j)
+ enddo
+ ind=ind+3
+ enddo
+
+ endif
+
+ return
+ end subroutine stochastic_force
+!-----------------------------------------------------------------------------
+ subroutine sdarea(gamvec)
+!
+! Scale the friction coefficients according to solvent accessible surface areas
+! Code adapted from TINKER
+! AL 9/3/04
+!
+ use energy_data
+! 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'
+ real(kind=8),dimension(2*nres) :: radius,gamvec !(maxres2)
+ real(kind=8),parameter :: twosix = 1.122462048309372981d0
+ logical :: lprn = .false.
+ real(kind=8) :: probe,area,ratio
+ integer :: i,j,ind,iti
+!
+! determine new friction coefficients every few SD steps
+!
+! set the atomic radii to estimates of sigma values
+!
+! print *,"Entered sdarea"
+ probe = 0.0d0
+
+ do i=1,2*nres
+ radius(i)=0.0d0
+ enddo
+! Load peptide group radii
+ do i=nnt,nct-1
+ radius(i)=pstok
+ enddo
+! Load side chain radii
+ do i=nnt,nct
+ iti=itype(i)
+ radius(i+nres)=restok(iti)
+ enddo
+! do i=1,2*nres
+! write (iout,*) "i",i," radius",radius(i)
+! enddo
+ do i = 1, 2*nres
+ radius(i) = radius(i) / twosix
+ if (radius(i) .ne. 0.0d0) radius(i) = radius(i) + probe
+ end do
+!
+! scale atomic friction coefficients by accessible area
+!
+ if (lprn) write (iout,*) &
+ "Original gammas, surface areas, scaling factors, new gammas, ",&
+ "std's of stochastic forces"
+ ind=0
+ do i=nnt,nct-1
+ if (radius(i).gt.0.0d0) then
+ call surfatom (i,area,radius)
+ ratio = dmax1(area/(4.0d0*pi*radius(i)**2),1.0d-1)
+ if (lprn) write (iout,'(i5,3f10.5,$)') &
+ i,gamvec(ind+1),area,ratio
+ do j=1,3
+ ind=ind+1
+ gamvec(ind) = ratio * gamvec(ind)
+ enddo
+ stdforcp(i)=stdfp*dsqrt(gamvec(ind))
+ if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcp(i)
+ endif
+ enddo
+ do i=nnt,nct
+ if (radius(i+nres).gt.0.0d0) then
+ call surfatom (i+nres,area,radius)
+ ratio = dmax1(area/(4.0d0*pi*radius(i+nres)**2),1.0d-1)
+ if (lprn) write (iout,'(i5,3f10.5,$)') &
+ i,gamvec(ind+1),area,ratio
+ do j=1,3
+ ind=ind+1
+ gamvec(ind) = ratio * gamvec(ind)
+ enddo
+ stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamvec(ind))
+ if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcsc(i)
+ endif
+ enddo
+
+ return
+ end subroutine sdarea
+!-----------------------------------------------------------------------------
+! surfatom.f
+!-----------------------------------------------------------------------------
+!
+!
+! ###################################################
+! ## COPYRIGHT (C) 1996 by Jay William Ponder ##
+! ## All Rights Reserved ##
+! ###################################################
+!
+! ################################################################
+! ## ##
+! ## subroutine surfatom -- exposed surface area of an atom ##
+! ## ##
+! ################################################################
+!
+!
+! "surfatom" performs an analytical computation of the surface
+! area of a specified atom; a simplified version of "surface"
+!
+! literature references:
+!
+! T. J. Richmond, "Solvent Accessible Surface Area and
+! Excluded Volume in Proteins", Journal of Molecular Biology,
+! 178, 63-89 (1984)
+!
+! L. Wesson and D. Eisenberg, "Atomic Solvation Parameters
+! Applied to Molecular Dynamics of Proteins in Solution",
+! Protein Science, 1, 227-235 (1992)
+!
+! variables and parameters:
+!
+! ir number of atom for which area is desired
+! area accessible surface area of the atom
+! radius radii of each of the individual atoms
+!
+!
+ subroutine surfatom(ir,area,radius)
+
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'sizes.i'
+! include 'COMMON.GEO'
+! include 'COMMON.IOUNITS'
+! integer :: nres,
+ integer :: nsup,nstart_sup
+! double precision c,dc,dc_old,d_c_work,xloc,xrot,dc_norm
+! common /chain/ c(3,maxres2+2),dc(3,0:maxres2),dc_old(3,0:maxres2),
+! & xloc(3,maxres),xrot(3,maxres),dc_norm(3,0:maxres2),
+! & dc_work(MAXRES6),nres,nres0
+ integer,parameter :: maxarc=300
+ integer :: i,j,k,m
+ integer :: ii,ib,jb
+ integer :: io,ir
+ integer :: mi,ni,narc
+ integer :: key(maxarc)
+ integer :: intag(maxarc)
+ integer :: intag1(maxarc)
+ real(kind=8) :: area,arcsum
+ real(kind=8) :: arclen,exang
+ real(kind=8) :: delta,delta2
+ real(kind=8) :: eps,rmove
+ real(kind=8) :: xr,yr,zr
+ real(kind=8) :: rr,rrsq
+ real(kind=8) :: rplus,rminus
+ real(kind=8) :: axx,axy,axz
+ real(kind=8) :: ayx,ayy
+ real(kind=8) :: azx,azy,azz
+ real(kind=8) :: uxj,uyj,uzj
+ real(kind=8) :: tx,ty,tz
+ real(kind=8) :: txb,tyb,td
+ real(kind=8) :: tr2,tr,txr,tyr
+ real(kind=8) :: tk1,tk2
+ real(kind=8) :: thec,the,t,tb
+ real(kind=8) :: txk,tyk,tzk
+ real(kind=8) :: t1,ti,tf,tt
+ real(kind=8) :: txj,tyj,tzj
+ real(kind=8) :: ccsq,cc,xysq
+ real(kind=8) :: bsqk,bk,cosine
+ real(kind=8) :: dsqj,gi,pix2
+ real(kind=8) :: therk,dk,gk
+ real(kind=8) :: risqk,rik
+ real(kind=8) :: radius(2*nres) !(maxatm) (maxatm=maxres2)
+ real(kind=8) :: ri(maxarc),risq(maxarc)
+ real(kind=8) :: ux(maxarc),uy(maxarc),uz(maxarc)
+ real(kind=8) :: xc(maxarc),yc(maxarc),zc(maxarc)
+ real(kind=8) :: xc1(maxarc),yc1(maxarc),zc1(maxarc)
+ real(kind=8) :: dsq(maxarc),bsq(maxarc)
+ real(kind=8) :: dsq1(maxarc),bsq1(maxarc)
+ real(kind=8) :: arci(maxarc),arcf(maxarc)
+ real(kind=8) :: ex(maxarc),lt(maxarc),gr(maxarc)
+ real(kind=8) :: b(maxarc),b1(maxarc),bg(maxarc)
+ real(kind=8) :: kent(maxarc),kout(maxarc)
+ real(kind=8) :: ther(maxarc)
+ logical :: moved,top
+ logical :: omit(maxarc)
+!
+! include 'sizes.i'
+ maxatm = 2*nres !maxres2 maxres2=2*maxres
+ maxlight = 8*maxatm
+ maxbnd = 2*maxatm
+ maxang = 3*maxatm
+ maxtors = 4*maxatm
+!
+! zero out the surface area for the sphere of interest
+!
+ area = 0.0d0
+! write (2,*) "ir",ir," radius",radius(ir)
+ if (radius(ir) .eq. 0.0d0) return
+!
+! set the overlap significance and connectivity shift
+!
+ pix2 = 2.0d0 * pi
+ delta = 1.0d-8
+ delta2 = delta * delta
+ eps = 1.0d-8
+ moved = .false.
+ rmove = 1.0d-8
+!
+! store coordinates and radius of the sphere of interest
+!
+ xr = c(1,ir)
+ yr = c(2,ir)
+ zr = c(3,ir)
+ rr = radius(ir)
+ rrsq = rr * rr
+!
+! initialize values of some counters and summations
+!
+ 10 continue
+ io = 0
+ jb = 0
+ ib = 0
+ arclen = 0.0d0
+ exang = 0.0d0
+!
+! test each sphere to see if it overlaps the sphere of interest
+!
+ do i = 1, 2*nres
+ if (i.eq.ir .or. radius(i).eq.0.0d0) goto 30
+ rplus = rr + radius(i)
+ tx = c(1,i) - xr
+ if (abs(tx) .ge. rplus) goto 30
+ ty = c(2,i) - yr
+ if (abs(ty) .ge. rplus) goto 30
+ tz = c(3,i) - zr
+ if (abs(tz) .ge. rplus) goto 30
+!
+! check for sphere overlap by testing distance against radii
+!
+ xysq = tx*tx + ty*ty
+ if (xysq .lt. delta2) then
+ tx = delta
+ ty = 0.0d0
+ xysq = delta2
+ end if
+ ccsq = xysq + tz*tz
+ cc = sqrt(ccsq)
+ if (rplus-cc .le. delta) goto 30
+ rminus = rr - radius(i)
+!
+! check to see if sphere of interest is completely buried
+!
+ if (cc-abs(rminus) .le. delta) then
+ if (rminus .le. 0.0d0) goto 170
+ goto 30
+ end if
+!
+! check for too many overlaps with sphere of interest
+!
+ if (io .ge. maxarc) then
+ write (iout,20)
+ 20 format (/,' SURFATOM -- Increase the Value of MAXARC')
+ stop
+ end if
+!
+! get overlap between current sphere and sphere of interest
+!
+ io = io + 1
+ xc1(io) = tx
+ yc1(io) = ty
+ zc1(io) = tz
+ dsq1(io) = xysq
+ bsq1(io) = ccsq
+ b1(io) = cc
+ gr(io) = (ccsq+rplus*rminus) / (2.0d0*rr*b1(io))
+ intag1(io) = i
+ omit(io) = .false.
+ 30 continue
+ end do
+!
+! case where no other spheres overlap the sphere of interest
+!
+ if (io .eq. 0) then
+ area = 4.0d0 * pi * rrsq
+ return
+ end if
+!
+! case where only one sphere overlaps the sphere of interest
+!
+ if (io .eq. 1) then
+ area = pix2 * (1.0d0 + gr(1))
+ area = mod(area,4.0d0*pi) * rrsq
+ return
+ end if
+!
+! case where many spheres intersect the sphere of interest;
+! sort the intersecting spheres by their degree of overlap
+!
+ call sort2 (io,gr,key)
+ do i = 1, io
+ k = key(i)
+ intag(i) = intag1(k)
+ xc(i) = xc1(k)
+ yc(i) = yc1(k)
+ zc(i) = zc1(k)
+ dsq(i) = dsq1(k)
+ b(i) = b1(k)
+ bsq(i) = bsq1(k)
+ end do
+!
+! get radius of each overlap circle on surface of the sphere
+!
+ do i = 1, io
+ gi = gr(i) * rr
+ bg(i) = b(i) * gi
+ risq(i) = rrsq - gi*gi
+ ri(i) = sqrt(risq(i))
+ ther(i) = 0.5d0*pi - asin(min(1.0d0,max(-1.0d0,gr(i))))
+ end do
+!
+! find boundary of inaccessible area on sphere of interest
+!
+ do k = 1, io-1
+ if (.not. omit(k)) then
+ txk = xc(k)
+ tyk = yc(k)
+ tzk = zc(k)
+ bk = b(k)
+ therk = ther(k)
+!
+! check to see if J circle is intersecting K circle;
+! get distance between circle centers and sum of radii
+!
+ do j = k+1, io
+ if (omit(j)) goto 60
+ cc = (txk*xc(j)+tyk*yc(j)+tzk*zc(j))/(bk*b(j))
+ cc = acos(min(1.0d0,max(-1.0d0,cc)))
+ td = therk + ther(j)
+!
+! check to see if circles enclose separate regions
+!
+ if (cc .ge. td) goto 60
+!
+! check for circle J completely inside circle K
+!
+ if (cc+ther(j) .lt. therk) goto 40
+!
+! check for circles that are essentially parallel
+!
+ if (cc .gt. delta) goto 50
+ 40 continue
+ omit(j) = .true.
+ goto 60
+!
+! check to see if sphere of interest is completely buried
+!
+ 50 continue
+ if (pix2-cc .le. td) goto 170
+ 60 continue
+ end do
+ end if
+ end do
+!
+! find T value of circle intersections
+!
+ do k = 1, io
+ if (omit(k)) goto 110
+ omit(k) = .true.
+ narc = 0
+ top = .false.
+ txk = xc(k)
+ tyk = yc(k)
+ tzk = zc(k)
+ dk = sqrt(dsq(k))
+ bsqk = bsq(k)
+ bk = b(k)
+ gk = gr(k) * rr
+ risqk = risq(k)
+ rik = ri(k)
+ therk = ther(k)
+!
+! rotation matrix elements
+!
+ t1 = tzk / (bk*dk)
+ axx = txk * t1
+ axy = tyk * t1
+ axz = dk / bk
+ ayx = tyk / dk
+ ayy = txk / dk
+ azx = txk / bk
+ azy = tyk / bk
+ azz = tzk / bk
+ do j = 1, io
+ if (.not. omit(j)) then
+ txj = xc(j)
+ tyj = yc(j)
+ tzj = zc(j)
+!
+! rotate spheres so K vector colinear with z-axis
+!
+ uxj = txj*axx + tyj*axy - tzj*axz
+ uyj = tyj*ayy - txj*ayx
+ uzj = txj*azx + tyj*azy + tzj*azz
+ cosine = min(1.0d0,max(-1.0d0,uzj/b(j)))
+ if (acos(cosine) .lt. therk+ther(j)) then
+ dsqj = uxj*uxj + uyj*uyj
+ tb = uzj*gk - bg(j)
+ txb = uxj * tb
+ tyb = uyj * tb
+ td = rik * dsqj
+ tr2 = risqk*dsqj - tb*tb
+ tr2 = max(eps,tr2)
+ tr = sqrt(tr2)
+ txr = uxj * tr
+ tyr = uyj * tr
+!
+! get T values of intersection for K circle
+!
+ tb = (txb+tyr) / td
+ tb = min(1.0d0,max(-1.0d0,tb))
+ tk1 = acos(tb)
+ if (tyb-txr .lt. 0.0d0) tk1 = pix2 - tk1
+ tb = (txb-tyr) / td
+ tb = min(1.0d0,max(-1.0d0,tb))
+ tk2 = acos(tb)
+ if (tyb+txr .lt. 0.0d0) tk2 = pix2 - tk2
+ thec = (rrsq*uzj-gk*bg(j)) / (rik*ri(j)*b(j))
+ if (abs(thec) .lt. 1.0d0) then
+ the = -acos(thec)
+ else if (thec .ge. 1.0d0) then
+ the = 0.0d0
+ else if (thec .le. -1.0d0) then
+ the = -pi
+ end if
+!
+! see if "tk1" is entry or exit point; check t=0 point;
+! "ti" is exit point, "tf" is entry point
+!
+ cosine = min(1.0d0,max(-1.0d0, &
+ (uzj*gk-uxj*rik)/(b(j)*rr)))
+ if ((acos(cosine)-ther(j))*(tk2-tk1) .le. 0.0d0) then
+ ti = tk2
+ tf = tk1
+ else
+ ti = tk2
+ tf = tk1
+ end if
+ narc = narc + 1
+ if (narc .ge. maxarc) then
+ write (iout,70)
+ 70 format (/,' SURFATOM -- Increase the Value',&
+ ' of MAXARC')
+ stop
+ end if
+ if (tf .le. ti) then
+ arcf(narc) = tf
+ arci(narc) = 0.0d0
+ tf = pix2
+ lt(narc) = j
+ ex(narc) = the
+ top = .true.
+ narc = narc + 1
+ end if
+ arcf(narc) = tf
+ arci(narc) = ti
+ lt(narc) = j
+ ex(narc) = the
+ ux(j) = uxj
+ uy(j) = uyj
+ uz(j) = uzj
+ end if
+ end if
+ end do
+ omit(k) = .false.
+!
+! special case; K circle without intersections
+!
+ if (narc .le. 0) goto 90
+!
+! general case; sum up arclength and set connectivity code
+!
+ call sort2 (narc,arci,key)
+ arcsum = arci(1)
+ mi = key(1)
+ t = arcf(mi)
+ ni = mi
+ if (narc .gt. 1) then
+ do j = 2, narc
+ m = key(j)
+ if (t .lt. arci(j)) then
+ arcsum = arcsum + arci(j) - t
+ exang = exang + ex(ni)
+ jb = jb + 1
+ if (jb .ge. maxarc) then
+ write (iout,80)
+ 80 format (/,' SURFATOM -- Increase the Value',&
+ ' of MAXARC')
+ stop
+ end if
+ i = lt(ni)
+ kent(jb) = maxarc*i + k
+ i = lt(m)
+ kout(jb) = maxarc*k + i
+ end if
+ tt = arcf(m)
+ if (tt .ge. t) then
+ t = tt
+ ni = m
+ end if
+ end do
+ end if
+ arcsum = arcsum + pix2 - t
+ if (.not. top) then
+ exang = exang + ex(ni)
+ jb = jb + 1
+ i = lt(ni)
+ kent(jb) = maxarc*i + k
+ i = lt(mi)
+ kout(jb) = maxarc*k + i
+ end if
+ goto 100
+ 90 continue
+ arcsum = pix2
+ ib = ib + 1
+ 100 continue
+ arclen = arclen + gr(k)*arcsum
+ 110 continue
+ end do
+ if (arclen .eq. 0.0d0) goto 170
+ if (jb .eq. 0) goto 150
+!
+! find number of independent boundaries and check connectivity
+!
+ j = 0
+ do k = 1, jb
+ if (kout(k) .ne. 0) then
+ i = k
+ 120 continue
+ m = kout(i)
+ kout(i) = 0
+ j = j + 1
+ do ii = 1, jb
+ if (m .eq. kent(ii)) then
+ if (ii .eq. k) then
+ ib = ib + 1
+ if (j .eq. jb) goto 150
+ goto 130
+ end if
+ i = ii
+ goto 120
+ end if
+ end do
+ 130 continue
+ end if
+ end do
+ ib = ib + 1
+!
+! attempt to fix connectivity error by moving atom slightly
+!
+ if (moved) then
+ write (iout,140) ir
+ 140 format (/,' SURFATOM -- Connectivity Error at Atom',i6)
+ else
+ moved = .true.
+ xr = xr + rmove
+ yr = yr + rmove
+ zr = zr + rmove
+ goto 10
+ end if
+!
+! compute the exposed surface area for the sphere of interest
+!
+ 150 continue
+ area = ib*pix2 + exang + arclen
+ area = mod(area,4.0d0*pi) * rrsq
+!
+! attempt to fix negative area by moving atom slightly
+!
+ if (area .lt. 0.0d0) then
+ if (moved) then
+ write (iout,160) ir
+ 160 format (/,' SURFATOM -- Negative Area at Atom',i6)
+ else
+ moved = .true.
+ xr = xr + rmove
+ yr = yr + rmove
+ zr = zr + rmove
+ goto 10
+ end if
+ end if
+ 170 continue
+ return
+ end subroutine surfatom
+!----------------------------------------------------------------
+!----------------------------------------------------------------
+ subroutine alloc_MD_arrays
+!EL Allocation of arrays used by MD module
+
+ integer :: nres2,nres6
+ nres2=nres*2
+ nres6=nres*6
+!----------------------
+#ifndef LANG0
+! commom.langevin
+! common /langforc/
+ allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
+ allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
+ if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
+ allocate(fricvec(nres2,nres2))
+ allocate(pfric_mat(nres2,nres2),vfric_mat(nres2,nres2))
+ allocate(afric_mat(nres2,nres2),prand_mat(nres2,nres2))
+ allocate(vrand_mat1(nres2,nres2),vrand_mat2(nres2,nres2)) !(MAXRES2,MAXRES2)
+ allocate(pfric0_mat(nres2,nres2,0:maxflag_stoch))
+ allocate(afric0_mat(nres2,nres2,0:maxflag_stoch))
+ allocate(vfric0_mat(nres2,nres2,0:maxflag_stoch))
+ allocate(prand0_mat(nres2,nres2,0:maxflag_stoch))
+ allocate(vrand0_mat1(nres2,nres2,0:maxflag_stoch))
+ allocate(vrand0_mat2(nres2,nres2,0:maxflag_stoch)) !(MAXRES2,MAXRES2,0:maxflag_stoch)
+ allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
+! common /langmat/
+ allocate(mt1(nres2,nres2),mt2(nres2,nres2),mt3(nres2,nres2)) !(maxres2,maxres2)
+!----------------------
+#else
+! commom.langevin.lang0
+! common /langforc/
+ allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
+ if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
+ allocate(fricvec(nres2,nres2)) !(MAXRES2,MAXRES2)
+ allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
+ allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
+#endif
+
+!el if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2))
+!----------------------
+! commom.hairpin in CSA module
+!----------------------
+! common.mce in MCM_MD module
+!----------------------
+! common.MD
+! common /mdgrad/ in module.energy
+! common /back_constr/ in module.energy
+! common /qmeas/ in module.energy
+! common /mdpar/
+! common /MDcalc/
+ allocate(potEcomp(0:n_ene+4)) !(0:n_ene+4)
+! common /lagrange/
+ allocate(d_t(3,0:nres2),d_a(3,0:nres2),d_t_old(3,0:nres2)) !(3,0:MAXRES2)
+ allocate(d_a_work(nres6)) !(6*MAXRES)
+ allocate(Gmat(nres2,nres2),A(nres2,nres2))
+ if(.not.allocated(Ginv)) allocate(Ginv(nres2,nres2)) !in control: ergastulum
+ allocate(Gsqrp(nres2,nres2),Gsqrm(nres2,nres2),Gvec(nres2,nres2)) !(maxres2,maxres2)
+ allocate(Geigen(nres2)) !(maxres2)
+ if(.not.allocated(vtot)) allocate(vtot(nres2)) !(maxres2)
+! common /inertia/ in io_conf: parmread
+! real(kind=8),dimension(:),allocatable :: ISC,msc !(ntyp+1)
+! common /langevin/in io read_MDpar
+! real(kind=8),dimension(:),allocatable :: gamsc !(ntyp1)
+! real(kind=8),dimension(:),allocatable :: stdfsc !(ntyp)
+! in io_conf: parmread
+! real(kind=8),dimension(:),allocatable :: restok !(ntyp+1)
+! common /mdpmpi/ in control: ergastulum
+ if(.not.allocated(ng_start)) allocate(ng_start(0:nfgtasks-1))
+ if(.not.allocated(ng_counts)) allocate(ng_counts(0:nfgtasks-1))
+ if(.not.allocated(nginv_counts)) allocate(nginv_counts(0:nfgtasks-1)) !(0:MaxProcs-1)
+ if(.not.allocated(nginv_start)) allocate(nginv_start(0:nfgtasks)) !(0:MaxProcs)
+!----------------------
+! common.muca in read_muca
+! common /double_muca/
+! real(kind=8) :: elow,ehigh,factor,hbin,factor_min
+! real(kind=8),dimension(:),allocatable :: emuca,nemuca,&
+! nemuca2,hist !(4*maxres)
+! common /integer_muca/
+! integer :: nmuca,imtime,muca_smooth
+! common /mucarem/
+! real(kind=8),dimension(:),allocatable :: elowi,ehighi !(maxprocs)
+!----------------------
+! common.MD
+! common /mdgrad/ in module.energy
+! common /back_constr/ in module.energy
+! common /qmeas/ in module.energy
+! common /mdpar/
+! common /MDcalc/
+! common /lagrange/
+ allocate(d_t_work(nres6),d_t_work_new(nres6),d_af_work(nres6))
+ allocate(d_as_work(nres6),kinetic_force(nres6)) !(MAXRES6)
+ allocate(d_t_new(3,0:nres2),d_a_old(3,0:nres2),d_a_short(3,0:nres2)) !,d_a !(3,0:MAXRES2)
+ allocate(stdforcp(nres),stdforcsc(nres)) !(MAXRES)
+!----------------------
+! COMMON /BANII/ D
+ allocate(D_ban(nres6)) !(MAXRES6) maxres6=6*maxres
+! common /stochcalc/ stochforcvec
+ allocate(stochforcvec(nres6)) !(MAXRES6) maxres6=6*maxres
+!----------------------
+ return
+ end subroutine alloc_MD_arrays
+!-----------------------------------------------------------------------------
+!-----------------------------------------------------------------------------
+ end module MDyn
projects / unres4.git / blobdiff