+++ /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