2 !-----------------------------------------------------------------------------
15 !-----------------------------------------------------------------------------
17 ! common /mdgrad/ in module.energy
18 ! common /back_constr/ in module.energy
19 ! common /qmeas/ in module.energy
23 real(kind=8),dimension(:),allocatable :: d_t_work,&
24 d_t_work_new,d_af_work,d_as_work,kinetic_force !(MAXRES6)
25 real(kind=8),dimension(:,:),allocatable :: d_t_new,&
26 d_a_old,d_a_short!,d_a !(3,0:MAXRES2)
27 ! real(kind=8),dimension(:),allocatable :: d_a_work !(6*MAXRES)
28 ! real(kind=8),dimension(:,:),allocatable :: Gmat,Ginv,A,&
29 ! Gsqrp,Gsqrm,Gvec !(maxres2,maxres2)
30 ! real(kind=8),dimension(:),allocatable :: Geigen !(maxres2)
31 ! integer :: dimen,dimen1,dimen3
32 ! integer :: lang,count_reset_moment,count_reset_vel
33 ! logical :: reset_moment,reset_vel,rattle,RESPA
36 ! real(kind=8) :: rwat,etawat,stdfp,cPoise
37 ! real(kind=8),dimension(:),allocatable :: gamsc !(ntyp1)
38 ! real(kind=8),dimension(:),allocatable :: stdfsc !(ntyp)
39 real(kind=8),dimension(:),allocatable :: stdforcp,stdforcsc !(MAXRES)
40 !-----------------------------------------------------------------------------
44 ! ###################################################
45 ! ## COPYRIGHT (C) 1992 by Jay William Ponder ##
46 ! ## All Rights Reserved ##
47 ! ###################################################
49 ! #############################################################
51 ! ## sizes.i -- parameter values to set array dimensions ##
53 ! #############################################################
56 ! "sizes.i" sets values for critical array dimensions used
57 ! throughout the software; these parameters will fix the size
58 ! of the largest systems that can be handled; values too large
59 ! for the computer's memory and/or swap space to accomodate
60 ! will result in poor performance or outright failure
62 ! parameter: maximum allowed number of:
64 ! maxatm atoms in the molecular system
65 ! maxval atoms directly bonded to an atom
66 ! maxgrp ! user-defined groups of atoms
67 ! maxtyp force field atom type definitions
68 ! maxclass force field atom class definitions
69 ! maxkey lines in the keyword file
70 ! maxrot bonds for torsional rotation
71 ! maxvar optimization variables (vector storage)
72 ! maxopt optimization variables (matrix storage)
73 ! maxhess off-diagonal Hessian elements
74 ! maxlight sites for method of lights neighbors
75 ! maxvib vibrational frequencies
76 ! maxgeo distance geometry points
77 ! maxcell unit cells in replicated crystal
78 ! maxring 3-, 4-, or 5-membered rings
79 ! maxfix geometric restraints
80 ! maxbio biopolymer atom definitions
81 ! maxres residues in the macromolecule
82 ! maxamino amino acid residue types
83 ! maxnuc nucleic acid residue types
84 ! maxbnd covalent bonds in molecular system
85 ! maxang bond angles in molecular system
86 ! maxtors torsional angles in molecular system
87 ! maxpi atoms in conjugated pisystem
88 ! maxpib covalent bonds involving pisystem
89 ! maxpit torsional angles involving pisystem
92 !el integer maxatm,maxval,maxgrp
93 !el integer maxtyp,maxclass,maxkey
94 !el integer maxrot,maxopt
95 !el integer maxhess,maxlight,maxvib
96 !el integer maxgeo,maxcell,maxring
97 !el integer maxfix,maxbio
98 !el integer maxamino,maxnuc,maxbnd
99 !el integer maxang,maxtors,maxpi
100 !el integer maxpib,maxpit
101 ! integer :: maxatm !=2*nres !maxres2 maxres2=2*maxres
102 ! integer,parameter :: maxval=8
103 ! integer,parameter :: maxgrp=1000
104 ! integer,parameter :: maxtyp=3000
105 ! integer,parameter :: maxclass=500
106 ! integer,parameter :: maxkey=10000
107 ! integer,parameter :: maxrot=1000
108 ! integer,parameter :: maxopt=1000
109 ! integer,parameter :: maxhess=1000000
110 ! integer :: maxlight !=8*maxatm
111 ! integer,parameter :: maxvib=1000
112 ! integer,parameter :: maxgeo=1000
113 ! integer,parameter :: maxcell=10000
114 ! integer,parameter :: maxring=10000
115 ! integer,parameter :: maxfix=10000
116 ! integer,parameter :: maxbio=10000
117 ! integer,parameter :: maxamino=31
118 ! integer,parameter :: maxnuc=12
119 ! integer :: maxbnd !=2*maxatm
120 ! integer :: maxang !=3*maxatm
121 ! integer :: maxtors !=4*maxatm
122 ! integer,parameter :: maxpi=100
123 ! integer,parameter :: maxpib=2*maxpi
124 ! integer,parameter :: maxpit=4*maxpi
125 !-----------------------------------------------------------------------------
126 ! Maximum number of seed
127 ! integer,parameter :: max_seed=1
128 !-----------------------------------------------------------------------------
129 real(kind=8),dimension(:),allocatable :: stochforcvec !(MAXRES6) maxres6=6*maxres
130 ! common /stochcalc/ stochforcvec
131 !-----------------------------------------------------------------------------
132 ! common /przechowalnia/ subroutines: rattle1,rattle2,rattle_brown
133 real(kind=8),dimension(:,:),allocatable :: GGinv !(2*nres,2*nres) maxres2=2*maxres
134 real(kind=8),dimension(:,:,:),allocatable :: gdc !(3,2*nres,2*nres) maxres2=2*maxres
135 real(kind=8),dimension(:,:),allocatable :: Cmat !(2*nres,2*nres) maxres2=2*maxres
136 !-----------------------------------------------------------------------------
137 ! common /syfek/ subroutines: friction_force,setup_fricmat
138 !el real(kind=8),dimension(:),allocatable :: gamvec !(MAXRES6) or (MAXRES2)
139 !-----------------------------------------------------------------------------
140 ! common /przechowalnia/ subroutines: friction_force,setup_fricmat
141 real(kind=8),dimension(:,:),allocatable :: ginvfric !(2*nres,2*nres) !maxres2=2*maxres
142 !-----------------------------------------------------------------------------
143 ! common /przechowalnia/ subroutine: setup_fricmat
144 real(kind=8),dimension(:,:),allocatable :: fcopy !(2*nres,2*nres)
145 !-----------------------------------------------------------------------------
148 !-----------------------------------------------------------------------------
150 !-----------------------------------------------------------------------------
152 !-----------------------------------------------------------------------------
153 subroutine brown_step(itime)
154 !------------------------------------------------
155 ! Perform a single Euler integration step of Brownian dynamics
156 !------------------------------------------------
157 ! implicit real*8 (a-h,o-z)
159 use control, only: tcpu
162 ! use io_conf, only:cartprint
163 ! include 'DIMENSIONS'
167 ! include 'COMMON.CONTROL'
168 ! include 'COMMON.VAR'
169 ! include 'COMMON.MD'
171 ! include 'COMMON.LANGEVIN'
173 ! include 'COMMON.LANGEVIN.lang0'
175 ! include 'COMMON.CHAIN'
176 ! include 'COMMON.DERIV'
177 ! include 'COMMON.GEO'
178 ! include 'COMMON.LOCAL'
179 ! include 'COMMON.INTERACT'
180 ! include 'COMMON.IOUNITS'
181 ! include 'COMMON.NAMES'
182 ! include 'COMMON.TIME1'
183 real(kind=8),dimension(6*nres) :: zapas !(MAXRES6) maxres6=6*maxres
184 integer :: rstcount !ilen,
186 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
187 real(kind=8),dimension(6*nres,2*nres) :: Bmat,GBmat,Tmat !(MAXRES6,MAXRES2) (maxres2=2*maxres,maxres6=6*maxres)
188 real(kind=8),dimension(2*nres,2*nres) :: Cmat_,Cinv !(maxres2,maxres2) maxres2=2*maxres
189 real(kind=8),dimension(6*nres,6*nres) :: Pmat !(maxres6,maxres6) maxres6=6*maxres
190 ! real(kind=8),dimension(:,:),allocatable :: Bmat,GBmat,Tmat !(MAXRES6,MAXRES2) (maxres2=2*maxres,maxres6=6*maxres)
191 ! real(kind=8),dimension(:,:),allocatable :: Cmat_,Cinv !(maxres2,maxres2) maxres2=2*maxres
192 ! real(kind=8),dimension(:,:),allocatable :: Pmat !(maxres6,maxres6) maxres6=6*maxres
193 real(kind=8),dimension(6*nres) :: Td !(maxres6) maxres6=6*maxres
194 real(kind=8),dimension(2*nres) :: ppvec !(maxres2) maxres2=2*maxres
195 !el common /stochcalc/ stochforcvec
196 !el real(kind=8),dimension(3) :: cm !el
197 !el common /gucio/ cm
199 logical :: lprn = .false.,lprn1 = .false.
200 integer :: maxiter = 5
201 real(kind=8) :: difftol = 1.0d-5
202 real(kind=8) :: xx,diffmax,blen2,diffbond,tt0
203 integer :: i,j,nbond,k,ind,ind1,iter
204 integer :: nres2,nres6
208 ! if (.not.allocated(Bmat)) allocate(Bmat(nres6,nres2))
209 ! if (.not.allocated(GBmat)) allocate (GBmat(nres6,nres2))
210 ! if (.not.allocated(Tmat)) allocate (Tmat(nres6,nres2))
211 ! if (.not.allocated(Cmat_)) allocate(Cmat_(nres2,nres2))
212 ! if (.not.allocated(Cinv)) allocate (Cinv(nres2,nres2))
213 ! if (.not.allocated(Pmat)) allocate(Pmat(6*nres,6*nres))
215 if (.not.allocated(stochforcvec)) allocate(stochforcvec(nres6)) !(MAXRES6) maxres6=6*maxres
219 if (itype(i,1).ne.10) nbond=nbond+1
223 write (iout,*) "Generalized inverse of fricmat"
224 call matout(dimen,dimen,nres6,nres6,fricmat)
236 Bmat(ind+j,ind1)=dC_norm(j,i)
241 if (itype(i,1).ne.10) then
244 Bmat(ind+j,ind1)=dC_norm(j,i+nres)
250 write (iout,*) "Matrix Bmat"
251 call MATOUT(nbond,dimen,nres6,nres6,Bmat)
257 GBmat(i,j)=GBmat(i,j)+fricmat(i,k)*Bmat(k,j)
262 write (iout,*) "Matrix GBmat"
263 call MATOUT(nbond,dimen,nres6,nres2,Gbmat)
269 Cmat_(i,j)=Cmat_(i,j)+Bmat(k,i)*GBmat(k,j)
274 write (iout,*) "Matrix Cmat"
275 call MATOUT(nbond,nbond,nres2,nres2,Cmat_)
277 call matinvert(nbond,nres2,Cmat_,Cinv,osob)
279 write (iout,*) "Matrix Cinv"
280 call MATOUT(nbond,nbond,nres2,nres2,Cinv)
286 Tmat(i,j)=Tmat(i,j)+GBmat(i,k)*Cinv(k,j)
291 write (iout,*) "Matrix Tmat"
292 call MATOUT(nbond,dimen,nres6,nres2,Tmat)
302 Pmat(i,j)=Pmat(i,j)-Tmat(i,k)*Bmat(j,k)
307 write (iout,*) "Matrix Pmat"
308 call MATOUT(dimen,dimen,nres6,nres6,Pmat)
315 Td(i)=Td(i)+vbl*Tmat(i,ind)
318 if (itype(k,1).ne.10) then
320 Td(i)=Td(i)+vbldsc0(1,itype(k,1))*Tmat(i,ind)
325 write (iout,*) "Vector Td"
327 write (iout,'(i5,f10.5)') i,Td(i)
330 call stochastic_force(stochforcvec)
332 write (iout,*) "stochforcvec"
334 write (iout,*) i,stochforcvec(i)
338 zapas(j)=-gcart(j,0)+stochforcvec(j)
340 dC_work(j)=dC_old(j,0)
346 zapas(ind)=-gcart(j,i)+stochforcvec(ind)
347 dC_work(ind)=dC_old(j,i)
351 if (itype(i,1).ne.10) then
354 zapas(ind)=-gxcart(j,i)+stochforcvec(ind)
355 dC_work(ind)=dC_old(j,i+nres)
361 write (iout,*) "Initial d_t_work"
363 write (iout,*) i,d_t_work(i)
370 d_t_work(i)=d_t_work(i)+fricmat(i,j)*zapas(j)
377 zapas(i)=zapas(i)+Pmat(i,j)*(dC_work(j)+d_t_work(j)*d_time)
381 write (iout,*) "Final d_t_work and zapas"
383 write (iout,*) i,d_t_work(i),zapas(i)
397 dc_work(ind+j)=dc(j,i)
403 d_t(j,i+nres)=d_t_work(ind+j)
404 dc(j,i+nres)=zapas(ind+j)
405 dc_work(ind+j)=dc(j,i+nres)
411 write (iout,*) "Before correction for rotational lengthening"
412 write (iout,*) "New coordinates",&
413 " and differences between actual and standard bond lengths"
418 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
422 if (itype(i,1).ne.10) then
424 xx=vbld(i+nres)-vbldsc0(1,itype(i,1))
425 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
426 i,(dC(j,i+nres),j=1,3),xx
430 ! Second correction (rotational lengthening)
436 blen2 = scalar(dc(1,i),dc(1,i))
437 ppvec(ind)=2*vbl**2-blen2
438 diffbond=dabs(vbl-dsqrt(blen2))
439 if (diffbond.gt.diffmax) diffmax=diffbond
440 if (ppvec(ind).gt.0.0d0) then
441 ppvec(ind)=dsqrt(ppvec(ind))
446 write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind)
450 if (itype(i,1).ne.10) then
452 blen2 = scalar(dc(1,i+nres),dc(1,i+nres))
453 ppvec(ind)=2*vbldsc0(1,itype(i,1))**2-blen2
454 diffbond=dabs(vbldsc0(1,itype(i,1))-dsqrt(blen2))
455 if (diffbond.gt.diffmax) diffmax=diffbond
456 if (ppvec(ind).gt.0.0d0) then
457 ppvec(ind)=dsqrt(ppvec(ind))
462 write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind)
466 if (lprn) write (iout,*) "iter",iter," diffmax",diffmax
467 if (diffmax.lt.difftol) goto 10
471 Td(i)=Td(i)+ppvec(j)*Tmat(i,j)
477 zapas(i)=zapas(i)+Pmat(i,j)*dc_work(j)
488 dc_work(ind+j)=zapas(ind+j)
493 if (itype(i,1).ne.10) then
495 dc(j,i+nres)=zapas(ind+j)
496 dc_work(ind+j)=zapas(ind+j)
501 ! Building the chain from the newly calculated coordinates
504 if (large.and. mod(itime,ntwe).eq.0) then
505 write (iout,*) "Cartesian and internal coordinates: step 1"
508 write (iout,'(a)') "Potential forces"
510 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(-gcart(j,i),j=1,3),&
513 write (iout,'(a)') "Stochastic forces"
515 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(stochforc(j,i),j=1,3),&
516 (stochforc(j,i+nres),j=1,3)
518 write (iout,'(a)') "Velocities"
520 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
521 (d_t(j,i+nres),j=1,3)
526 write (iout,*) "After correction for rotational lengthening"
527 write (iout,*) "New coordinates",&
528 " and differences between actual and standard bond lengths"
533 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
537 if (itype(i,1).ne.10) then
539 xx=vbld(i+nres)-vbldsc0(1,itype(i,1))
540 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
541 i,(dC(j,i+nres),j=1,3),xx
546 ! write (iout,*) "Too many attempts at correcting the bonds"
554 ! Calculate energy and forces
556 call etotal(potEcomp)
557 potE=potEcomp(0)-potEcomp(51)
561 ! Calculate the kinetic and total energy and the kinetic temperature
564 t_enegrad=t_enegrad+MPI_Wtime()-tt0
566 t_enegrad=t_enegrad+tcpu()-tt0
569 kinetic_T=2.0d0/(dimen*Rb)*EK
571 end subroutine brown_step
572 !-----------------------------------------------------------------------------
574 !-----------------------------------------------------------------------------
575 subroutine gauss(RO,AP,MT,M,N,*)
577 ! CALCULATES (RO**(-1))*AP BY GAUSS ELIMINATION
578 ! RO IS A SQUARE MATRIX
579 ! THE CALCULATED PRODUCT IS STORED IN AP
580 ! ABNORMAL EXIT IF RO IS SINGULAR
582 integer :: MT, M, N, M1,I,J,IM,&
584 real(kind=8) :: RO(MT,M),AP(MT,N),X,RM,PR,Y
590 if(dabs(X).le.1.0D-13) return 1
602 if(DABS(RO(J,I)).LE.RM) goto 2
616 if(dabs(X).le.1.0E-13) return 1
625 8 AP(J,K)=AP(J,K)-Y*AP(I,K)
627 9 RO(J,K)=RO(J,K)-Y*RO(I,K)
631 if(dabs(X).le.1.0E-13) return 1
641 15 X=X-AP(K,J)*RO(MI,K)
646 !-----------------------------------------------------------------------------
649 subroutine kinetic(KE_total)
650 !c----------------------------------------------------------------
651 !c This subroutine calculates the total kinetic energy of the chain
652 !c-----------------------------------------------------------------
653 !c 3/5/2020 AL Corrected for multichain systems, no fake peptide groups
654 !c inside, implemented with five-diagonal inertia matrix
657 real(kind=8):: KE_total,KEt_p,KEt_sc,KEr_p,KEr_sc,mag1,mag2
658 integer i,j,k,iti,mnum
659 real(kind=8),dimension(3) :: incr,v
665 !c write (iout,*) "ISC",(isc(itype(i)),i=1,nres)
666 !c The translational part for peptide virtual bonds
670 do i=nnt,nct-1 !czy na pewno nct-1??
672 !c write (iout,*) "Kinetic trp:",i,(incr(j),j=1,3
673 !c Skip dummy peptide groups
674 if (itype(i,mnum).ne.ntyp1_molec(mnum)&
675 .and. itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
677 v(j)=incr(j)+0.5d0*d_t(j,i)
679 if (mnum.eq.5) mp(mnum)=0.0d0
680 ! if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
681 !c write (iout,*) "Kinetic trp:",i,(v(j),j=1,3)
682 vtot(i)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
683 KEt_p=KEt_p+mp(mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
686 incr(j)=incr(j)+d_t(j,i)
689 !c write(iout,*) 'KEt_p', KEt_p
690 !c The translational part for the side chain virtual bond
691 !c Only now we can initialize incr with zeros. It must be equal
692 !c to the velocities of the first Calpha.
698 iti=iabs(itype(i,mnum))
699 if (mnum.eq.5) iti=itype(i,mnum)
700 if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
707 v(j)=incr(j)+d_t(j,nres+i)
710 ! if (mnum.ne.5) then
711 ! write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3)
712 ! write (iout,*) "i",i," msc",msc(iti,mnum)," v",(v(j),j=1,3)
713 KEt_sc=KEt_sc+msc(iti,mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
714 vtot(i+nres)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
717 incr(j)=incr(j)+d_t(j,i)
721 ! write(iout,*) 'KEt_sc', KEt_sc
722 ! The part due to stretching and rotation of the peptide groups
725 if (itype(i,mnum).ne.ntyp1_molec(mnum)&
726 .and.itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
727 if (mnum.eq.5) Ip(mnum)=Isc(itype(i,mnum),mnum)
728 ! write (iout,*) "i",i
729 ! write (iout,*) "i",i," mag1",mag1," mag2",mag2
733 !c write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3)
734 KEr_p=KEr_p+Ip(mnum)*(incr(1)*incr(1)+incr(2)*incr(2) &
739 !c write(iout,*) 'KEr_p', KEr_p
740 !c The rotational part of the side chain virtual bond
743 iti=iabs(itype(i,mnum))
744 if (itype(i,1).ne.10.and.itype(i,mnum).ne.ntyp1_molec(mnum)&
747 incr(j)=d_t(j,nres+i)
749 ! write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3)
750 KEr_sc=KEr_sc+Isc(iti,mnum)*(incr(1)*incr(1)+incr(2)*incr(2)+&
754 !c The total kinetic energy
756 ! write(iout,*) ' KEt_p',KEt_p,' KEt_sc',KEt_sc,' KEr_p',KEr_p, &
758 KE_total=0.5d0*(KEt_p+KEt_sc+0.25d0*KEr_p+KEr_sc)
759 !c write (iout,*) "KE_total",KE_total
761 end subroutine kinetic
764 !-----------------------------------------------------------------------------
765 subroutine kinetic(KE_total)
766 !----------------------------------------------------------------
767 ! This subroutine calculates the total kinetic energy of the chain
768 !-----------------------------------------------------------------
770 ! implicit real*8 (a-h,o-z)
771 ! include 'DIMENSIONS'
772 ! include 'COMMON.VAR'
773 ! include 'COMMON.CHAIN'
774 ! include 'COMMON.DERIV'
775 ! include 'COMMON.GEO'
776 ! include 'COMMON.LOCAL'
777 ! include 'COMMON.INTERACT'
778 ! include 'COMMON.MD'
779 ! include 'COMMON.IOUNITS'
780 real(kind=8) :: KE_total,mscab
782 integer :: i,j,k,iti,mnum,term
783 real(kind=8) :: KEt_p,KEt_sc,KEr_p,KEr_sc,incr(3),&
786 write (iout,*) "Velocities, kietic"
788 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
789 (d_t(j,i+nres),j=1,3)
794 ! write (iout,*) "ISC",(isc(itype(i,1)),i=1,nres)
795 ! The translational part for peptide virtual bonds
800 ! if (molnum(nct).gt.3) term=nct
803 if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
804 ! write (iout,*) "Kinetic trp:",i,(incr(j),j=1,3),mp(mnum)
807 v(j)=incr(j)+0.5d0*d_t(j,i)
811 v(j)=incr(j)+0.5d0*d_t(j,i)
814 vtot(i)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
815 KEt_p=KEt_p+mp(mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
817 incr(j)=incr(j)+d_t(j,i)
820 ! write(iout,*) 'KEt_p', KEt_p
821 ! The translational part for the side chain virtual bond
822 ! Only now we can initialize incr with zeros. It must be equal
823 ! to the velocities of the first Calpha.
829 iti=iabs(itype(i,mnum))
830 ! if (mnum.ge.4) then
835 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
836 if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
837 .or.(mnum.ge.4)) then
843 v(j)=incr(j)+d_t(j,nres+i)
846 ! write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3)
847 ! write (iout,*) "i",i," msc",msc(iti,mnum)," v",(v(j),j=1,3)
848 KEt_sc=KEt_sc+mscab*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
849 vtot(i+nres)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
851 incr(j)=incr(j)+d_t(j,i)
855 ! write(iout,*) 'KEt_sc', KEt_sc
856 ! The part due to stretching and rotation of the peptide groups
860 ! write (iout,*) "i",i
861 ! write (iout,*) "i",i," mag1",mag1," mag2",mag2
865 ! write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3)
866 KEr_p=KEr_p+Ip(mnum)*(incr(1)*incr(1)+incr(2)*incr(2) &
870 ! write(iout,*) 'KEr_p', KEr_p
871 ! The rotational part of the side chain virtual bond
875 iti=iabs(itype(i,mnum))
876 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
877 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
878 .and.(mnum.lt.4)) then
880 incr(j)=d_t(j,nres+i)
882 ! write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3)
883 KEr_sc=KEr_sc+Isc(iti,mnum)*(incr(1)*incr(1)+incr(2)*incr(2)+ &
887 ! The total kinetic energy
889 ! write(iout,*) 'KEr_sc', KEr_sc
890 KE_total=0.5d0*(KEt_p+KEt_sc+0.25d0*KEr_p+KEr_sc)
891 ! write (iout,*) "KE_total",KE_total
893 end subroutine kinetic
894 !-----------------------------------------------------------------------------
896 subroutine kinetic_CASC(KE_total)
897 !c----------------------------------------------------------------
898 !c Compute the kinetic energy of the system using the Calpha-SC
900 !c-----------------------------------------------------------------
902 double precision KE_total
904 integer i,j,k,iti,ichain,innt,inct,mnum
905 double precision KEt_p,KEt_sc,KEr_p,KEr_sc,incr(3),&
912 !c write (iout,*) "ISC",(isc(itype(i)),i=1,nres)
913 !c The translational part for peptide virtual bonds
916 innt=chain_border(1,ichain)
917 inct=chain_border(2,ichain)
918 !c write (iout,*) "Kinetic_CASC chain",ichain," innt",innt,
923 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
924 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
925 !c write (iout,*) i,(d_t(j,i),j=1,3),(d_t(j,i+1),j=1,3)
927 v(j)=0.5d0*(d_t(j,i)+d_t(j,i+1))
929 !c write (iout,*) "Kinetic trp i",i," v",(v(j),j=1,3)
930 KEt_p=KEt_p+mp(mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
932 !c write(iout,*) 'KEt_p', KEt_p
933 !c The translational part for the side chain virtual bond
934 !c Only now we can initialize incr with zeros. It must be equal
935 !c to the velocities of the first Calpha.
938 iti=iabs(itype(i,mnum))
939 if (itype(i,1).eq.10.or.mnum.ge.3.or. itype(i,mnum).eq.ntyp1_molec(mnum)) then
940 !c write (iout,*) i,iti,(d_t(j,i),j=1,3)
945 !c write (iout,*) i,iti,(d_t(j,nres+i),j=1,3)
950 !c write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3)
951 !c write (iout,*) "i",i," msc",msc(iti)," v",(v(j),j=1,3)
952 KEt_sc=KEt_sc+msc(iti,mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
955 !c write(iout,*) 'KEt_sc', KEt_sc
956 !c The part due to stretching and rotation of the peptide groups
960 incr(j)=d_t(j,i+1)-d_t(j,i)
962 if (mnum.eq.5) Ip(mnum)=Isc(itype(i,mnum),mnum)
963 !c write (iout,*) i,(incr(j),j=1,3)
964 !c write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3)
965 KEr_p=KEr_p+Ip(mnum)*(incr(1)*incr(1)+incr(2)*incr(2)&
969 !c write(iout,*) 'KEr_p', KEr_p
970 !c The rotational part of the side chain virtual bond
973 iti=iabs(itype(i,mnum))
974 ! if (iti.ne.10.and.mnum.lt.3) then
975 if (itype(i,1).ne.10.and.mnum.lt.3.and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
977 incr(j)=d_t(j,nres+i)-d_t(j,i)
979 !c write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3)
980 KEr_sc=KEr_sc+Isc(iti,mnum)*(incr(1)*incr(1)+incr(2)*incr(2)+&
986 !c The total kinetic energy
988 !c write(iout,*) ' KEt_p',KEt_p,' KEt_sc',KEt_sc,' KEr_p',KEr_p,
989 !c & ' KEr_sc', KEr_sc
990 KE_total=0.5d0*(KEt_p+KEt_sc+0.25d0*KEr_p+KEr_sc)
991 !c write (iout,*) "KE_total",KE_tota
993 write (iout,*) "Need to compile with -DFIVEDIAG to use this sub!"
997 end subroutine kinetic_CASC
1000 !-----------------------------------------------------------------------------
1002 !------------------------------------------------
1003 ! The driver for molecular dynamics subroutines
1004 !------------------------------------------------
1007 use control, only:tcpu,ovrtim
1008 ! use io_comm, only:ilen
1010 use compare, only:secondary2,hairpin
1011 use io, only:cartout,statout
1012 ! implicit real*8 (a-h,o-z)
1013 ! include 'DIMENSIONS'
1016 integer :: IERROR,ERRCODE
1018 ! include 'COMMON.SETUP'
1019 ! include 'COMMON.CONTROL'
1020 ! include 'COMMON.VAR'
1021 ! include 'COMMON.MD'
1023 ! include 'COMMON.LANGEVIN'
1025 ! include 'COMMON.LANGEVIN.lang0'
1027 ! include 'COMMON.CHAIN'
1028 ! include 'COMMON.DERIV'
1029 ! include 'COMMON.GEO'
1030 ! include 'COMMON.LOCAL'
1031 ! include 'COMMON.INTERACT'
1032 ! include 'COMMON.IOUNITS'
1033 ! include 'COMMON.NAMES'
1034 ! include 'COMMON.TIME1'
1035 ! include 'COMMON.HAIRPIN'
1036 real(kind=8),dimension(3) :: L,vcm
1038 real(kind=8),dimension(6*nres) :: v_work,v_transf !(maxres6) maxres6=6*maxres
1040 integer :: rstcount !ilen,
1042 character(len=50) :: tytul
1043 !el common /gucio/ cm
1044 integer :: i,j,nharp
1045 integer,dimension(4,nres) :: iharp !(4,nres/3)(4,maxres/3)
1047 real(kind=8) :: tt0,scalfac
1048 integer :: nres2,itime
1053 print *,"MY tmpdir",tmpdir,ilen(tmpdir)
1054 if (ilen(tmpdir).gt.0) &
1055 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_" &
1056 //liczba(:ilen(liczba))//'.rst')
1058 if (ilen(tmpdir).gt.0) &
1059 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"//'.rst')
1066 write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
1072 print *,"just befor setup matix",nres
1073 ! Determine the inverse of the inertia matrix.
1074 call setup_MD_matrices
1076 print *,"AFTER SETUP MATRICES"
1078 print *,"AFTER INIT MD"
1081 t_MDsetup = MPI_Wtime()-tt0
1083 t_MDsetup = tcpu()-tt0
1086 ! Entering the MD loop
1092 if (lang.eq.2 .or. lang.eq.3) then
1096 call sd_verlet_p_setup
1098 call sd_verlet_ciccotti_setup
1102 pfric0_mat(i,j,0)=pfric_mat(i,j)
1103 afric0_mat(i,j,0)=afric_mat(i,j)
1104 vfric0_mat(i,j,0)=vfric_mat(i,j)
1105 prand0_mat(i,j,0)=prand_mat(i,j)
1106 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
1107 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
1110 flag_stoch(0)=.true.
1111 do i=1,maxflag_stoch
1112 flag_stoch(i)=.false.
1116 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
1118 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1122 else if (lang.eq.1 .or. lang.eq.4) then
1123 print *,"before setup_fricmat"
1125 print *,"after setup_fricmat"
1128 t_langsetup=MPI_Wtime()-tt0
1131 t_langsetup=tcpu()-tt0
1134 do itime=1,n_timestep
1135 if (large) print *,itime,ntwe
1137 if (large.and. mod(itime,ntwe).eq.0) &
1138 write (iout,*) "itime",itime
1140 if (lang.gt.0 .and. surfarea .and. &
1141 mod(itime,reset_fricmat).eq.0) then
1142 if (lang.eq.2 .or. lang.eq.3) then
1146 call sd_verlet_p_setup
1148 call sd_verlet_ciccotti_setup
1152 pfric0_mat(i,j,0)=pfric_mat(i,j)
1153 afric0_mat(i,j,0)=afric_mat(i,j)
1154 vfric0_mat(i,j,0)=vfric_mat(i,j)
1155 prand0_mat(i,j,0)=prand_mat(i,j)
1156 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
1157 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
1160 flag_stoch(0)=.true.
1161 do i=1,maxflag_stoch
1162 flag_stoch(i)=.false.
1165 else if (lang.eq.1 .or. lang.eq.4) then
1166 print *,"before setup_fricmat"
1168 print *,"after setup_fricmat"
1170 write (iout,'(a,i10)') &
1171 "Friction matrix reset based on surface area, itime",itime
1173 if (reset_vel .and. tbf .and. lang.eq.0 &
1174 .and. mod(itime,count_reset_vel).eq.0) then
1176 write(iout,'(a,f20.2)') &
1177 "Velocities reset to random values, time",totT
1180 d_t_old(j,i)=d_t(j,i)
1184 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
1188 d_t(j,0)=d_t(j,0)-vcm(j)
1191 kinetic_T=2.0d0/(dimen3*Rb)*EK
1192 scalfac=dsqrt(T_bath/kinetic_T)
1193 write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
1196 d_t_old(j,i)=scalfac*d_t(j,i)
1202 ! Time-reversible RESPA algorithm
1203 ! (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
1204 call RESPA_step(itime)
1206 ! Variable time step algorithm.
1207 if (large) print *,"before verlet_step"
1208 call velverlet_step(itime)
1209 if (large) print *,"after verlet_step"
1213 call brown_step(itime)
1215 print *,"Brown dynamics not here!"
1217 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1224 if (mod(itime,ntwe).eq.0) then
1228 ! call check_ecartint
1238 v_work(ind)=d_t(j,i)
1243 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum).and.mnum.lt.4) then
1246 v_work(ind)=d_t(j,i+nres)
1251 write (66,'(80f10.5)') &
1252 ((d_t(j,i),j=1,3),i=0,nres-1),((d_t(j,i+nres),j=1,3),i=1,nres)
1256 v_transf(i)=v_transf(i)+gvec(j,i)*v_work(j)
1258 v_transf(i)= v_transf(i)*dsqrt(geigen(i))
1260 write (67,'(80f10.5)') (v_transf(i),i=1,ind)
1263 if (mod(itime,ntwx).eq.0) then
1265 call enerprint(potEcomp)
1267 write (tytul,'("time",f8.2)') totT
1269 call hairpin(.true.,nharp,iharp)
1270 call secondary2(.true.)
1271 call pdbout(potE,tytul,ipdb)
1272 ! call enerprint(potEcomp)
1277 write(iout,*) "starting fodstep"
1278 call fodstep(nfodstep)
1279 write(iout,*) "after fodstep"
1282 call hairpin(.true.,nharp,iharp)
1283 call secondary2(.true.)
1284 call pdbout(potE,tytul,ipdb)
1291 if (rstcount.eq.1000.or.itime.eq.n_timestep) then
1292 open(irest2,file=rest2name,status='unknown')
1293 write(irest2,*) totT,EK,potE,totE,t_bath
1295 ! AL 4/17/17: Now writing d_t(0,:) too
1297 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
1299 ! AL 4/17/17: Now writing d_c(0,:) too
1301 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
1309 t_MD=MPI_Wtime()-tt0
1313 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))') &
1315 'MD calculations setup:',t_MDsetup,&
1316 'Energy & gradient evaluation:',t_enegrad,&
1317 'Stochastic MD setup:',t_langsetup,&
1318 'Stochastic MD step setup:',t_sdsetup,&
1320 write (iout,'(/28(1h=),a25,27(1h=))') &
1321 ' End of MD calculation '
1323 write (iout,*) "time for etotal",t_etotal," elong",t_elong,&
1325 write (iout,*) "time_fric",time_fric," time_stoch",time_stoch,&
1326 " time_fricmatmult",time_fricmatmult," time_fsample ",&
1331 !-----------------------------------------------------------------------------
1332 subroutine velverlet_step(itime)
1333 !-------------------------------------------------------------------------------
1334 ! Perform a single velocity Verlet step; the time step can be rescaled if
1335 ! increments in accelerations exceed the threshold
1336 !-------------------------------------------------------------------------------
1337 ! implicit real*8 (a-h,o-z)
1338 ! include 'DIMENSIONS'
1340 use control, only:tcpu
1342 use minimm, only:minim_dc
1345 integer :: ierror,ierrcode
1346 real(kind=8) :: errcode
1348 ! include 'COMMON.SETUP'
1349 ! include 'COMMON.VAR'
1350 ! include 'COMMON.MD'
1352 ! include 'COMMON.LANGEVIN'
1354 ! include 'COMMON.LANGEVIN.lang0'
1356 ! include 'COMMON.CHAIN'
1357 ! include 'COMMON.DERIV'
1358 ! include 'COMMON.GEO'
1359 ! include 'COMMON.LOCAL'
1360 ! include 'COMMON.INTERACT'
1361 ! include 'COMMON.IOUNITS'
1362 ! include 'COMMON.NAMES'
1363 ! include 'COMMON.TIME1'
1364 ! include 'COMMON.MUCA'
1365 real(kind=8),dimension(3) :: vcm,incr
1366 real(kind=8),dimension(3) :: L
1367 integer :: count,rstcount !ilen,
1369 character(len=50) :: tytul
1370 integer :: maxcount_scale = 30
1371 !el common /gucio/ cm
1372 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1373 !el common /stochcalc/ stochforcvec
1374 integer :: icount_scale,itime_scal,i,j,ifac_time,iretcode,itime
1376 real(kind=8) :: epdrift,tt0,fac_time
1378 if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres)) !(MAXRES6) maxres6=6*maxres
1384 if (large) print *,"after sddir_precalc"
1385 else if (lang.eq.2 .or. lang.eq.3) then
1387 call stochastic_force(stochforcvec)
1390 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
1392 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1399 icount_scale=icount_scale+1
1400 ! write(iout,*) "icount_scale",icount_scale,scalel
1401 if (icount_scale.gt.maxcount_scale) then
1403 "ERROR: too many attempts at scaling down the time step. ",&
1404 "amax=",amax,"epdrift=",epdrift,&
1405 "damax=",damax,"edriftmax=",edriftmax,&
1409 call MPI_Abort(MPI_COMM_WORLD,IERROR,IERRCODE)
1413 ! First step of the velocity Verlet algorithm
1418 else if (lang.eq.3) then
1420 call sd_verlet1_ciccotti
1422 else if (lang.eq.1) then
1427 ! Build the chain from the newly calculated coordinates
1428 call chainbuild_cart
1429 if (rattle) call rattle1
1431 if (large) then !.and. mod(itime,ntwe).eq.0) then
1432 write (iout,*) "Cartesian and internal coordinates: step 1"
1437 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),&
1438 (dc(j,i+nres),j=1,3)
1440 write (iout,*) "Accelerations"
1442 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1443 (d_a(j,i+nres),j=1,3)
1445 write (iout,*) "Velocities, step 1"
1447 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1448 (d_t(j,i+nres),j=1,3)
1457 ! Calculate energy and forces
1459 call etotal(potEcomp)
1460 ! AL 4/17/17: Reduce the steps if NaNs occurred.
1461 if (potEcomp(0).gt.0.99e18 .or. isnan(potEcomp(0)).gt.0) then
1462 call enerprint(potEcomp)
1464 if (icount_scale.gt.15) then
1465 write (iout,*) "Tu jest problem",potEcomp(0),d_time
1466 ! call gen_rand_conf(1,*335)
1467 ! call minim_dc(potEcomp(0),iretcode,100)
1470 ! call etotal(potEcomp)
1471 ! write(iout,*) "needed to repara,",potEcomp
1474 ! 335 write(iout,*) "Failed genrand"
1478 if (large.and. mod(itime,ntwe).eq.0) &
1479 call enerprint(potEcomp)
1482 t_etotal=t_etotal+MPI_Wtime()-tt0
1484 t_etotal=t_etotal+tcpu()-tt0
1487 potE=potEcomp(0)-potEcomp(51)
1489 ! Get the new accelerations
1492 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1494 t_enegrad=t_enegrad+tcpu()-tt0
1496 ! Determine maximum acceleration and scale down the timestep if needed
1498 amax=amax/(itime_scal+1)**2
1499 call predict_edrift(epdrift)
1500 ! write(iout,*) "amax=",amax,damax,epdrift,edriftmax,amax/(itime_scal+1)
1502 ! write (iout,*) "before enter if",scalel,icount_scale
1503 if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then
1504 ! write(iout,*) "I enter if"
1505 ! Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
1507 ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax)) &
1509 itime_scal=itime_scal+ifac_time
1510 ! fac_time=dmin1(damax/amax,0.5d0)
1511 fac_time=0.5d0**ifac_time
1512 d_time=d_time*fac_time
1513 if (lang.eq.2 .or. lang.eq.3) then
1515 ! write (iout,*) "Calling sd_verlet_setup: 1"
1516 ! Rescale the stochastic forces and recalculate or restore
1517 ! the matrices of tinker integrator
1518 if (itime_scal.gt.maxflag_stoch) then
1519 if (large) write (iout,'(a,i5,a)') &
1520 "Calculate matrices for stochastic step;",&
1521 " itime_scal ",itime_scal
1523 call sd_verlet_p_setup
1525 call sd_verlet_ciccotti_setup
1527 write (iout,'(2a,i3,a,i3,1h.)') &
1528 "Warning: cannot store matrices for stochastic",&
1529 " integration because the index",itime_scal,&
1530 " is greater than",maxflag_stoch
1531 write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",&
1532 " integration Langevin algorithm for better efficiency."
1533 else if (flag_stoch(itime_scal)) then
1534 if (large) write (iout,'(a,i5,a,l1)') &
1535 "Restore matrices for stochastic step; itime_scal ",&
1536 itime_scal," flag ",flag_stoch(itime_scal)
1539 pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
1540 afric_mat(i,j)=afric0_mat(i,j,itime_scal)
1541 vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
1542 prand_mat(i,j)=prand0_mat(i,j,itime_scal)
1543 vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
1544 vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
1548 if (large) write (iout,'(2a,i5,a,l1)') &
1549 "Calculate & store matrices for stochastic step;",&
1550 " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
1552 call sd_verlet_p_setup
1554 call sd_verlet_ciccotti_setup
1556 flag_stoch(ifac_time)=.true.
1559 pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
1560 afric0_mat(i,j,itime_scal)=afric_mat(i,j)
1561 vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
1562 prand0_mat(i,j,itime_scal)=prand_mat(i,j)
1563 vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
1564 vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
1568 fac_time=1.0d0/dsqrt(fac_time)
1570 stochforcvec(i)=fac_time*stochforcvec(i)
1573 else if (lang.eq.1) then
1574 ! Rescale the accelerations due to stochastic forces
1575 fac_time=1.0d0/dsqrt(fac_time)
1577 d_as_work(i)=d_as_work(i)*fac_time
1580 if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)') &
1581 "itime",itime," Timestep scaled down to ",&
1582 d_time," ifac_time",ifac_time," itime_scal",itime_scal
1584 ! Second step of the velocity Verlet algorithm
1589 else if (lang.eq.3) then
1591 call sd_verlet2_ciccotti
1593 else if (lang.eq.1) then
1598 if (rattle) call rattle2
1601 if (d_time.ne.d_time0) then
1604 if (lang.eq.2 .or. lang.eq.3) then
1605 if (large) write (iout,'(a)') &
1606 "Restore original matrices for stochastic step"
1607 ! write (iout,*) "Calling sd_verlet_setup: 2"
1608 ! Restore the matrices of tinker integrator if the time step has been restored
1611 pfric_mat(i,j)=pfric0_mat(i,j,0)
1612 afric_mat(i,j)=afric0_mat(i,j,0)
1613 vfric_mat(i,j)=vfric0_mat(i,j,0)
1614 prand_mat(i,j)=prand0_mat(i,j,0)
1615 vrand_mat1(i,j)=vrand0_mat1(i,j,0)
1616 vrand_mat2(i,j)=vrand0_mat2(i,j,0)
1624 ! Calculate the kinetic and the total energy and the kinetic temperature
1628 ! call kinetic1(EK1)
1629 ! write (iout,*) "step",itime," EK",EK," EK1",EK1
1631 ! Couple the system to Berendsen bath if needed
1632 if (tbf .and. lang.eq.0) then
1635 kinetic_T=2.0d0/(dimen3*Rb)*EK
1636 ! Backup the coordinates, velocities, and accelerations
1640 d_t_old(j,i)=d_t(j,i)
1641 d_a_old(j,i)=d_a(j,i)
1645 if (mod(itime,ntwe).eq.0 .and. large) then
1646 write (iout,*) "Velocities, step 2"
1648 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1649 (d_t(j,i+nres),j=1,3)
1654 end subroutine velverlet_step
1655 !-----------------------------------------------------------------------------
1656 subroutine RESPA_step(itime)
1657 !-------------------------------------------------------------------------------
1658 ! Perform a single RESPA step.
1659 !-------------------------------------------------------------------------------
1660 ! implicit real*8 (a-h,o-z)
1661 ! include 'DIMENSIONS'
1665 use control, only:tcpu
1667 ! use io_conf, only:cartprint
1670 integer :: IERROR,ERRCODE
1672 ! include 'COMMON.SETUP'
1673 ! include 'COMMON.CONTROL'
1674 ! include 'COMMON.VAR'
1675 ! include 'COMMON.MD'
1677 ! include 'COMMON.LANGEVIN'
1679 ! include 'COMMON.LANGEVIN.lang0'
1681 ! include 'COMMON.CHAIN'
1682 ! include 'COMMON.DERIV'
1683 ! include 'COMMON.GEO'
1684 ! include 'COMMON.LOCAL'
1685 ! include 'COMMON.INTERACT'
1686 ! include 'COMMON.IOUNITS'
1687 ! include 'COMMON.NAMES'
1688 ! include 'COMMON.TIME1'
1689 real(kind=8),dimension(0:n_ene) :: energia_short,energia_long
1690 real(kind=8),dimension(3) :: L,vcm,incr
1691 real(kind=8),dimension(3,0:2*nres) :: dc_old0,d_t_old0,d_a_old0 !(3,0:maxres2) maxres2=2*maxres
1692 logical :: PRINT_AMTS_MSG = .false.
1693 integer :: count,rstcount !ilen,
1695 character(len=50) :: tytul
1696 integer :: maxcount_scale = 10
1697 !el common /gucio/ cm
1698 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1699 !el common /stochcalc/ stochforcvec
1700 integer :: itt,i,j,itsplit,itime
1702 !el common /cipiszcze/ itt
1704 real(kind=8) :: epdrift,tt0,epdriftmax
1707 if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres)) !(MAXRES6) maxres6=6*maxres
1711 if (large.and. mod(itime,ntwe).eq.0) then
1712 write (iout,*) "***************** RESPA itime",itime
1713 write (iout,*) "Cartesian and internal coordinates: step 0"
1715 call pdbout(0.0d0,"cipiszcze",iout)
1717 write (iout,*) "Accelerations from long-range forces"
1719 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1720 (d_a(j,i+nres),j=1,3)
1722 write (iout,*) "Velocities, step 0"
1724 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1725 (d_t(j,i+nres),j=1,3)
1730 ! Perform the initial RESPA step (increment velocities)
1731 ! write (iout,*) "*********************** RESPA ini"
1734 if (mod(itime,ntwe).eq.0 .and. large) then
1735 write (iout,*) "Velocities, end"
1737 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1738 (d_t(j,i+nres),j=1,3)
1742 ! Compute the short-range forces
1748 ! 7/2/2009 commented out
1750 ! call etotal_short(energia_short)
1753 ! 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
1756 d_a(j,i)=d_a_short(j,i)
1760 if (large.and. mod(itime,ntwe).eq.0) then
1761 write (iout,*) "energia_short",energia_short(0)
1762 write (iout,*) "Accelerations from short-range forces"
1764 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1765 (d_a(j,i+nres),j=1,3)
1770 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1772 t_enegrad=t_enegrad+tcpu()-tt0
1777 d_t_old(j,i)=d_t(j,i)
1778 d_a_old(j,i)=d_a(j,i)
1781 ! 6/30/08 A-MTS: attempt at increasing the split number
1784 dc_old0(j,i)=dc_old(j,i)
1785 d_t_old0(j,i)=d_t_old(j,i)
1786 d_a_old0(j,i)=d_a_old(j,i)
1789 if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
1790 if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
1797 ! write (iout,*) "itime",itime," ntime_split",ntime_split
1798 ! Split the time step
1799 d_time=d_time0/ntime_split
1800 ! Perform the short-range RESPA steps (velocity Verlet increments of
1801 ! positions and velocities using short-range forces)
1802 ! write (iout,*) "*********************** RESPA split"
1803 do itsplit=1,ntime_split
1806 else if (lang.eq.2 .or. lang.eq.3) then
1808 call stochastic_force(stochforcvec)
1811 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
1813 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1818 ! First step of the velocity Verlet algorithm
1823 else if (lang.eq.3) then
1825 call sd_verlet1_ciccotti
1827 else if (lang.eq.1) then
1832 ! Build the chain from the newly calculated coordinates
1833 call chainbuild_cart
1834 if (rattle) call rattle1
1836 if (large.and. mod(itime,ntwe).eq.0) then
1837 write (iout,*) "***** ITSPLIT",itsplit
1838 write (iout,*) "Cartesian and internal coordinates: step 1"
1839 call pdbout(0.0d0,"cipiszcze",iout)
1842 write (iout,*) "Velocities, step 1"
1844 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1845 (d_t(j,i+nres),j=1,3)
1854 ! Calculate energy and forces
1856 call etotal_short(energia_short)
1857 ! AL 4/17/17: Exit itime_split loop when energy goes infinite
1858 if (energia_short(0).gt.0.99e20 .or. isnan(energia_short(0)) ) then
1859 if (PRINT_AMTS_MSG) &
1860 write (iout,*) "Infinities/NaNs in energia_short",energia_short(0),"; increasing ntime_split to",ntime_split
1861 ntime_split=ntime_split*2
1862 if (ntime_split.gt.maxtime_split) then
1865 "Cannot rescue the run; aborting job. Retry with a smaller time step"
1867 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1870 "Cannot rescue the run; terminating. Retry with a smaller time step"
1876 if (large.and. mod(itime,ntwe).eq.0) &
1877 call enerprint(energia_short)
1880 t_eshort=t_eshort+MPI_Wtime()-tt0
1882 t_eshort=t_eshort+tcpu()-tt0
1886 ! Get the new accelerations
1888 ! 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
1891 d_a_short(j,i)=d_a(j,i)
1895 if (large.and. mod(itime,ntwe).eq.0) then
1896 write (iout,*)"energia_short",energia_short(0)
1897 write (iout,*) "Accelerations from short-range forces"
1899 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1900 (d_a(j,i+nres),j=1,3)
1905 ! Determine maximum acceleration and scale down the timestep if needed
1907 amax=amax/ntime_split**2
1908 call predict_edrift(epdrift)
1909 if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0) &
1910 write (iout,*) "amax",amax," damax",damax,&
1911 " epdrift",epdrift," epdriftmax",epdriftmax
1912 ! Exit loop and try with increased split number if the change of
1913 ! acceleration is too big
1914 if (amax.gt.damax .or. epdrift.gt.edriftmax) then
1915 if (ntime_split.lt.maxtime_split) then
1917 ntime_split=ntime_split*2
1918 ! AL 4/17/17: We should exit the itime_split loop when acceleration change is too big
1922 dc_old(j,i)=dc_old0(j,i)
1923 d_t_old(j,i)=d_t_old0(j,i)
1924 d_a_old(j,i)=d_a_old0(j,i)
1927 if (PRINT_AMTS_MSG) then
1928 write (iout,*) "acceleration/energy drift too large",amax,&
1929 epdrift," split increased to ",ntime_split," itime",itime,&
1935 "Uh-hu. Bumpy landscape. Maximum splitting number",&
1937 " already reached!!! Trying to carry on!"
1941 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1943 t_enegrad=t_enegrad+tcpu()-tt0
1945 ! Second step of the velocity Verlet algorithm
1950 else if (lang.eq.3) then
1952 call sd_verlet2_ciccotti
1954 else if (lang.eq.1) then
1959 if (rattle) call rattle2
1960 ! Backup the coordinates, velocities, and accelerations
1964 d_t_old(j,i)=d_t(j,i)
1965 d_a_old(j,i)=d_a(j,i)
1972 ! Restore the time step
1974 ! Compute long-range forces
1981 call etotal_long(energia_long)
1982 if (energia_long(0).gt.0.99e20 .or. isnan(energia_long(0))) then
1985 "Infinitied/NaNs in energia_long, Aborting MPI job."
1987 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1989 write (iout,*) "Infinitied/NaNs in energia_long, terminating."
1993 if (large.and. mod(itime,ntwe).eq.0) &
1994 call enerprint(energia_long)
1997 t_elong=t_elong+MPI_Wtime()-tt0
1999 t_elong=t_elong+tcpu()-tt0
2002 potE=potEcomp(0)-potEcomp(51)
2006 t_enegrad=t_enegrad+MPI_Wtime()-tt0
2008 t_enegrad=t_enegrad+tcpu()-tt0
2010 ! Compute accelerations from long-range forces
2012 if (large.and. mod(itime,ntwe).eq.0) then
2013 write (iout,*) "energia_long",energia_long(0)
2014 write (iout,*) "Cartesian and internal coordinates: step 2"
2016 call pdbout(0.0d0,"cipiszcze",iout)
2018 write (iout,*) "Accelerations from long-range forces"
2020 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
2021 (d_a(j,i+nres),j=1,3)
2023 write (iout,*) "Velocities, step 2"
2025 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2026 (d_t(j,i+nres),j=1,3)
2030 ! Compute the final RESPA step (increment velocities)
2031 ! write (iout,*) "*********************** RESPA fin"
2033 ! Compute the complete potential energy
2035 potEcomp(i)=energia_short(i)+energia_long(i)
2037 potE=potEcomp(0)-potEcomp(51)
2038 ! potE=energia_short(0)+energia_long(0)
2041 ! Calculate the kinetic and the total energy and the kinetic temperature
2044 ! Couple the system to Berendsen bath if needed
2045 if (tbf .and. lang.eq.0) then
2048 kinetic_T=2.0d0/(dimen3*Rb)*EK
2049 ! Backup the coordinates, velocities, and accelerations
2051 if (mod(itime,ntwe).eq.0 .and. large) then
2052 write (iout,*) "Velocities, end"
2054 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2055 (d_t(j,i+nres),j=1,3)
2060 end subroutine RESPA_step
2061 !-----------------------------------------------------------------------------
2062 subroutine RESPA_vel
2063 ! First and last RESPA step (incrementing velocities using long-range
2066 ! implicit real*8 (a-h,o-z)
2067 ! include 'DIMENSIONS'
2068 ! include 'COMMON.CONTROL'
2069 ! include 'COMMON.VAR'
2070 ! include 'COMMON.MD'
2071 ! include 'COMMON.CHAIN'
2072 ! include 'COMMON.DERIV'
2073 ! include 'COMMON.GEO'
2074 ! include 'COMMON.LOCAL'
2075 ! include 'COMMON.INTERACT'
2076 ! include 'COMMON.IOUNITS'
2077 ! include 'COMMON.NAMES'
2078 integer :: i,j,inres,mnum
2081 d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
2085 d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
2090 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2091 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2092 .and.(mnum.lt.4)) then
2095 d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
2100 end subroutine RESPA_vel
2101 !-----------------------------------------------------------------------------
2103 ! Applying velocity Verlet algorithm - step 1 to coordinates
2105 ! implicit real*8 (a-h,o-z)
2106 ! include 'DIMENSIONS'
2107 ! include 'COMMON.CONTROL'
2108 ! include 'COMMON.VAR'
2109 ! include 'COMMON.MD'
2110 ! include 'COMMON.CHAIN'
2111 ! include 'COMMON.DERIV'
2112 ! include 'COMMON.GEO'
2113 ! include 'COMMON.LOCAL'
2114 ! include 'COMMON.INTERACT'
2115 ! include 'COMMON.IOUNITS'
2116 ! include 'COMMON.NAMES'
2117 real(kind=8) :: adt,adt2
2118 integer :: i,j,inres,mnum
2121 write (iout,*) "VELVERLET1 START: DC"
2123 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),&
2124 (dc(j,i+nres),j=1,3)
2128 adt=d_a_old(j,0)*d_time
2130 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
2131 d_t_new(j,0)=d_t_old(j,0)+adt2
2132 d_t(j,0)=d_t_old(j,0)+adt
2136 adt=d_a_old(j,i)*d_time
2138 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
2139 d_t_new(j,i)=d_t_old(j,i)+adt2
2140 d_t(j,i)=d_t_old(j,i)+adt
2145 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2146 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2147 .and.(mnum.lt.4)) then
2150 adt=d_a_old(j,inres)*d_time
2152 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
2153 d_t_new(j,inres)=d_t_old(j,inres)+adt2
2154 d_t(j,inres)=d_t_old(j,inres)+adt
2159 write (iout,*) "VELVERLET1 END: DC"
2161 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),&
2162 (dc(j,i+nres),j=1,3)
2166 end subroutine verlet1
2167 !-----------------------------------------------------------------------------
2169 ! Step 2 of the velocity Verlet algorithm: update velocities
2171 ! implicit real*8 (a-h,o-z)
2172 ! include 'DIMENSIONS'
2173 ! include 'COMMON.CONTROL'
2174 ! include 'COMMON.VAR'
2175 ! include 'COMMON.MD'
2176 ! include 'COMMON.CHAIN'
2177 ! include 'COMMON.DERIV'
2178 ! include 'COMMON.GEO'
2179 ! include 'COMMON.LOCAL'
2180 ! include 'COMMON.INTERACT'
2181 ! include 'COMMON.IOUNITS'
2182 ! include 'COMMON.NAMES'
2183 integer :: i,j,inres,mnum
2186 d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
2190 d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
2195 ! iti=iabs(itype(i,mnum))
2196 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2197 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2198 .and.(mnum.lt.4)) then
2201 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
2206 end subroutine verlet2
2207 !-----------------------------------------------------------------------------
2208 subroutine sddir_precalc
2209 ! Applying velocity Verlet algorithm - step 1 to coordinates
2210 ! implicit real*8 (a-h,o-z)
2211 ! include 'DIMENSIONS'
2217 ! include 'COMMON.CONTROL'
2218 ! include 'COMMON.VAR'
2219 ! include 'COMMON.MD'
2221 ! include 'COMMON.LANGEVIN'
2223 ! include 'COMMON.LANGEVIN.lang0'
2225 ! include 'COMMON.CHAIN'
2226 ! include 'COMMON.DERIV'
2227 ! include 'COMMON.GEO'
2228 ! include 'COMMON.LOCAL'
2229 ! include 'COMMON.INTERACT'
2230 ! include 'COMMON.IOUNITS'
2231 ! include 'COMMON.NAMES'
2232 ! include 'COMMON.TIME1'
2233 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
2234 !el common /stochcalc/ stochforcvec
2235 real(kind=8) :: time00
2238 ! Compute friction and stochastic forces
2242 if (large) print *,"before friction_force"
2244 if (large) print *,"after friction_force"
2245 time_fric=time_fric+MPI_Wtime()-time00
2247 call stochastic_force(stochforcvec)
2248 time_stoch=time_stoch+MPI_Wtime()-time00
2251 ! Compute the acceleration due to friction forces (d_af_work) and stochastic
2252 ! forces (d_as_work)
2254 ! call ginv_mult(fric_work, d_af_work)
2255 ! call ginv_mult(stochforcvec, d_as_work)
2257 call fivediaginv_mult(dimen,fric_work, d_af_work)
2258 call fivediaginv_mult(dimen,stochforcvec, d_as_work)
2260 write(iout,*),"dimen",dimen
2262 write(iout,*) "fricwork",fric_work(i), d_af_work(i)
2263 write(iout,*) "stochforcevec", stochforcvec(i), d_as_work(i)
2267 call ginv_mult(fric_work, d_af_work)
2268 call ginv_mult(stochforcvec, d_as_work)
2272 end subroutine sddir_precalc
2273 !-----------------------------------------------------------------------------
2274 subroutine sddir_verlet1
2275 ! Applying velocity Verlet algorithm - step 1 to velocities
2278 ! implicit real*8 (a-h,o-z)
2279 ! include 'DIMENSIONS'
2280 ! include 'COMMON.CONTROL'
2281 ! include 'COMMON.VAR'
2282 ! include 'COMMON.MD'
2284 ! include 'COMMON.LANGEVIN'
2286 ! include 'COMMON.LANGEVIN.lang0'
2288 ! include 'COMMON.CHAIN'
2289 ! include 'COMMON.DERIV'
2290 ! include 'COMMON.GEO'
2291 ! include 'COMMON.LOCAL'
2292 ! include 'COMMON.INTERACT'
2293 ! include 'COMMON.IOUNITS'
2294 ! include 'COMMON.NAMES'
2295 ! Revised 3/31/05 AL: correlation between random contributions to
2296 ! position and velocity increments included.
2297 real(kind=8) :: sqrt13 = 0.57735026918962576451d0 ! 1/sqrt(3)
2298 real(kind=8) :: adt,adt2
2299 integer :: i,j,ind,inres,mnum
2301 ! Add the contribution from BOTH friction and stochastic force to the
2302 ! coordinates, but ONLY the contribution from the friction forces to velocities
2305 adt=(d_a_old(j,0)+d_af_work(j))*d_time
2306 adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
2307 ! write(iout,*) i,"adt",adt,"ads",adt2,d_a_old(j,0),d_af_work(j),d_time
2308 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
2309 d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
2310 d_t(j,0)=d_t_old(j,0)+adt
2315 adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
2316 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
2317 ! write(iout,*) i,"adt",adt,"ads",adt2,d_a_old(j,i),d_af_work(ind+j)
2318 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
2319 d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
2320 d_t(j,i)=d_t_old(j,i)+adt
2326 ! iti=iabs(itype(i,mnum))
2327 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2328 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2329 .and.(mnum.lt.4)) then
2332 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
2333 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
2334 ! write(iout,*) i,"adt",adt,"ads",adt2,d_a_old(j,inres),d_af_work(ind+j)
2335 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
2336 d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
2337 d_t(j,inres)=d_t_old(j,inres)+adt
2344 end subroutine sddir_verlet1
2345 !-----------------------------------------------------------------------------
2346 subroutine sddir_verlet2
2347 ! Calculating the adjusted velocities for accelerations
2350 ! implicit real*8 (a-h,o-z)
2351 ! include 'DIMENSIONS'
2352 ! include 'COMMON.CONTROL'
2353 ! include 'COMMON.VAR'
2354 ! include 'COMMON.MD'
2356 ! include 'COMMON.LANGEVIN'
2358 ! include 'COMMON.LANGEVIN.lang0'
2360 ! include 'COMMON.CHAIN'
2361 ! include 'COMMON.DERIV'
2362 ! include 'COMMON.GEO'
2363 ! include 'COMMON.LOCAL'
2364 ! include 'COMMON.INTERACT'
2365 ! include 'COMMON.IOUNITS'
2366 ! include 'COMMON.NAMES'
2367 real(kind=8),dimension(6*nres) :: stochforcvec,d_as_work1 !(MAXRES6) maxres6=6*maxres
2368 real(kind=8) :: cos60 = 0.5d0, sin60 = 0.86602540378443864676d0
2369 integer :: i,j,ind,inres,mnum
2370 ! Revised 3/31/05 AL: correlation between random contributions to
2371 ! position and velocity increments included.
2372 ! The correlation coefficients are calculated at low-friction limit.
2373 ! Also, friction forces are now not calculated with new velocities.
2375 ! call friction_force
2376 call stochastic_force(stochforcvec)
2378 ! Compute the acceleration due to friction forces (d_af_work) and stochastic
2379 ! forces (d_as_work)
2382 call fivediaginv_mult(6*nres,stochforcvec, d_as_work1)
2384 call ginv_mult(stochforcvec, d_as_work1)
2391 d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j)) &
2392 +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
2397 d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j)) &
2398 +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
2404 ! iti=iabs(itype(i,mnum))
2405 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2406 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2407 .and.(mnum.lt.4)) then
2410 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres) &
2411 +d_af_work(ind+j))+sin60*d_as_work(ind+j) &
2412 +cos60*d_as_work1(ind+j))*d_time
2418 end subroutine sddir_verlet2
2419 !-----------------------------------------------------------------------------
2420 subroutine max_accel
2422 ! Find the maximum difference in the accelerations of the the sites
2423 ! at the beginning and the end of the time step.
2426 ! implicit real*8 (a-h,o-z)
2427 ! include 'DIMENSIONS'
2428 ! include 'COMMON.CONTROL'
2429 ! include 'COMMON.VAR'
2430 ! include 'COMMON.MD'
2431 ! include 'COMMON.CHAIN'
2432 ! include 'COMMON.DERIV'
2433 ! include 'COMMON.GEO'
2434 ! include 'COMMON.LOCAL'
2435 ! include 'COMMON.INTERACT'
2436 ! include 'COMMON.IOUNITS'
2437 real(kind=8),dimension(3) :: aux,accel,accel_old
2438 real(kind=8) :: dacc
2442 ! aux(j)=d_a(j,0)-d_a_old(j,0)
2443 accel_old(j)=d_a_old(j,0)
2450 ! 7/3/08 changed to asymmetric difference
2452 ! accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
2453 accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
2454 accel(j)=accel(j)+0.5d0*d_a(j,i)
2455 ! if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
2456 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
2457 dacc=dabs(accel(j)-accel_old(j))
2458 ! write (iout,*) i,dacc
2459 if (dacc.gt.amax) amax=dacc
2467 accel_old(j)=d_a_old(j,0)
2472 accel_old(j)=accel_old(j)+d_a_old(j,1)
2473 accel(j)=accel(j)+d_a(j,1)
2478 ! iti=iabs(itype(i,mnum))
2479 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2480 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2481 .and.(mnum.lt.4)) then
2483 ! accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
2484 accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
2485 accel(j)=accel(j)+d_a(j,i+nres)
2489 ! if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
2490 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
2491 dacc=dabs(accel(j)-accel_old(j))
2492 ! write (iout,*) "side-chain",i,dacc
2493 if (dacc.gt.amax) amax=dacc
2497 accel_old(j)=accel_old(j)+d_a_old(j,i)
2498 accel(j)=accel(j)+d_a(j,i)
2499 ! aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
2503 end subroutine max_accel
2504 !-----------------------------------------------------------------------------
2505 subroutine predict_edrift(epdrift)
2507 ! Predict the drift of the potential energy
2510 use control_data, only: lmuca
2511 ! implicit real*8 (a-h,o-z)
2512 ! include 'DIMENSIONS'
2513 ! include 'COMMON.CONTROL'
2514 ! include 'COMMON.VAR'
2515 ! include 'COMMON.MD'
2516 ! include 'COMMON.CHAIN'
2517 ! include 'COMMON.DERIV'
2518 ! include 'COMMON.GEO'
2519 ! include 'COMMON.LOCAL'
2520 ! include 'COMMON.INTERACT'
2521 ! include 'COMMON.IOUNITS'
2522 ! include 'COMMON.MUCA'
2523 real(kind=8) :: epdrift,epdriftij
2525 ! Drift of the potential energy
2531 epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
2532 if (lmuca) epdriftij=epdriftij*factor
2533 ! write (iout,*) "back",i,j,epdriftij
2534 if (epdriftij.gt.epdrift) epdrift=epdriftij
2538 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1.and.&
2539 molnum(i).lt.4) then
2542 dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
2543 if (lmuca) epdriftij=epdriftij*factor
2544 ! write (iout,*) "side",i,j,epdriftij
2545 if (epdriftij.gt.epdrift) epdrift=epdriftij
2549 epdrift=0.5d0*epdrift*d_time*d_time
2550 ! write (iout,*) "epdrift",epdrift
2552 end subroutine predict_edrift
2553 !-----------------------------------------------------------------------------
2554 subroutine verlet_bath
2556 ! Coupling to the thermostat by using the Berendsen algorithm
2559 ! implicit real*8 (a-h,o-z)
2560 ! include 'DIMENSIONS'
2561 ! include 'COMMON.CONTROL'
2562 ! include 'COMMON.VAR'
2563 ! include 'COMMON.MD'
2564 ! include 'COMMON.CHAIN'
2565 ! include 'COMMON.DERIV'
2566 ! include 'COMMON.GEO'
2567 ! include 'COMMON.LOCAL'
2568 ! include 'COMMON.INTERACT'
2569 ! include 'COMMON.IOUNITS'
2570 ! include 'COMMON.NAMES'
2571 real(kind=8) :: T_half,fact
2572 integer :: i,j,inres,mnum
2574 T_half=2.0d0/(dimen3*Rb)*EK
2575 fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
2576 ! write(iout,*) "T_half", T_half
2577 ! write(iout,*) "EK", EK
2578 ! write(iout,*) "fact", fact
2580 d_t(j,0)=fact*d_t(j,0)
2584 d_t(j,i)=fact*d_t(j,i)
2589 ! iti=iabs(itype(i,mnum))
2590 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2591 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2592 .and.(mnum.lt.4)) then
2595 d_t(j,inres)=fact*d_t(j,inres)
2600 end subroutine verlet_bath
2601 !-----------------------------------------------------------------------------
2603 ! Set up the initial conditions of a MD simulation
2606 use control, only:tcpu
2607 !el use io_basic, only:ilen
2610 use minimm, only:minim_dc,minimize,sc_move
2611 use io_config, only:readrst
2612 use io, only:statout
2613 use random, only: iran_num
2614 ! implicit real*8 (a-h,o-z)
2615 ! include 'DIMENSIONS'
2618 character(len=16) :: form
2619 integer :: IERROR,ERRCODE
2621 integer :: iranmin,itrial,itmp,n_model_try,k, &
2623 integer, dimension(:),allocatable :: list_model_try
2624 integer, dimension(0:nodes-1) :: i_start_models
2625 ! include 'COMMON.SETUP'
2626 ! include 'COMMON.CONTROL'
2627 ! include 'COMMON.VAR'
2628 ! include 'COMMON.MD'
2630 ! include 'COMMON.LANGEVIN'
2632 ! include 'COMMON.LANGEVIN.lang0'
2634 ! include 'COMMON.CHAIN'
2635 ! include 'COMMON.DERIV'
2636 ! include 'COMMON.GEO'
2637 ! include 'COMMON.LOCAL'
2638 ! include 'COMMON.INTERACT'
2639 ! include 'COMMON.IOUNITS'
2640 ! include 'COMMON.NAMES'
2641 ! include 'COMMON.REMD'
2642 real(kind=8),dimension(0:n_ene) :: energia_long,energia_short,energia
2643 real(kind=8),dimension(3) :: vcm,incr,L
2644 real(kind=8) :: xv,sigv,lowb,highb
2645 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
2646 character(len=256) :: qstr
2649 character(len=50) :: tytul
2650 logical :: file_exist
2651 !el common /gucio/ cm
2652 integer :: i,j,ipos,iq,iw,nft_sc,iretcode,nfun,ierr,mnum,itime
2653 real(kind=8) :: etot,tt0
2657 ! write(iout,*) "d_time", d_time
2658 ! Compute the standard deviations of stochastic forces for Langevin dynamics
2659 ! if the friction coefficients do not depend on surface area
2660 if (lang.gt.0 .and. .not.surfarea) then
2663 stdforcp(i)=stdfp(mnum)*dsqrt(gamp(mnum))
2667 stdforcsc(i)=stdfsc(iabs(itype(i,mnum)),mnum) &
2668 *dsqrt(gamsc(iabs(itype(i,mnum)),mnum))
2672 ! Open the pdb file for snapshotshots
2675 if (ilen(tmpdir).gt.0) &
2676 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2677 liczba(:ilen(liczba))//".pdb")
2679 file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2683 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) &
2684 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2685 liczba(:ilen(liczba))//".x")
2686 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2689 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) &
2690 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2691 liczba(:ilen(liczba))//".cx")
2692 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2698 if (ilen(tmpdir).gt.0) &
2699 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
2700 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
2702 if (ilen(tmpdir).gt.0) &
2703 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
2704 cartname=prefix(:ilen(prefix))//"_MD.cx"
2708 write (qstr,'(256(1h ))')
2711 iq = qinfrag(i,iset)*10
2712 iw = wfrag(i,iset)/100
2714 if(me.eq.king.or..not.out1file) &
2715 write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
2716 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
2721 iq = qinpair(i,iset)*10
2722 iw = wpair(i,iset)/100
2724 if(me.eq.king.or..not.out1file) &
2725 write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
2726 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
2730 ! pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
2732 ! cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
2734 ! cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
2736 ! statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
2740 if (restart1file) then
2742 inquire(file=mremd_rst_name,exist=file_exist)
2743 write (*,*) me," Before broadcast: file_exist",file_exist
2745 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,&
2748 write (*,*) me," After broadcast: file_exist",file_exist
2749 ! inquire(file=mremd_rst_name,exist=file_exist)
2750 if(me.eq.king.or..not.out1file) &
2751 write(iout,*) "Initial state read by master and distributed"
2753 if (ilen(tmpdir).gt.0) &
2754 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_' &
2755 //liczba(:ilen(liczba))//'.rst')
2756 inquire(file=rest2name,exist=file_exist)
2759 if(.not.restart1file) then
2760 if(me.eq.king.or..not.out1file) &
2761 write(iout,*) "Initial state will be read from file ",&
2762 rest2name(:ilen(rest2name))
2765 call rescale_weights(t_bath)
2767 if(me.eq.king.or..not.out1file)then
2768 if (restart1file) then
2769 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),&
2772 write(iout,*) "File ",rest2name(:ilen(rest2name)),&
2775 write(iout,*) "Initial velocities randomly generated"
2782 ! Generate initial velocities
2783 if(me.eq.king.or..not.out1file) &
2784 write(iout,*) "Initial velocities randomly generated"
2789 ! rest2name = prefix(:ilen(prefix))//'.rst'
2790 if(me.eq.king.or..not.out1file)then
2791 write (iout,*) "Initial velocities"
2793 write (iout,'(i6,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2794 (d_t(j,i+nres),j=1,3)
2796 ! Zeroing the total angular momentum of the system
2797 write(iout,*) "Calling the zero-angular momentum subroutine"
2800 ! Getting the potential energy and forces and velocities and accelerations
2802 ! write (iout,*) "velocity of the center of the mass:"
2803 ! write (iout,*) (vcm(j),j=1,3)
2805 d_t(j,0)=d_t(j,0)-vcm(j)
2807 ! Removing the velocity of the center of mass
2809 if(me.eq.king.or..not.out1file)then
2810 write (iout,*) "vcm right after adjustment:"
2811 write (iout,*) (vcm(j),j=1,3)
2818 if ((.not.rest).or.(forceminim)) then
2819 if (forceminim) call chainbuild_cart
2821 if(iranconf.ne.0 .or.indpdb.gt.0.and..not.unres_pdb .or.preminim) then
2823 print *, 'Calling OVERLAP_SC'
2824 call overlap_sc(fail)
2825 print *,'after OVERLAP'
2828 print *,'call SC_MOVE'
2829 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
2830 print *,'SC_move',nft_sc,etot
2831 if(me.eq.king.or..not.out1file) &
2832 write(iout,*) 'SC_move',nft_sc,etot
2836 print *, 'Calling MINIM_DC'
2837 call minim_dc(etot,iretcode,nfun)
2839 call geom_to_var(nvar,varia)
2840 print *,'Calling MINIMIZE.'
2841 call minimize(etot,varia,iretcode,nfun)
2842 call var_to_geom(nvar,varia)
2844 write(iout,*) "just before minimin"
2846 if(me.eq.king.or..not.out1file) &
2847 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2849 write(iout,*) "just after minimin"
2851 if(iranconf.ne.0) then
2852 !c 8/22/17 AL Loop to produce a low-energy random conformation
2855 if(me.eq.king.or..not.out1file) &
2856 write (iout,*) 'Calling OVERLAP_SC'
2857 call overlap_sc(fail)
2858 endif !endif overlap
2861 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
2862 print *,'SC_move',nft_sc,etot
2863 if(me.eq.king.or..not.out1file) &
2864 write(iout,*) 'SC_move',nft_sc,etot
2868 print *, 'Calling MINIM_DC'
2869 call minim_dc(etot,iretcode,nfun)
2870 call int_from_cart1(.false.)
2872 call geom_to_var(nvar,varia)
2873 print *,'Calling MINIMIZE.'
2874 call minimize(etot,varia,iretcode,nfun)
2875 call var_to_geom(nvar,varia)
2877 if(me.eq.king.or..not.out1file) &
2878 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2879 write(iout,*) "just after minimin"
2881 if (isnan(etot) .or. etot.gt.4.0d6) then
2882 write (iout,*) "Energy too large",etot, &
2883 " trying another random conformation"
2886 call gen_rand_conf(itmp,*30)
2888 30 write (iout,*) 'Failed to generate random conformation', &
2890 write (*,*) 'Processor:',me, &
2891 ' Failed to generate random conformation',&
2900 write (iout,'(a,i3,a)') 'Processor:',me, &
2901 ' error in generating random conformation.'
2902 write (*,'(a,i3,a)') 'Processor:',me, &
2903 ' error in generating random conformation.'
2906 ! call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
2907 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
2917 write (iout,'(a,i3,a)') 'Processor:',me, &
2918 ' failed to generate a low-energy random conformation.'
2919 write (*,'(a,i3,a,f10.3)') 'Processor:',me, &
2920 ' failed to generate a low-energy random conformation.',etot
2924 ! call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
2925 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
2930 else if (preminim) then
2931 if (start_from_model) then
2935 do while (fail .and. n_model_try.lt.nmodel_start)
2936 write (iout,*) "n_model_try",n_model_try
2938 i_model=iran_num(1,nmodel_start)
2940 if (i_model.eq.list_model_try(k)) exit
2942 if (k.gt.n_model_try) exit
2944 n_model_try=n_model_try+1
2945 list_model_try(n_model_try)=i_model
2946 if (me.eq.king .or. .not. out1file) &
2947 write (iout,*) 'Trying to start from model ',&
2948 pdbfiles_chomo(i_model)(:ilen(pdbfiles_chomo(i_model)))
2951 c(j,i)=chomo(j,i,i_model)
2954 call int_from_cart(.true.,.false.)
2955 call sc_loc_geom(.false.)
2959 dc(j,i)=c(j,i+1)-c(j,i)
2960 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
2965 dc(j,i+nres)=c(j,i+nres)-c(j,i)
2966 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
2969 if (me.eq.king.or..not.out1file) then
2970 write (iout,*) "Energies before removing overlaps"
2971 call etotal(energia(0))
2972 call enerprint(energia(0))
2974 ! Remove SC overlaps if requested
2976 write (iout,*) 'Calling OVERLAP_SC'
2977 call overlap_sc(fail)
2980 "Failed to remove overlap from model",i_model
2984 if (me.eq.king.or..not.out1file) then
2985 write (iout,*) "Energies after removing overlaps"
2986 call etotal(energia(0))
2987 call enerprint(energia(0))
2990 ! Search for better SC rotamers if requested
2992 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
2993 print *,'SC_move',nft_sc,etot
2994 if (me.eq.king.or..not.out1file)&
2995 write(iout,*) 'SC_move',nft_sc,etot
2997 call etotal(energia(0))
3000 call MPI_Gather(i_model,1,MPI_INTEGER,i_start_models(0),&
3001 1,MPI_INTEGER,king,CG_COMM,IERROR)
3002 if (n_model_try.gt.nmodel_start .and.&
3003 (me.eq.king .or. out1file)) then
3005 "All models have irreparable overlaps. Trying randoms starts."
3007 i_model=nmodel_start+1
3011 ! Remove SC overlaps if requested
3013 write (iout,*) 'Calling OVERLAP_SC'
3014 call overlap_sc(fail)
3017 "Failed to remove overlap"
3020 if (me.eq.king.or..not.out1file) then
3021 write (iout,*) "Energies after removing overlaps"
3022 call etotal(energia(0))
3023 call enerprint(energia(0))
3026 ! 8/22/17 AL Minimize initial structure
3028 if (me.eq.king.or..not.out1file) write(iout,*)&
3029 'Minimizing initial PDB structure: Calling MINIM_DC'
3030 call minim_dc(etot,iretcode,nfun)
3032 call geom_to_var(nvar,varia)
3033 if(me.eq.king.or..not.out1file) write (iout,*)&
3034 'Minimizing initial PDB structure: Calling MINIMIZE.'
3035 call minimize(etot,varia,iretcode,nfun)
3036 call var_to_geom(nvar,varia)
3038 if (me.eq.king.or..not.out1file)&
3039 write(iout,*) 'LBFGS return code is ',status,' eval ',nfun
3040 if(me.eq.king.or..not.out1file)&
3041 write(iout,*) 'LBFGS return code is ',status,' eval ',nfun
3043 if (me.eq.king.or..not.out1file)&
3044 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
3045 if(me.eq.king.or..not.out1file)&
3046 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
3050 if (nmodel_start.gt.0 .and. me.eq.king) then
3051 write (iout,'(a)') "Task Starting model"
3053 if (i_start_models(i).gt.nmodel_start) then
3054 write (iout,'(i4,2x,a)') i,"RANDOM STRUCTURE"
3056 write(iout,'(i4,2x,a)')i,pdbfiles_chomo(i_start_models(i)) &
3057 (:ilen(pdbfiles_chomo(i_start_models(i))))
3062 call chainbuild_cart
3064 write(iout,*) "just after kinetic"
3069 kinetic_T=2.0d0/(dimen3*Rb)*EK
3070 if(me.eq.king.or..not.out1file)then
3071 write(iout,*) "just after verlet_bath"
3081 write(iout,*) "before ETOTAL"
3082 call etotal(potEcomp)
3083 if (large) call enerprint(potEcomp)
3086 t_etotal=t_etotal+MPI_Wtime()-tt0
3088 t_etotal=t_etotal+tcpu()-tt0
3093 write(iout,*) "before lagrangian"
3095 write(iout,*) "before max_accel"
3097 if (amax*d_time .gt. dvmax) then
3098 d_time=d_time*dvmax/amax
3099 if(me.eq.king.or..not.out1file) write (iout,*) &
3100 "Time step reduced to",d_time,&
3101 " because of too large initial acceleration."
3103 if(me.eq.king.or..not.out1file)then
3104 write(iout,*) "Potential energy and its components"
3105 call enerprint(potEcomp)
3106 ! write(iout,*) (potEcomp(i),i=0,n_ene)
3108 potE=potEcomp(0)-potEcomp(51)
3112 if (ntwe.ne.0) call statout(itime)
3113 if(me.eq.king.or..not.out1file) &
3114 write (iout,'(/a/3(a25,1pe14.5/))') "Initial:", &
3115 " Kinetic energy",EK," Potential energy",potE, &
3116 " Total energy",totE," Maximum acceleration ", &
3119 write (iout,*) "Initial coordinates"
3121 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),&
3124 write (iout,*) "Initial dC"
3126 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),&
3127 (dc(j,i+nres),j=1,3)
3129 write (iout,*) "Initial velocities"
3130 write (iout,"(13x,' backbone ',23x,' side chain')")
3132 write (iout,'(i6,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
3133 (d_t(j,i+nres),j=1,3)
3135 write (iout,*) "Initial accelerations"
3137 ! write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
3138 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),&
3139 (d_a(j,i+nres),j=1,3)
3145 d_t_old(j,i)=d_t(j,i)
3146 d_a_old(j,i)=d_a(j,i)
3148 ! write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
3157 call etotal_short(energia_short)
3158 if (large) call enerprint(potEcomp)
3161 t_eshort=t_eshort+MPI_Wtime()-tt0
3163 t_eshort=t_eshort+tcpu()-tt0
3168 if(.not.out1file .and. large) then
3169 write (iout,*) "energia_long",energia_long(0),&
3170 " energia_short",energia_short(0),&
3171 " total",energia_long(0)+energia_short(0)
3172 write (iout,*) "Initial fast-force accelerations"
3174 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
3175 (d_a(j,i+nres),j=1,3)
3178 ! 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
3181 d_a_short(j,i)=d_a(j,i)
3190 call etotal_long(energia_long)
3191 if (large) call enerprint(potEcomp)
3194 t_elong=t_elong+MPI_Wtime()-tt0
3196 t_elong=t_elong+tcpu()-tt0
3201 if(.not.out1file .and. large) then
3202 write (iout,*) "energia_long",energia_long(0)
3203 write (iout,*) "Initial slow-force accelerations"
3205 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
3206 (d_a(j,i+nres),j=1,3)
3210 t_enegrad=t_enegrad+MPI_Wtime()-tt0
3212 t_enegrad=t_enegrad+tcpu()-tt0
3216 end subroutine init_MD
3217 !-----------------------------------------------------------------------------
3218 subroutine random_vel
3220 ! implicit real*8 (a-h,o-z)
3222 use random, only:anorm_distr
3224 ! include 'DIMENSIONS'
3225 ! include 'COMMON.CONTROL'
3226 ! include 'COMMON.VAR'
3227 ! include 'COMMON.MD'
3229 ! include 'COMMON.LANGEVIN'
3231 ! include 'COMMON.LANGEVIN.lang0'
3233 ! include 'COMMON.CHAIN'
3234 ! include 'COMMON.DERIV'
3235 ! include 'COMMON.GEO'
3236 ! include 'COMMON.LOCAL'
3237 ! include 'COMMON.INTERACT'
3238 ! include 'COMMON.IOUNITS'
3239 ! include 'COMMON.NAMES'
3240 ! include 'COMMON.TIME1'
3241 real(kind=8) :: xv,sigv,lowb,highb ,Ek1
3243 integer ichain,n,innt,inct,ibeg,ierr
3244 real(kind=8) ,allocatable, dimension(:):: work
3245 integer,allocatable,dimension(:) :: iwork
3246 ! double precision Ghalf(mmaxres2_chain),Geigen(maxres2_chain),&
3247 ! Gvec(maxres2_chain,maxres2_chain)
3248 ! common /przechowalnia/Ghalf,Geigen,Gvec
3250 ! double precision inertia(maxres2_chain,maxres2_chain)
3255 real(kind=8) ,allocatable, dimension(:) :: xsolv,DML,rs
3256 real(kind=8) :: sumx,Ek2,Ek3,aux
3258 real(kind=8) ,allocatable, dimension(:) :: rsold
3259 real (kind=8),allocatable,dimension(:,:) :: matold,inertia
3263 integer :: i,j,ii,k,mark,imark,mnum,nres2
3264 integer(kind=8) :: ind
3265 ! Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
3267 ! First generate velocities in the eigenspace of the G matrix
3268 ! write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
3271 if(.not.allocated(work)) then
3272 allocate(work(48*nres))
3273 allocate(iwork(6*nres))
3275 print *,"IN RANDOM VEL"
3277 ! print *,size(ghalf)
3280 write (iout,*) "Random_vel, fivediag"
3282 allocate(inertia(2*nres,2*nres))
3289 write(iout,*), nchain
3293 if(.not.allocated(ghalf)) print *,"COCO"
3294 if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2))
3296 n=dimen_chain(ichain)
3297 innt=iposd_chain(ichain)
3300 write (iout,*) "Chain",ichain," n",n," start",innt
3302 if (i.lt.inct-1) then
3303 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i),DU1orig(i),&
3305 else if (i.eq.inct-1) then
3306 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i),DU1orig(i)
3308 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i)
3313 ghalf(ind+1)=dmorig(innt)
3314 ghalf(ind+2)=du1orig(innt)
3315 ghalf(ind+3)=dmorig(innt+1)
3319 write (iout,*) "i",i," ind",ind," indu2",innt+i-2,&
3320 " indu1",innt+i-1," indm",innt+i
3321 ghalf(ind+1)=du2orig(innt-1+i-2)
3322 ghalf(ind+2)=du1orig(innt-1+i-1)
3323 ghalf(ind+3)=dmorig(innt-1+i)
3324 !c write (iout,'(3(a,i2,1x))') "DU2",innt-1+i-2,
3325 !c & "DU1",innt-1+i-1,"DM ",innt-1+i
3333 inertia(i,j)=ghalf(ind)
3334 inertia(j,i)=ghalf(ind)
3339 write (iout,*) "Chain ",ichain," ind",ind," dim",n*(n+1)/2
3340 write (iout,*) "Five-diagonal inertia matrix, lower triangle"
3341 ! call matoutr(n,ghalf)
3343 call gldiag(nres*2,n,n,Ghalf,work,Geigen,Gvec,ierr,iwork)
3345 write (iout,'(//a,i3)')&
3346 "Eigenvectors and eigenvalues of the G matrix chain",ichain
3347 call eigout(n,n,nres*2,nres*2,Gvec,Geigen)
3350 !c check diagonalization
3356 aux=aux+gvec(k,i)*gvec(l,j)*inertia(k,l)
3360 write (iout,*) i,j,aux,geigen(i)
3362 write (iout,*) i,j,aux
3372 sigv=dsqrt((Rb*t_bath)/geigen(i))
3375 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
3376 EK=EK+0.5d0*geigen(i)*d_t_work_new(ii)**2
3377 !c write (iout,*) "i",i," ii",ii," geigen",geigen(i),
3378 !c & " d_t_work_new",d_t_work_new(ii)
3386 d_t_work(ind)=d_t_work(ind)&
3387 +Gvec(i,j)*d_t_work_new((j-1)*3+k)
3389 !c write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
3398 aux=aux+inertia(i,j)*d_t_work(3*(i-1)+k)*d_t_work(3*(j-1)+k)
3404 !c Transfer to the d_t vector
3405 innt=chain_border(1,ichain)
3406 inct=chain_border(2,ichain)
3408 !c write (iout,*) "ichain",ichain," innt",innt," inct",inct
3412 d_t(j,i)=d_t_work(ind)
3415 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum).and.mnum.le.2) then
3418 d_t(j,i+nres)=d_t_work(ind)
3425 write (iout,*) "Random velocities in the Calpha,SC space"
3428 write (iout,'(a3,1h(,i5,1h),3f10.5,3x,3f10.5)')&
3429 restyp(itype(i,mnum),mnum),i,(d_t(j,i),j=1,3),(d_t(j,i+nres),j=1,3)
3432 call kinetic_CASC(Ek1)
3434 ! Transform the velocities to virtual-bond space
3443 if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum).or.mnum.ge.3) then
3445 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
3449 d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
3450 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
3461 !c d_a(:,0)=d_a(:,1)
3463 !c write (iout,*) "Shifting accelerations"
3465 !c write (iout,*) "ichain",chain_border1(1,ichain)-1,
3466 !c & chain_border1(1,ichain)
3467 d_t(:,chain_border1(1,ichain)-1)=d_t(:,chain_border1(1,ichain))
3468 d_t(:,chain_border1(1,ichain))=0.0d0
3470 !c write (iout,*) "Adding accelerations"
3472 !c write (iout,*) "chain",ichain,chain_border1(1,ichain)-1,
3473 !c & chain_border(2,ichain-1)
3474 d_t(:,chain_border1(1,ichain)-1)=&
3475 d_t(:,chain_border1(1,ichain)-1)+d_t(:,chain_border(2,ichain-1))
3476 d_t(:,chain_border(2,ichain-1))=0.0d0
3479 write (iout,*) "chain",ichain,chain_border1(1,ichain)-1,&
3480 chain_border(2,ichain-1)
3481 d_t(:,chain_border1(1,ichain)-1)=&
3482 d_t(:,chain_border1(1,ichain)-1)+d_t(:,chain_border(2,ichain-1))
3483 d_t(:,chain_border(2,ichain-1))=0.0d0
3488 !c d_t(j,0)=d_t(j,nnt)
3491 innt=chain_border(1,ichain)
3492 inct=chain_border(2,ichain)
3493 !c write (iout,*) "ichain",ichain," innt",innt," inct",inct
3494 !c write (iout,*) "ibeg",ibeg
3496 d_t(j,ibeg)=d_t(j,innt)
3501 if (iabs(itype(i,1).eq.10).or.mnum.ge.3) then
3502 !c write (iout,*) "i",i,(d_t(j,i),j=1,3),(d_t(j,i+1),j=1,3)
3504 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
3508 d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
3509 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
3518 "Random velocities in the virtual-bond-vector space"
3519 write (iout,'(3hORG,1h(,i5,1h),3f10.5)') 0,(d_t(j,0),j=1,3)
3521 write (iout,'(a3,1h(,i5,1h),3f10.5,3x,3f10.5)')&
3522 restyp(itype(i,mnum),mnum),i,(d_t(j,i),j=1,3),(d_t(j,i+nres),j=1,3)
3525 write (iout,*) "Kinetic energy from inertia matrix eigenvalues",&
3528 "Kinetic temperatures from inertia matrix eigenvalues",&
3531 write (iout,*) "Kinetic energy from inertia matrix",Ek3
3532 write (iout,*) "Kinetic temperatures from inertia",&
3535 write (iout,*) "Kinetic energy from velocities in CA-SC space",&
3538 "Kinetic temperatures from velovities in CA-SC space",&
3542 "Kinetic energy from virtual-bond-vector velocities",Ek1
3544 "Kinetic temperature from virtual-bond-vector velocities ",&
3553 sigv=dsqrt((Rb*t_bath)/geigen(i))
3556 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
3558 write (iout,*) "i",i," ii",ii," geigen",geigen(i),&
3559 " d_t_work_new",d_t_work_new(ii)
3570 Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
3573 write (iout,*) "Ek from eigenvectors",Ek1
3574 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
3583 d_t_work(ind)=d_t_work(ind) &
3584 +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
3586 ! write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
3590 ! Transfer to the d_t vector
3592 d_t(j,0)=d_t_work(j)
3598 d_t(j,i)=d_t_work(ind)
3603 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3604 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3605 .and.(mnum.lt.4)) then
3608 d_t(j,i+nres)=d_t_work(ind)
3614 ! write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",&
3615 ! 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
3617 ! write(iout,*) "end init MD"
3620 end subroutine random_vel
3621 !-----------------------------------------------------------------------------
3623 subroutine sd_verlet_p_setup
3624 ! Sets up the parameters of stochastic Verlet algorithm
3625 ! implicit real*8 (a-h,o-z)
3626 ! include 'DIMENSIONS'
3627 use control, only: tcpu
3632 ! include 'COMMON.CONTROL'
3633 ! include 'COMMON.VAR'
3634 ! include 'COMMON.MD'
3636 ! include 'COMMON.LANGEVIN'
3638 ! include 'COMMON.LANGEVIN.lang0'
3640 ! include 'COMMON.CHAIN'
3641 ! include 'COMMON.DERIV'
3642 ! include 'COMMON.GEO'
3643 ! include 'COMMON.LOCAL'
3644 ! include 'COMMON.INTERACT'
3645 ! include 'COMMON.IOUNITS'
3646 ! include 'COMMON.NAMES'
3647 ! include 'COMMON.TIME1'
3648 real(kind=8),dimension(6*nres) :: emgdt !(MAXRES6) maxres6=6*maxres
3649 real(kind=8) :: pterm,vterm,rho,rhoc,vsig
3650 real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
3651 prand_vec,vrand_vec1,vrand_vec2 !(MAXRES6) maxres6=6*maxres
3652 logical :: lprn = .false.
3653 real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
3654 real(kind=8) :: ktm,gdt,egdt,gdt2,gdt3,gdt4,gdt5,gdt6,gdt7,gdt8,&
3656 integer :: i,maxres2
3663 ! AL 8/17/04 Code adapted from tinker
3665 ! Get the frictional and random terms for stochastic dynamics in the
3666 ! eigenspace of mass-scaled UNRES friction matrix
3670 gdt = fricgam(i) * d_time
3672 ! Stochastic dynamics reduces to simple MD for zero friction
3674 if (gdt .le. zero) then
3675 pfric_vec(i) = 1.0d0
3676 vfric_vec(i) = d_time
3677 afric_vec(i) = 0.5d0 * d_time * d_time
3678 prand_vec(i) = 0.0d0
3679 vrand_vec1(i) = 0.0d0
3680 vrand_vec2(i) = 0.0d0
3682 ! Analytical expressions when friction coefficient is large
3685 if (gdt .ge. gdt_radius) then
3688 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
3689 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
3690 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
3691 vterm = 1.0d0 - egdt**2
3692 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
3694 ! Use series expansions when friction coefficient is small
3705 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0 &
3706 - gdt5/120.0d0 + gdt6/720.0d0 &
3707 - gdt7/5040.0d0 + gdt8/40320.0d0 &
3708 - gdt9/362880.0d0) / fricgam(i)**2
3709 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
3710 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
3711 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0 &
3712 + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0 &
3713 + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0 &
3714 + 127.0d0*gdt9/90720.0d0
3715 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0 &
3716 - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0 &
3717 - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0 &
3718 - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
3719 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0 &
3720 - 17.0d0*gdt2/1280.0d0 &
3721 + 17.0d0*gdt3/6144.0d0 &
3722 + 40967.0d0*gdt4/34406400.0d0 &
3723 - 57203.0d0*gdt5/275251200.0d0 &
3724 - 1429487.0d0*gdt6/13212057600.0d0)
3727 ! Compute the scaling factors of random terms for the nonzero friction case
3729 ktm = 0.5d0*d_time/fricgam(i)
3730 psig = dsqrt(ktm*pterm) / fricgam(i)
3731 vsig = dsqrt(ktm*vterm)
3732 rhoc = dsqrt(1.0d0 - rho*rho)
3734 vrand_vec1(i) = vsig * rho
3735 vrand_vec2(i) = vsig * rhoc
3740 "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",&
3743 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),&
3744 afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
3748 ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
3751 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
3752 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
3753 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
3754 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
3755 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
3756 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
3759 t_sdsetup=t_sdsetup+MPI_Wtime()
3761 t_sdsetup=t_sdsetup+tcpu()-tt0
3764 end subroutine sd_verlet_p_setup
3765 !-----------------------------------------------------------------------------
3766 subroutine eigtransf1(n,ndim,ab,d,c)
3770 real(kind=8) :: ab(ndim,ndim,n),c(ndim,n),d(ndim)
3776 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
3781 end subroutine eigtransf1
3782 !-----------------------------------------------------------------------------
3783 subroutine eigtransf(n,ndim,a,b,d,c)
3787 real(kind=8) :: a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
3793 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
3798 end subroutine eigtransf
3799 !-----------------------------------------------------------------------------
3800 subroutine sd_verlet1
3802 ! Applying stochastic velocity Verlet algorithm - step 1 to velocities
3804 ! implicit real*8 (a-h,o-z)
3805 ! include 'DIMENSIONS'
3806 ! include 'COMMON.CONTROL'
3807 ! include 'COMMON.VAR'
3808 ! include 'COMMON.MD'
3810 ! include 'COMMON.LANGEVIN'
3812 ! include 'COMMON.LANGEVIN.lang0'
3814 ! include 'COMMON.CHAIN'
3815 ! include 'COMMON.DERIV'
3816 ! include 'COMMON.GEO'
3817 ! include 'COMMON.LOCAL'
3818 ! include 'COMMON.INTERACT'
3819 ! include 'COMMON.IOUNITS'
3820 ! include 'COMMON.NAMES'
3821 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
3822 !el common /stochcalc/ stochforcvec
3823 logical :: lprn = .false.
3824 real(kind=8) :: ddt1,ddt2
3825 integer :: i,j,ind,inres
3827 ! write (iout,*) "dc_old"
3829 ! write (iout,'(i5,3f10.5,5x,3f10.5)')
3830 ! & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3833 dc_work(j)=dc_old(j,0)
3834 d_t_work(j)=d_t_old(j,0)
3835 d_a_work(j)=d_a_old(j,0)
3840 dc_work(ind+j)=dc_old(j,i)
3841 d_t_work(ind+j)=d_t_old(j,i)
3842 d_a_work(ind+j)=d_a_old(j,i)
3848 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3849 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3850 .and.(mnum.lt.4)) then
3852 dc_work(ind+j)=dc_old(j,i+nres)
3853 d_t_work(ind+j)=d_t_old(j,i+nres)
3854 d_a_work(ind+j)=d_a_old(j,i+nres)
3862 "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",&
3866 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
3867 vfric_mat(i,j),afric_mat(i,j),&
3868 prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3876 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) &
3877 +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3878 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3879 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3881 d_t_work_new(i)=ddt1+0.5d0*ddt2
3882 d_t_work(i)=ddt1+ddt2
3887 d_t(j,0)=d_t_work(j)
3892 dc(j,i)=dc_work(ind+j)
3893 d_t(j,i)=d_t_work(ind+j)
3899 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3900 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3901 .and.(mnum.lt.4)) then
3904 dc(j,inres)=dc_work(ind+j)
3905 d_t(j,inres)=d_t_work(ind+j)
3911 end subroutine sd_verlet1
3912 !-----------------------------------------------------------------------------
3913 subroutine sd_verlet2
3915 ! Calculating the adjusted velocities for accelerations
3917 ! implicit real*8 (a-h,o-z)
3918 ! include 'DIMENSIONS'
3919 ! include 'COMMON.CONTROL'
3920 ! include 'COMMON.VAR'
3921 ! include 'COMMON.MD'
3923 ! include 'COMMON.LANGEVIN'
3925 ! include 'COMMON.LANGEVIN.lang0'
3927 ! include 'COMMON.CHAIN'
3928 ! include 'COMMON.DERIV'
3929 ! include 'COMMON.GEO'
3930 ! include 'COMMON.LOCAL'
3931 ! include 'COMMON.INTERACT'
3932 ! include 'COMMON.IOUNITS'
3933 ! include 'COMMON.NAMES'
3934 !el real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV !(MAXRES6) maxres6=6*maxres
3935 real(kind=8),dimension(6*nres) :: stochforcvecV !(MAXRES6) maxres6=6*maxres
3936 !el common /stochcalc/ stochforcvec
3938 real(kind=8) :: ddt1,ddt2
3939 integer :: i,j,ind,inres
3940 ! Compute the stochastic forces which contribute to velocity change
3942 call stochastic_force(stochforcvecV)
3949 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3950 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+ &
3951 vrand_mat2(i,j)*stochforcvecV(j)
3953 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
3957 d_t(j,0)=d_t_work(j)
3962 d_t(j,i)=d_t_work(ind+j)
3968 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3969 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3970 .and.(mnum.lt.4)) then
3973 d_t(j,inres)=d_t_work(ind+j)
3979 end subroutine sd_verlet2
3980 !-----------------------------------------------------------------------------
3981 subroutine sd_verlet_ciccotti_setup
3983 ! Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
3985 ! implicit real*8 (a-h,o-z)
3986 ! include 'DIMENSIONS'
3987 use control, only: tcpu
3992 ! include 'COMMON.CONTROL'
3993 ! include 'COMMON.VAR'
3994 ! include 'COMMON.MD'
3996 ! include 'COMMON.LANGEVIN'
3998 ! include 'COMMON.LANGEVIN.lang0'
4000 ! include 'COMMON.CHAIN'
4001 ! include 'COMMON.DERIV'
4002 ! include 'COMMON.GEO'
4003 ! include 'COMMON.LOCAL'
4004 ! include 'COMMON.INTERACT'
4005 ! include 'COMMON.IOUNITS'
4006 ! include 'COMMON.NAMES'
4007 ! include 'COMMON.TIME1'
4008 real(kind=8),dimension(6*nres) :: emgdt !(MAXRES6) maxres6=6*maxres
4009 real(kind=8) :: pterm,vterm,rho,rhoc,vsig
4010 real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
4011 prand_vec,vrand_vec1,vrand_vec2 !(MAXRES6) maxres6=6*maxres
4012 logical :: lprn = .false.
4013 real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
4014 real(kind=8) :: ktm,gdt,egdt,tt0
4015 integer :: i,maxres2
4022 ! AL 8/17/04 Code adapted from tinker
4024 ! Get the frictional and random terms for stochastic dynamics in the
4025 ! eigenspace of mass-scaled UNRES friction matrix
4029 write (iout,*) "i",i," fricgam",fricgam(i)
4030 gdt = fricgam(i) * d_time
4032 ! Stochastic dynamics reduces to simple MD for zero friction
4034 if (gdt .le. zero) then
4035 pfric_vec(i) = 1.0d0
4036 vfric_vec(i) = d_time
4037 afric_vec(i) = 0.5d0*d_time*d_time
4038 prand_vec(i) = afric_vec(i)
4039 vrand_vec2(i) = vfric_vec(i)
4041 ! Analytical expressions when friction coefficient is large
4046 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
4047 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
4048 prand_vec(i) = afric_vec(i)
4049 vrand_vec2(i) = vfric_vec(i)
4051 ! Compute the scaling factors of random terms for the nonzero friction case
4053 ! ktm = 0.5d0*d_time/fricgam(i)
4054 ! psig = dsqrt(ktm*pterm) / fricgam(i)
4055 ! vsig = dsqrt(ktm*vterm)
4056 ! prand_vec(i) = psig*afric_vec(i)
4057 ! vrand_vec2(i) = vsig*vfric_vec(i)
4062 "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",&
4065 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),&
4066 afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
4070 ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
4072 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
4073 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
4074 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
4075 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
4076 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
4078 t_sdsetup=t_sdsetup+MPI_Wtime()
4080 t_sdsetup=t_sdsetup+tcpu()-tt0
4083 end subroutine sd_verlet_ciccotti_setup
4084 !-----------------------------------------------------------------------------
4085 subroutine sd_verlet1_ciccotti
4087 ! Applying stochastic velocity Verlet algorithm - step 1 to velocities
4088 ! implicit real*8 (a-h,o-z)
4090 ! include 'DIMENSIONS'
4094 ! include 'COMMON.CONTROL'
4095 ! include 'COMMON.VAR'
4096 ! include 'COMMON.MD'
4098 ! include 'COMMON.LANGEVIN'
4100 ! include 'COMMON.LANGEVIN.lang0'
4102 ! include 'COMMON.CHAIN'
4103 ! include 'COMMON.DERIV'
4104 ! include 'COMMON.GEO'
4105 ! include 'COMMON.LOCAL'
4106 ! include 'COMMON.INTERACT'
4107 ! include 'COMMON.IOUNITS'
4108 ! include 'COMMON.NAMES'
4109 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
4110 !el common /stochcalc/ stochforcvec
4111 logical :: lprn = .false.
4112 real(kind=8) :: ddt1,ddt2
4113 integer :: i,j,ind,inres
4114 ! write (iout,*) "dc_old"
4116 ! write (iout,'(i5,3f10.5,5x,3f10.5)')
4117 ! & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
4120 dc_work(j)=dc_old(j,0)
4121 d_t_work(j)=d_t_old(j,0)
4122 d_a_work(j)=d_a_old(j,0)
4127 dc_work(ind+j)=dc_old(j,i)
4128 d_t_work(ind+j)=d_t_old(j,i)
4129 d_a_work(ind+j)=d_a_old(j,i)
4134 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
4136 dc_work(ind+j)=dc_old(j,i+nres)
4137 d_t_work(ind+j)=d_t_old(j,i+nres)
4138 d_a_work(ind+j)=d_a_old(j,i+nres)
4147 "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",&
4151 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
4152 vfric_mat(i,j),afric_mat(i,j),&
4153 prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
4161 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) &
4162 +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
4163 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
4164 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
4166 d_t_work_new(i)=ddt1+0.5d0*ddt2
4167 d_t_work(i)=ddt1+ddt2
4172 d_t(j,0)=d_t_work(j)
4177 dc(j,i)=dc_work(ind+j)
4178 d_t(j,i)=d_t_work(ind+j)
4184 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
4185 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4186 .and.(mnum.lt.4)) then
4189 dc(j,inres)=dc_work(ind+j)
4190 d_t(j,inres)=d_t_work(ind+j)
4196 end subroutine sd_verlet1_ciccotti
4197 !-----------------------------------------------------------------------------
4198 subroutine sd_verlet2_ciccotti
4200 ! Calculating the adjusted velocities for accelerations
4202 ! implicit real*8 (a-h,o-z)
4203 ! include 'DIMENSIONS'
4204 ! include 'COMMON.CONTROL'
4205 ! include 'COMMON.VAR'
4206 ! include 'COMMON.MD'
4208 ! include 'COMMON.LANGEVIN'
4210 ! include 'COMMON.LANGEVIN.lang0'
4212 ! include 'COMMON.CHAIN'
4213 ! include 'COMMON.DERIV'
4214 ! include 'COMMON.GEO'
4215 ! include 'COMMON.LOCAL'
4216 ! include 'COMMON.INTERACT'
4217 ! include 'COMMON.IOUNITS'
4218 ! include 'COMMON.NAMES'
4219 !el real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV !(MAXRES6) maxres6=6*maxres
4220 real(kind=8),dimension(6*nres) :: stochforcvecV !(MAXRES6) maxres6=6*maxres
4221 !el common /stochcalc/ stochforcvec
4222 real(kind=8) :: ddt1,ddt2
4223 integer :: i,j,ind,inres
4225 ! Compute the stochastic forces which contribute to velocity change
4227 call stochastic_force(stochforcvecV)
4234 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
4235 ! ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
4236 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
4238 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
4242 d_t(j,0)=d_t_work(j)
4247 d_t(j,i)=d_t_work(ind+j)
4253 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4255 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
4258 d_t(j,inres)=d_t_work(ind+j)
4264 end subroutine sd_verlet2_ciccotti
4266 !-----------------------------------------------------------------------------
4268 !-----------------------------------------------------------------------------
4270 subroutine inertia_tensor
4273 real(kind=8) Im(3,3),Imcp(3,3),pr(3),M_SC,&
4274 eigvec(3,3),Id(3,3),eigval(3),L(3),vp(3),vrot(3),&
4275 vpp(3,0:MAXRES),vs_p(3),pr1(3,3),&
4276 pr2(3,3),pp(3),incr(3),v(3),mag,mag2,M_PEP
4277 integer iti,inres,i,j,k,mnum,mnum1
4290 !c caulating the center of the mass of the protein
4294 if (itype(i,mnum).eq.ntyp1_molec(mnum)&
4295 .or. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) cycle
4296 ! if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
4297 if (mnum.ge.5) mp(mnum)=0.0d0
4298 M_PEP=M_PEP+mp(mnum)
4301 cm(j)=cm(j)+(c(j,i)+0.5d0*dc(j,i))*mp(mnum)
4311 iti=iabs(itype(i,mnum))
4312 if (iti.eq.ntyp1_molec(mnum)) cycle
4313 M_SC=M_SC+msc(iabs(iti),mnum)
4316 cm(j)=cm(j)+msc(iabs(iti),mnum)*c(j,inres)
4320 cm(j)=cm(j)/(M_SC+M_PEP)
4325 ! if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
4326 if (mnum.ge.5) mp(mnum)=0.0d0
4327 if (itype(i,mnum).eq.ntyp1_molec(mnum)&
4328 .or. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) cycle
4330 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
4332 Im(1,1)=Im(1,1)+mp(mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
4333 Im(1,2)=Im(1,2)-mp(mnum)*pr(1)*pr(2)
4334 Im(1,3)=Im(1,3)-mp(mnum)*pr(1)*pr(3)
4335 Im(2,3)=Im(2,3)-mp(mnum)*pr(2)*pr(3)
4336 Im(2,2)=Im(2,2)+mp(mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
4337 Im(3,3)=Im(3,3)+mp(mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
4342 iti=iabs(itype(i,mnum))
4343 if (iti.eq.ntyp1_molec(mnum)) cycle
4346 pr(j)=c(j,inres)-cm(j)
4348 Im(1,1)=Im(1,1)+msc(iabs(iti),mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
4349 Im(1,2)=Im(1,2)-msc(iabs(iti),mnum)*pr(1)*pr(2)
4350 Im(1,3)=Im(1,3)-msc(iabs(iti),mnum)*pr(1)*pr(3)
4351 Im(2,3)=Im(2,3)-msc(iabs(iti),mnum)*pr(2)*pr(3)
4352 Im(2,2)=Im(2,2)+msc(iabs(iti),mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
4353 Im(3,3)=Im(3,3)+msc(iabs(iti),mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
4358 if (itype(i,mnum).eq.ntyp1_molec(mnum)&
4359 .or. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) cycle
4360 Im(1,1)=Im(1,1)+Ip(mnum)*(1-dc_norm(1,i)*dc_norm(1,i))*&
4361 vbld(i+1)*vbld(i+1)*0.25d0
4362 Im(1,2)=Im(1,2)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(2,i))*&
4363 vbld(i+1)*vbld(i+1)*0.25d0
4364 Im(1,3)=Im(1,3)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(3,i))*&
4365 vbld(i+1)*vbld(i+1)*0.25d0
4366 Im(2,3)=Im(2,3)+Ip(mnum)*(-dc_norm(2,i)*dc_norm(3,i))*&
4367 vbld(i+1)*vbld(i+1)*0.25d0
4368 Im(2,2)=Im(2,2)+Ip(mnum)*(1-dc_norm(2,i)*dc_norm(2,i))*&
4369 vbld(i+1)*vbld(i+1)*0.25d0
4370 Im(3,3)=Im(3,3)+Ip(mnum)*(1-dc_norm(3,i)*dc_norm(3,i))*&
4371 vbld(i+1)*vbld(i+1)*0.25d0
4376 iti=iabs(itype(i,mnum))
4377 if (iti.ne.10 .and. iti.ne.ntyp1_molec(mnum).and.mnum.le.2) then
4379 Im(1,1)=Im(1,1)+Isc(iti,mnum)*(1-dc_norm(1,inres)*&
4380 dc_norm(1,inres))*vbld(inres)*vbld(inres)
4381 Im(1,2)=Im(1,2)-Isc(iti,mnum)*(dc_norm(1,inres)*&
4382 dc_norm(2,inres))*vbld(inres)*vbld(inres)
4383 Im(1,3)=Im(1,3)-Isc(iti,mnum)*(dc_norm(1,inres)*&
4384 dc_norm(3,inres))*vbld(inres)*vbld(inres)
4385 Im(2,3)=Im(2,3)-Isc(iti,mnum)*(dc_norm(2,inres)*&
4386 dc_norm(3,inres))*vbld(inres)*vbld(inres)
4387 Im(2,2)=Im(2,2)+Isc(iti,mnum)*(1-dc_norm(2,inres)*&
4388 dc_norm(2,inres))*vbld(inres)*vbld(inres)
4389 Im(3,3)=Im(3,3)+Isc(iti,mnum)*(1-dc_norm(3,inres)*&
4390 dc_norm(3,inres))*vbld(inres)*vbld(inres)
4399 !c Copng the Im matrix for the djacob subroutine
4406 !c Finding the eigenvectors and eignvalues of the inertia tensor
4407 call djacob(3,3,10000,1.0d-10,Imcp,eigvec,eigval)
4409 if (dabs(eigval(i)).gt.1.0d-15) then
4410 Id(i,i)=1.0d0/eigval(i)
4417 Imcp(i,j)=eigvec(j,i)
4423 pr1(i,j)=pr1(i,j)+Id(i,k)*Imcp(k,j)
4430 pr2(i,j)=pr2(i,j)+eigvec(i,k)*pr1(k,j)
4434 !c Calculating the total rotational velocity of the molecule
4437 vrot(i)=vrot(i)+pr2(i,j)*L(j)
4443 ! write (iout,*) itype(i+1,mnum1),itype(i,mnum)
4444 if (itype(i+1,mnum1).ne.ntyp1_molec(mnum1) &
4445 .and. itype(i,mnum).eq.ntyp1_molec(mnum) .or.&
4446 itype(i,mnum).ne.ntyp1_molec(mnum) &
4447 .and. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) cycle
4448 call vecpr(vrot(1),dc(1,i),vp)
4450 d_t(j,i)=d_t(j,i)-vp(j)
4455 if(itype(i,1).ne.10 .and.itype(i,mnum).ne.ntyp1_molec(mnum)&
4456 .and.mnum.le.2) then
4458 call vecpr(vrot(1),dc(1,inres),vp)
4460 d_t(j,inres)=d_t(j,inres)-vp(j)
4466 end subroutine inertia_tensor
4467 !------------------------------------------------------------
4468 subroutine angmom(cm,L)
4470 double precision L(3),cm(3),pr(3),vp(3),vrot(3),incr(3),v(3),&
4472 integer iti,inres,i,j,mnum,mnum1
4473 !c Calculate the angular momentum
4483 ! if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
4484 if (mnum.ge.5) mp(mnum)=0.0d0
4485 if (itype(i,mnum).eq.ntyp1_molec(mnum)&
4486 .or. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) cycle
4488 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
4491 v(j)=incr(j)+0.5d0*d_t(j,i)
4494 incr(j)=incr(j)+d_t(j,i)
4496 call vecpr(pr(1),v(1),vp)
4498 L(j)=L(j)+mp(mnum)*vp(j)
4502 pp(j)=0.5d0*d_t(j,i)
4504 call vecpr(pr(1),pp(1),vp)
4506 L(j)=L(j)+Ip(mnum)*vp(j)
4514 iti=iabs(itype(i,mnum))
4517 pr(j)=c(j,inres)-cm(j)
4519 if (itype(i,1).ne.10 .and.itype(i,mnum).ne.ntyp1_molec(mnum)&
4520 .and.mnum.le.2) then
4522 v(j)=incr(j)+d_t(j,inres)
4529 call vecpr(pr(1),v(1),vp)
4530 !c write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),
4531 !c & " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
4532 ! if (mnum.gt.4) then
4538 L(j)=L(j)+mscab*vp(j)
4540 !c write (iout,*) "L",(l(j),j=1,3)
4541 if (itype(i,1).ne.10 .and.itype(i,mnum).ne.ntyp1_molec(mnum)&
4542 .and.mnum.le.2) then
4544 v(j)=incr(j)+d_t(j,inres)
4546 call vecpr(dc(1,inres),d_t(1,inres),vp)
4548 L(j)=L(j)+Isc(iti,mnum)*vp(j)
4552 incr(j)=incr(j)+d_t(j,i)
4556 end subroutine angmom
4557 !---------------------------------------------------------------
4558 subroutine vcm_vel(vcm)
4559 double precision vcm(3),vv(3),summas,amas
4560 integer i,j,mnum,mnum1
4568 if ((mnum.ge.5).or.(mnum.eq.3))&
4569 mp(mnum)=msc(itype(i,mnum),mnum)
4571 summas=summas+mp(mnum)
4573 vcm(j)=vcm(j)+mp(mnum)*(vv(j)+0.5d0*d_t(j,i))
4577 amas=msc(iabs(itype(i,mnum)),mnum)
4581 ! amas=msc(iabs(itype(i)))
4583 ! if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
4584 if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4585 .and.(mnum.lt.4)) then
4587 vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres))
4591 vcm(j)=vcm(j)+amas*vv(j)
4595 vv(j)=vv(j)+d_t(j,i)
4598 !c write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas
4600 vcm(j)=vcm(j)/summas
4603 end subroutine vcm_vel
4605 subroutine inertia_tensor
4607 ! Calculating the intertia tensor for the entire protein in order to
4608 ! remove the perpendicular components of velocity matrix which cause
4609 ! the molecule to rotate.
4612 ! implicit real*8 (a-h,o-z)
4613 ! include 'DIMENSIONS'
4614 ! include 'COMMON.CONTROL'
4615 ! include 'COMMON.VAR'
4616 ! include 'COMMON.MD'
4617 ! include 'COMMON.CHAIN'
4618 ! include 'COMMON.DERIV'
4619 ! include 'COMMON.GEO'
4620 ! include 'COMMON.LOCAL'
4621 ! include 'COMMON.INTERACT'
4622 ! include 'COMMON.IOUNITS'
4623 ! include 'COMMON.NAMES'
4625 real(kind=8),dimension(3,3) :: Im,Imcp,eigvec,Id
4626 real(kind=8),dimension(3) :: pr,eigval,L,vp,vrot
4627 real(kind=8) :: M_SC,mag,mag2,M_PEP
4628 real(kind=8),dimension(3,0:nres) :: vpp !(3,0:MAXRES)
4629 real(kind=8),dimension(3) :: vs_p,pp,incr,v
4630 real(kind=8),dimension(3,3) :: pr1,pr2
4632 !el common /gucio/ cm
4633 integer :: iti,inres,i,j,k,mnum
4644 ! calculating the center of the mass of the protein
4648 if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
4649 write(iout,*) "WTF",itype(i,mnum),i,mnum,mp(mnum)
4650 ! if (itype(i,mnum).eq.ntyp1_molec(mnum)) cycle
4651 M_PEP=M_PEP+mp(mnum)
4654 cm(j)=cm(j)+(c(j,i)+0.5d0*dc(j,i))*mp(mnum)
4663 if (mnum.ge.5) cycle
4664 iti=iabs(itype(i,mnum))
4665 M_SC=M_SC+msc(iabs(iti),mnum)
4668 cm(j)=cm(j)+msc(iabs(iti),mnum)*c(j,inres)
4672 cm(j)=cm(j)/(M_SC+M_PEP)
4674 ! write(iout,*) "Center of mass:",cm
4677 if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
4679 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
4681 Im(1,1)=Im(1,1)+mp(mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
4682 Im(1,2)=Im(1,2)-mp(mnum)*pr(1)*pr(2)
4683 Im(1,3)=Im(1,3)-mp(mnum)*pr(1)*pr(3)
4684 Im(2,3)=Im(2,3)-mp(mnum)*pr(2)*pr(3)
4685 Im(2,2)=Im(2,2)+mp(mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
4686 Im(3,3)=Im(3,3)+mp(mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
4689 ! write(iout,*) "The angular momentum before msc add"
4691 ! write (iout,*) (Im(i,j),j=1,3)
4696 iti=iabs(itype(i,mnum))
4697 if (mnum.ge.5) cycle
4700 pr(j)=c(j,inres)-cm(j)
4702 Im(1,1)=Im(1,1)+msc(iabs(iti),mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
4703 Im(1,2)=Im(1,2)-msc(iabs(iti),mnum)*pr(1)*pr(2)
4704 Im(1,3)=Im(1,3)-msc(iabs(iti),mnum)*pr(1)*pr(3)
4705 Im(2,3)=Im(2,3)-msc(iabs(iti),mnum)*pr(2)*pr(3)
4706 Im(2,2)=Im(2,2)+msc(iabs(iti),mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
4707 Im(3,3)=Im(3,3)+msc(iabs(iti),mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
4709 ! write(iout,*) "The angular momentum before Ip add"
4711 ! write (iout,*) (Im(i,j),j=1,3)
4716 Im(1,1)=Im(1,1)+Ip(mnum)*(1-dc_norm(1,i)*dc_norm(1,i))* &
4717 vbld(i+1)*vbld(i+1)*0.25d0
4718 Im(1,2)=Im(1,2)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(2,i))* &
4719 vbld(i+1)*vbld(i+1)*0.25d0
4720 Im(1,3)=Im(1,3)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(3,i))* &
4721 vbld(i+1)*vbld(i+1)*0.25d0
4722 Im(2,3)=Im(2,3)+Ip(mnum)*(-dc_norm(2,i)*dc_norm(3,i))* &
4723 vbld(i+1)*vbld(i+1)*0.25d0
4724 Im(2,2)=Im(2,2)+Ip(mnum)*(1-dc_norm(2,i)*dc_norm(2,i))* &
4725 vbld(i+1)*vbld(i+1)*0.25d0
4726 Im(3,3)=Im(3,3)+Ip(mnum)*(1-dc_norm(3,i)*dc_norm(3,i))* &
4727 vbld(i+1)*vbld(i+1)*0.25d0
4729 ! write(iout,*) "The angular momentum before Isc add"
4731 ! write (iout,*) (Im(i,j),j=1,3)
4737 ! if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
4738 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4739 .and.(mnum.lt.4)) then
4740 iti=iabs(itype(i,mnum))
4742 Im(1,1)=Im(1,1)+Isc(iti,mnum)*(1-dc_norm(1,inres)* &
4743 dc_norm(1,inres))*vbld(inres)*vbld(inres)
4744 Im(1,2)=Im(1,2)-Isc(iti,mnum)*(dc_norm(1,inres)* &
4745 dc_norm(2,inres))*vbld(inres)*vbld(inres)
4746 Im(1,3)=Im(1,3)-Isc(iti,mnum)*(dc_norm(1,inres)* &
4747 dc_norm(3,inres))*vbld(inres)*vbld(inres)
4748 Im(2,3)=Im(2,3)-Isc(iti,mnum)*(dc_norm(2,inres)* &
4749 dc_norm(3,inres))*vbld(inres)*vbld(inres)
4750 Im(2,2)=Im(2,2)+Isc(iti,mnum)*(1-dc_norm(2,inres)* &
4751 dc_norm(2,inres))*vbld(inres)*vbld(inres)
4752 Im(3,3)=Im(3,3)+Isc(iti,mnum)*(1-dc_norm(3,inres)* &
4753 dc_norm(3,inres))*vbld(inres)*vbld(inres)
4757 ! write(iout,*) "The angular momentum before agnom:"
4759 ! write (iout,*) (Im(i,j),j=1,3)
4763 ! write(iout,*) "The angular momentum before adjustment:"
4764 ! write(iout,*) (L(j),j=1,3)
4766 ! write (iout,*) (Im(i,j),j=1,3)
4772 ! Copying the Im matrix for the djacob subroutine
4780 ! Finding the eigenvectors and eignvalues of the inertia tensor
4781 call djacob(3,3,10000,1.0d-10,Imcp,eigvec,eigval)
4782 ! write (iout,*) "Eigenvalues & Eigenvectors"
4783 ! write (iout,'(5x,3f10.5)') (eigval(i),i=1,3)
4786 ! write (iout,'(i5,3f10.5)') i,(eigvec(i,j),j=1,3)
4788 ! Constructing the diagonalized matrix
4790 if (dabs(eigval(i)).gt.1.0d-15) then
4791 Id(i,i)=1.0d0/eigval(i)
4798 Imcp(i,j)=eigvec(j,i)
4804 pr1(i,j)=pr1(i,j)+Id(i,k)*Imcp(k,j)
4811 pr2(i,j)=pr2(i,j)+eigvec(i,k)*pr1(k,j)
4815 ! Calculating the total rotational velocity of the molecule
4818 vrot(i)=vrot(i)+pr2(i,j)*L(j)
4821 ! Resetting the velocities
4823 call vecpr(vrot(1),dc(1,i),vp)
4825 ! print *,"HERE2",d_t(j,i),vp(j)
4826 d_t(j,i)=d_t(j,i)-vp(j)
4827 ! print *,"HERE2",d_t(j,i),vp(j)
4832 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4833 .and.(mnum.lt.4)) then
4834 ! if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
4835 ! if(itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
4837 call vecpr(vrot(1),dc(1,inres),vp)
4839 d_t(j,inres)=d_t(j,inres)-vp(j)
4844 ! write(iout,*) "The angular momentum after adjustment:"
4845 ! write(iout,*) (L(j),j=1,3)
4848 end subroutine inertia_tensor
4849 !-----------------------------------------------------------------------------
4850 subroutine angmom(cm,L)
4853 ! implicit real*8 (a-h,o-z)
4854 ! include 'DIMENSIONS'
4855 ! include 'COMMON.CONTROL'
4856 ! include 'COMMON.VAR'
4857 ! include 'COMMON.MD'
4858 ! include 'COMMON.CHAIN'
4859 ! include 'COMMON.DERIV'
4860 ! include 'COMMON.GEO'
4861 ! include 'COMMON.LOCAL'
4862 ! include 'COMMON.INTERACT'
4863 ! include 'COMMON.IOUNITS'
4864 ! include 'COMMON.NAMES'
4865 real(kind=8) :: mscab
4866 real(kind=8),dimension(3) :: L,cm,pr,vp,vrot,incr,v,pp
4867 integer :: iti,inres,i,j,mnum
4868 ! Calculate the angular momentum
4877 if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
4879 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
4882 v(j)=incr(j)+0.5d0*d_t(j,i)
4885 incr(j)=incr(j)+d_t(j,i)
4887 call vecpr(pr(1),v(1),vp)
4889 L(j)=L(j)+mp(mnum)*vp(j)
4890 ! print *,"HERE3",J,i,L(j),mp(mnum),Ip(mnum),mnum
4894 pp(j)=0.5d0*d_t(j,i)
4896 call vecpr(pr(1),pp(1),vp)
4899 L(j)=L(j)+Ip(mnum)*vp(j)
4907 iti=iabs(itype(i,mnum))
4915 pr(j)=c(j,inres)-cm(j)
4918 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4919 .and.(mnum.lt.4)) then
4921 v(j)=incr(j)+d_t(j,inres)
4928 ! print *,i,pr,"pr",v
4929 call vecpr(pr(1),v(1),vp)
4930 ! write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),&
4931 ! " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
4933 L(j)=L(j)+mscab*vp(j)
4935 ! write (iout,*) "L",(l(j),j=1,3)
4936 ! print *,"L",(l(j),j=1,3),i,vp(1)
4938 if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4939 .and.(mnum.lt.4)) then
4941 v(j)=incr(j)+d_t(j,inres)
4943 call vecpr(dc(1,inres),d_t(1,inres),vp)
4945 L(j)=L(j)+Isc(iti,mnum)*vp(j)
4949 incr(j)=incr(j)+d_t(j,i)
4953 end subroutine angmom
4954 !-----------------------------------------------------------------------------
4955 subroutine vcm_vel(vcm)
4958 ! implicit real*8 (a-h,o-z)
4959 ! include 'DIMENSIONS'
4960 ! include 'COMMON.VAR'
4961 ! include 'COMMON.MD'
4962 ! include 'COMMON.CHAIN'
4963 ! include 'COMMON.DERIV'
4964 ! include 'COMMON.GEO'
4965 ! include 'COMMON.LOCAL'
4966 ! include 'COMMON.INTERACT'
4967 ! include 'COMMON.IOUNITS'
4968 real(kind=8),dimension(3) :: vcm,vv
4969 real(kind=8) :: summas,amas
4979 if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
4981 summas=summas+mp(mnum)
4983 vcm(j)=vcm(j)+mp(mnum)*(vv(j)+0.5d0*d_t(j,i))
4984 ! print *,i,j,vv(j),d_t(j,i)
4988 amas=msc(iabs(itype(i,mnum)),mnum)
4993 if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4994 .and.(mnum.lt.4)) then
4996 vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres))
5000 vcm(j)=vcm(j)+amas*vv(j)
5004 vv(j)=vv(j)+d_t(j,i)
5007 ! write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas
5009 vcm(j)=vcm(j)/summas
5012 end subroutine vcm_vel
5014 !-----------------------------------------------------------------------------
5016 !-----------------------------------------------------------------------------
5018 ! RATTLE algorithm for velocity Verlet - step 1, UNRES
5022 ! implicit real*8 (a-h,o-z)
5023 ! include 'DIMENSIONS'
5025 ! include 'COMMON.CONTROL'
5026 ! include 'COMMON.VAR'
5027 ! include 'COMMON.MD'
5029 ! include 'COMMON.LANGEVIN'
5031 ! include 'COMMON.LANGEVIN.lang0'
5033 ! include 'COMMON.CHAIN'
5034 ! include 'COMMON.DERIV'
5035 ! include 'COMMON.GEO'
5036 ! include 'COMMON.LOCAL'
5037 ! include 'COMMON.INTERACT'
5038 ! include 'COMMON.IOUNITS'
5039 ! include 'COMMON.NAMES'
5040 ! include 'COMMON.TIME1'
5041 !el real(kind=8) :: gginv(2*nres,2*nres),&
5042 !el gdc(3,2*nres,2*nres)
5043 real(kind=8) :: dC_uncor(3,2*nres) !,&
5044 !el real(kind=8) :: Cmat(2*nres,2*nres)
5045 real(kind=8) :: x(2*nres),xcorr(3,2*nres) !maxres2=2*maxres
5046 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
5047 !el common /przechowalnia/ nbond
5048 integer :: max_rattle = 5
5049 logical :: lprn = .false., lprn1 = .false., not_done
5050 real(kind=8) :: tol_rattle = 1.0d-5
5052 integer :: ii,i,j,jj,l,ind,ind1,nres2
5055 !el /common/ przechowalnia
5057 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
5058 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
5059 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
5061 if (lprn) write (iout,*) "RATTLE1"
5065 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5066 .and.(mnum.lt.4)) nbond=nbond+1
5068 ! Make a folded form of the Ginv-matrix
5081 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
5086 if (itype(k,1).ne.10 .and. itype(k,mnum).ne.ntyp1_molec(mnum)&
5087 .and.(mnum.lt.4)) then
5091 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
5099 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5110 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
5114 if (itype(k,1).ne.10) then
5118 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
5126 write (iout,*) "Matrix GGinv"
5127 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
5133 if (iter.gt.max_rattle) then
5134 write (iout,*) "Error - too many iterations in RATTLE."
5137 ! Calculate the matrix C = GG**(-1) dC_old o dC
5142 dC_uncor(j,ind1)=dC(j,i)
5146 if (itype(i,1).ne.10) then
5149 dC_uncor(j,ind1)=dC(j,i+nres)
5158 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
5162 if (itype(k,1).ne.10) then
5165 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
5170 ! Calculate deviations from standard virtual-bond lengths
5174 x(ind)=vbld(i+1)**2-vbl**2
5177 if (itype(i,1).ne.10) then
5179 x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
5183 write (iout,*) "Coordinates and violations"
5185 write(iout,'(i5,3f10.5,5x,e15.5)') &
5186 i,(dC_uncor(j,i),j=1,3),x(i)
5188 write (iout,*) "Velocities and violations"
5192 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5193 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
5197 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5198 .and.(mnum.lt.4)) then
5201 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5202 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
5203 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
5206 ! write (iout,*) "gdc"
5208 ! write (iout,*) "i",i
5210 ! write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
5216 if (dabs(x(i)).gt.xmax) then
5220 if (xmax.lt.tol_rattle) then
5224 ! Calculate the matrix of the system of equations
5229 Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
5234 write (iout,*) "Matrix Cmat"
5235 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
5237 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
5238 ! Add constraint term to positions
5245 xx = xx+x(ii)*gdc(j,ind,ii)
5249 d_t_new(j,i)=d_t_new(j,i)-xx/d_time
5253 if (itype(i,1).ne.10) then
5258 xx = xx+x(ii)*gdc(j,ind,ii)
5261 dC(j,i+nres)=dC(j,i+nres)-xx
5262 d_t_new(j,i+nres)=d_t_new(j,i+nres)-xx/d_time
5266 ! Rebuild the chain using the new coordinates
5267 call chainbuild_cart
5269 write (iout,*) "New coordinates, Lagrange multipliers,",&
5270 " and differences between actual and standard bond lengths"
5274 xx=vbld(i+1)**2-vbl**2
5275 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
5276 i,(dC(j,i),j=1,3),x(ind),xx
5280 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5283 xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
5284 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
5285 i,(dC(j,i+nres),j=1,3),x(ind),xx
5288 write (iout,*) "Velocities and violations"
5292 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5293 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
5296 if (itype(i,1).ne.10) then
5298 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5299 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
5300 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
5307 10 write (iout,*) "Error - singularity in solving the system",&
5308 " of equations for Lagrange multipliers."
5312 "RATTLE inactive; use -DRATTLE switch at compile time."
5315 end subroutine rattle1
5316 !-----------------------------------------------------------------------------
5318 ! RATTLE algorithm for velocity Verlet - step 2, UNRES
5322 ! implicit real*8 (a-h,o-z)
5323 ! include 'DIMENSIONS'
5325 ! include 'COMMON.CONTROL'
5326 ! include 'COMMON.VAR'
5327 ! include 'COMMON.MD'
5329 ! include 'COMMON.LANGEVIN'
5331 ! include 'COMMON.LANGEVIN.lang0'
5333 ! include 'COMMON.CHAIN'
5334 ! include 'COMMON.DERIV'
5335 ! include 'COMMON.GEO'
5336 ! include 'COMMON.LOCAL'
5337 ! include 'COMMON.INTERACT'
5338 ! include 'COMMON.IOUNITS'
5339 ! include 'COMMON.NAMES'
5340 ! include 'COMMON.TIME1'
5341 !el real(kind=8) :: gginv(2*nres,2*nres),&
5342 !el gdc(3,2*nres,2*nres)
5343 real(kind=8) :: dC_uncor(3,2*nres) !,&
5344 !el Cmat(2*nres,2*nres)
5345 real(kind=8) :: x(2*nres) !maxres2=2*maxres
5346 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
5347 !el common /przechowalnia/ nbond
5348 integer :: max_rattle = 5
5349 logical :: lprn = .false., lprn1 = .false., not_done
5350 real(kind=8) :: tol_rattle = 1.0d-5
5354 !el /common/ przechowalnia
5355 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
5356 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
5357 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
5359 if (lprn) write (iout,*) "RATTLE2"
5360 if (lprn) write (iout,*) "Velocity correction"
5361 ! Calculate the matrix G dC
5367 gdc(j,i,ind)=GGinv(i,ind)*dC(j,k)
5372 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5373 .and.(mnum.lt.4)) then
5376 gdc(j,i,ind)=GGinv(i,ind)*dC(j,k+nres)
5382 ! write (iout,*) "gdc"
5384 ! write (iout,*) "i",i
5386 ! write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
5390 ! Calculate the matrix of the system of equations
5397 Cmat(ind,j)=Cmat(ind,j)+dC(k,i)*gdc(k,ind,j)
5403 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5404 .and.(mnum.lt.4)) then
5409 Cmat(ind,j)=Cmat(ind,j)+dC(k,i+nres)*gdc(k,ind,j)
5414 ! Calculate the scalar product dC o d_t_new
5418 x(ind)=scalar(d_t(1,i),dC(1,i))
5422 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5423 .and.(mnum.lt.4)) then
5425 x(ind)=scalar(d_t(1,i+nres),dC(1,i+nres))
5429 write (iout,*) "Velocities and violations"
5433 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5434 i,ind,(d_t(j,i),j=1,3),x(ind)
5438 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5439 .and.(mnum.lt.4)) then
5441 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5442 i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind)
5448 if (dabs(x(i)).gt.xmax) then
5452 if (xmax.lt.tol_rattle) then
5457 write (iout,*) "Matrix Cmat"
5458 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
5460 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
5461 ! Add constraint term to velocities
5468 xx = xx+x(ii)*gdc(j,ind,ii)
5470 d_t(j,i)=d_t(j,i)-xx
5475 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5476 .and.(mnum.lt.4)) then
5481 xx = xx+x(ii)*gdc(j,ind,ii)
5483 d_t(j,i+nres)=d_t(j,i+nres)-xx
5489 "New velocities, Lagrange multipliers violations"
5493 if (lprn) write (iout,'(2i5,3f10.5,5x,2e15.5)') &
5494 i,ind,(d_t(j,i),j=1,3),x(ind),scalar(d_t(1,i),dC(1,i))
5498 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5501 write (iout,'(2i5,3f10.5,5x,2e15.5)') &
5502 i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind),&
5503 scalar(d_t(1,i+nres),dC(1,i+nres))
5509 10 write (iout,*) "Error - singularity in solving the system",&
5510 " of equations for Lagrange multipliers."
5514 "RATTLE inactive; use -DRATTLE option at compile time."
5517 end subroutine rattle2
5518 !-----------------------------------------------------------------------------
5519 subroutine rattle_brown
5520 ! RATTLE/LINCS algorithm for Brownian dynamics, UNRES
5524 ! implicit real*8 (a-h,o-z)
5525 ! include 'DIMENSIONS'
5527 ! include 'COMMON.CONTROL'
5528 ! include 'COMMON.VAR'
5529 ! include 'COMMON.MD'
5531 ! include 'COMMON.LANGEVIN'
5533 ! include 'COMMON.LANGEVIN.lang0'
5535 ! include 'COMMON.CHAIN'
5536 ! include 'COMMON.DERIV'
5537 ! include 'COMMON.GEO'
5538 ! include 'COMMON.LOCAL'
5539 ! include 'COMMON.INTERACT'
5540 ! include 'COMMON.IOUNITS'
5541 ! include 'COMMON.NAMES'
5542 ! include 'COMMON.TIME1'
5543 !el real(kind=8) :: gginv(2*nres,2*nres),&
5544 !el gdc(3,2*nres,2*nres)
5545 real(kind=8) :: dC_uncor(3,2*nres) !,&
5546 !el real(kind=8) :: Cmat(2*nres,2*nres)
5547 real(kind=8) :: x(2*nres) !maxres2=2*maxres
5548 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
5549 !el common /przechowalnia/ nbond
5550 integer :: max_rattle = 5
5551 logical :: lprn = .false., lprn1 = .false., not_done
5552 real(kind=8) :: tol_rattle = 1.0d-5
5556 !el /common/ przechowalnia
5557 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
5558 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
5559 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
5562 if (lprn) write (iout,*) "RATTLE_BROWN"
5565 if (itype(i,1).ne.10) nbond=nbond+1
5567 ! Make a folded form of the Ginv-matrix
5580 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
5584 if (itype(k,1).ne.10) then
5588 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
5595 if (itype(i,1).ne.10) then
5605 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
5609 if (itype(k,1).ne.10) then
5613 if (j.eq.1 .and. l.eq.1)GGinv(ii,jj)=fricmat(ind,ind1)
5621 write (iout,*) "Matrix GGinv"
5622 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
5628 if (iter.gt.max_rattle) then
5629 write (iout,*) "Error - too many iterations in RATTLE."
5632 ! Calculate the matrix C = GG**(-1) dC_old o dC
5637 dC_uncor(j,ind1)=dC(j,i)
5641 if (itype(i,1).ne.10) then
5644 dC_uncor(j,ind1)=dC(j,i+nres)
5653 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
5657 if (itype(k,1).ne.10) then
5660 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
5665 ! Calculate deviations from standard virtual-bond lengths
5669 x(ind)=vbld(i+1)**2-vbl**2
5672 if (itype(i,1).ne.10) then
5674 x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
5678 write (iout,*) "Coordinates and violations"
5680 write(iout,'(i5,3f10.5,5x,e15.5)') &
5681 i,(dC_uncor(j,i),j=1,3),x(i)
5683 write (iout,*) "Velocities and violations"
5687 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5688 i,ind,(d_t(j,i),j=1,3),scalar(d_t(1,i),dC_old(1,i))
5691 if (itype(i,1).ne.10) then
5693 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5694 i+nres,ind,(d_t(j,i+nres),j=1,3),&
5695 scalar(d_t(1,i+nres),dC_old(1,i+nres))
5698 write (iout,*) "gdc"
5700 write (iout,*) "i",i
5702 write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
5708 if (dabs(x(i)).gt.xmax) then
5712 if (xmax.lt.tol_rattle) then
5716 ! Calculate the matrix of the system of equations
5721 Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
5726 write (iout,*) "Matrix Cmat"
5727 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
5729 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
5730 ! Add constraint term to positions
5737 xx = xx+x(ii)*gdc(j,ind,ii)
5740 d_t(j,i)=d_t(j,i)+xx/d_time
5745 if (itype(i,1).ne.10) then
5750 xx = xx+x(ii)*gdc(j,ind,ii)
5753 d_t(j,i+nres)=d_t(j,i+nres)+xx/d_time
5754 dC(j,i+nres)=dC(j,i+nres)+xx
5758 ! Rebuild the chain using the new coordinates
5759 call chainbuild_cart
5761 write (iout,*) "New coordinates, Lagrange multipliers,",&
5762 " and differences between actual and standard bond lengths"
5766 xx=vbld(i+1)**2-vbl**2
5767 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
5768 i,(dC(j,i),j=1,3),x(ind),xx
5771 if (itype(i,1).ne.10) then
5773 xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
5774 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
5775 i,(dC(j,i+nres),j=1,3),x(ind),xx
5778 write (iout,*) "Velocities and violations"
5782 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5783 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
5786 if (itype(i,1).ne.10) then
5788 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5789 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
5790 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
5797 10 write (iout,*) "Error - singularity in solving the system",&
5798 " of equations for Lagrange multipliers."
5802 "RATTLE inactive; use -DRATTLE option at compile time"
5805 end subroutine rattle_brown
5806 !-----------------------------------------------------------------------------
5808 !-----------------------------------------------------------------------------
5809 subroutine friction_force
5814 ! implicit real*8 (a-h,o-z)
5815 ! include 'DIMENSIONS'
5816 ! include 'COMMON.VAR'
5817 ! include 'COMMON.CHAIN'
5818 ! include 'COMMON.DERIV'
5819 ! include 'COMMON.GEO'
5820 ! include 'COMMON.LOCAL'
5821 ! include 'COMMON.INTERACT'
5822 ! include 'COMMON.MD'
5824 ! include 'COMMON.LANGEVIN'
5826 ! include 'COMMON.LANGEVIN.lang0'
5828 ! include 'COMMON.IOUNITS'
5829 !el real(kind=8),dimension(6*nres) :: gamvec !(MAXRES6) maxres6=6*maxres
5830 !el common /syfek/ gamvec
5832 integer iposc,ichain,n,innt,inct
5833 double precision rs(nres*2)
5834 real(kind=8) ::v_work(3,6*nres),vvec(2*nres)
5836 real(kind=8) :: v_work(6*nres)
5839 real(kind=8) :: vv(3),vvtot(3,nres)!,v_work(6*nres) !,&
5840 !el ginvfric(2*nres,2*nres) !maxres2=2*maxres
5841 !el common /przechowalnia/ ginvfric
5843 logical :: lprn, checkmode
5844 integer :: i,j,ind,k,nres2,nres6,mnum
5849 ! if (large) lprn=.true.
5850 ! if (large) checkmode=.true.
5852 !c Here accelerations due to friction forces are computed right after forces.
5853 d_t_work(:6*nres)=0.0d0
5855 v_work(j,1)=d_t(j,0)
5856 v_work(j,nnt)=d_t(j,0)
5860 v_work(j,i)=v_work(j,i-1)+d_t(j,i-1)
5865 if (iabs(itype(i,1)).ne.10 .and. iabs(itype(i,mnum)).ne.ntyp1_molec(mnum).and.mnum.lt.3) then
5867 v_work(j,i+nres)=v_work(j,i)+d_t(j,i+nres)
5872 write (iout,*) "v_work"
5874 write (iout,'(i5,3f10.5)') i,(v_work(j,i),j=1,3)
5880 n=dimen_chain(ichain)
5881 iposc=iposd_chain(ichain)
5882 !c write (iout,*) "friction_force j",j," ichain",ichain,
5883 !c & " n",n," iposc",iposc,iposc+n-1
5884 innt=chain_border(1,ichain)
5885 inct=chain_border(2,ichain)
5887 !c innt=chain_border(1,1)
5888 !c inct=chain_border(2,1)
5890 vvec(ind+1)=v_work(j,i)
5892 ! if (iabs(itype(i)).ne.10) then
5893 if (iabs(itype(i,1)).ne.10 .and. iabs(itype(i,mnum)).ne.ntyp1_molec(mnum).and.mnum.lt.3) then
5894 vvec(ind+1)=v_work(j,i+nres)
5899 write (iout,*) "vvec ind",ind," n",n
5900 write (iout,'(f10.5)') (vvec(i),i=iposc,ind)
5902 !c write (iout,*) "chain",i," ind",ind," n",n
5910 ! if (large) print *,"before fivediagmult"
5911 call fivediagmult(n,DMfric(iposc),DU1fric(iposc),&
5912 DU2fric(iposc),vvec(iposc),rs)
5913 ! if (large) print *,"after fivediagmult"
5917 time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
5919 time_fricmatmult=time_fricmatmult+tcpu()-time01
5924 write (iout,'(f10.5)') (rs(i),i=1,n)
5926 do i=iposc,iposc+n-1
5927 ! write (iout,*) "ichain",ichain," i",i," j",j,&
5928 ! "index",3*(i-1)+j,"rs",rs(i-iposc+1)
5929 fric_work(3*(i-1)+j)=-rs(i-iposc+1)
5934 write (iout,*) "Vector fric_work dimen3",dimen3
5935 write (iout,'(3f10.5)') (fric_work(j),j=1,dimen3)
5938 if(.not.allocated(gamvec)) allocate(gamvec(nres6)) !(MAXRES6)
5939 if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
5947 d_t_work(j)=d_t(j,0)
5952 d_t_work(ind+j)=d_t(j,i)
5958 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5959 .and.(mnum.lt.4)) then
5961 d_t_work(ind+j)=d_t(j,i+nres)
5967 call fricmat_mult(d_t_work,fric_work)
5969 if (.not.checkmode) return
5972 write (iout,*) "d_t_work and fric_work"
5974 write (iout,'(i3,2e15.5)') i,d_t_work(i),fric_work(i)
5978 friction(j,0)=fric_work(j)
5983 friction(j,i)=fric_work(ind+j)
5989 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5990 .and.(mnum.lt.4)) then
5991 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
5993 friction(j,i+nres)=fric_work(ind+j)
5999 write(iout,*) "Friction backbone"
6001 write(iout,'(i5,3e15.5,5x,3e15.5)') &
6002 i,(friction(j,i),j=1,3),(d_t(j,i),j=1,3)
6004 write(iout,*) "Friction side chain"
6006 write(iout,'(i5,3e15.5,5x,3e15.5)') &
6007 i,(friction(j,i+nres),j=1,3),(d_t(j,i+nres),j=1,3)
6016 vvtot(j,i)=vv(j)+0.5d0*d_t(j,i)
6017 vvtot(j,i+nres)=vv(j)+d_t(j,i+nres)
6018 vv(j)=vv(j)+d_t(j,i)
6021 write (iout,*) "vvtot backbone and sidechain"
6023 write (iout,'(i5,3e15.5,5x,3e15.5)') i,(vvtot(j,i),j=1,3),&
6024 (vvtot(j,i+nres),j=1,3)
6029 v_work(ind+j)=vvtot(j,i)
6035 v_work(ind+j)=vvtot(j,i+nres)
6039 write (iout,*) "v_work gamvec and site-based friction forces"
6041 write (iout,'(i5,3e15.5)') i,v_work(i),gamvec(i),&
6045 ! fric_work1(i)=0.0d0
6047 ! fric_work1(i)=fric_work1(i)-A(j,i)*gamvec(j)*v_work(j)
6050 ! write (iout,*) "fric_work and fric_work1"
6052 ! write (iout,'(i5,2e15.5)') i,fric_work(i),fric_work1(i)
6058 ginvfric(i,j)=ginvfric(i,j)+ginv(i,k)*fricmat(k,j)
6062 write (iout,*) "ginvfric"
6064 write (iout,'(i5,100f8.3)') i,(ginvfric(i,j),j=1,dimen)
6066 write (iout,*) "symmetry check"
6069 write (iout,*) i,j,ginvfric(i,j)-ginvfric(j,i)
6075 end subroutine friction_force
6076 !-----------------------------------------------------------------------------
6077 subroutine setup_fricmat
6081 use control_data, only:time_Bcast
6082 use control, only:tcpu
6084 ! implicit real*8 (a-h,o-z)
6088 real(kind=8) :: time00
6090 ! include 'DIMENSIONS'
6091 ! include 'COMMON.VAR'
6092 ! include 'COMMON.CHAIN'
6093 ! include 'COMMON.DERIV'
6094 ! include 'COMMON.GEO'
6095 ! include 'COMMON.LOCAL'
6096 ! include 'COMMON.INTERACT'
6097 ! include 'COMMON.MD'
6098 ! include 'COMMON.SETUP'
6099 ! include 'COMMON.TIME1'
6100 ! integer licznik /0/
6103 ! include 'COMMON.LANGEVIN'
6105 ! include 'COMMON.LANGEVIN.lang0'
6107 ! include 'COMMON.IOUNITS'
6109 integer :: i,j,ind,ind1,m,ichain,innt,inct
6110 logical :: lprn = .false.
6111 real(kind=8) :: dtdi !el ,gamvec(2*nres)
6112 !el real(kind=8),dimension(2*nres,2*nres) :: ginvfric,fcopy
6113 ! real(kind=8),allocatable,dimension(:,:) :: fcopy
6114 !el real(kind=8),dimension(2*nres*(2*nres+1)/2) :: Ghalf !(mmaxres2) (mmaxres2=(maxres2*(maxres2+1)/2))
6115 !el common /syfek/ gamvec
6116 real(kind=8) :: work(8*2*nres)
6117 integer :: iwork(2*nres)
6118 !el common /przechowalnia/ ginvfric,Ghalf,fcopy
6119 integer :: ii,iti,k,l,nzero,nres2,nres6,ierr,mnum
6124 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
6125 if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
6126 if (fg_rank.ne.king) goto 10
6132 if(.not.allocated(gamvec)) allocate(gamvec(nres2)) !(MAXRES2)
6134 if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
6135 if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
6136 !el allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
6137 if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !maxres2=2*maxres
6139 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
6142 if (.not.allocated(DMfric)) allocate(DMfric(nres2))
6143 if (.not.allocated(DU1fric)) allocate(DU1fric(nres2))
6144 if (.not.allocated(DU2fric)) allocate(DU2fric(nres2))
6145 ! Load the friction coefficients corresponding to peptide groups
6150 gamvec(ind1)=gamp(mnum)
6153 if (molnum(nct).eq.5) then
6156 gamvec(ind1)=gamp(mnum)
6158 ! Load the friction coefficients corresponding to side chains
6162 gamsc(ntyp1_molec(j),j)=1.0d0
6169 gamvec(ii)=gamsc(iabs(iti),mnum)
6171 if (surfarea) call sdarea(gamvec)
6177 innt=chain_border(1,ichain)
6178 inct=chain_border(2,ichain)
6179 !c write (iout,*) "ichain",ichain," innt",innt," inct",inct
6182 DMfric(ind)=gamvec(innt-nnt+1)/4
6183 if (iabs(itype(innt,1)).eq.10.or.mnum.gt.2) then
6184 DMfric(ind)=DMfric(ind)+gamvec(m+innt-nnt+1)
6187 DMfric(ind+1)=gamvec(m+innt-nnt+1)
6190 !c write (iout,*) "DMfric init ind",ind
6194 DMfric(ind)=gamvec(i-nnt+1)/2
6195 if (iabs(itype(i,1)).eq.10.or.mnum.gt.2) then
6196 DMfric(ind)=DMfric(ind)+gamvec(m+i-nnt+1)
6199 DMfric(ind+1)=gamvec(m+i-nnt+1)
6203 !c write (iout,*) "DMfric endloop ind",ind
6204 if (inct.gt.innt) then
6205 DMfric(ind)=gamvec(inct-1-nnt+1)/4
6207 if (iabs(itype(inct,1)).eq.10.or.mnum.gt.2) then
6208 DMfric(ind)=DMfric(ind)+gamvec(inct+m-nnt+1)
6211 DMfric(ind+1)=gamvec(inct+m-nnt+1)
6215 !c write (iout,*) "DMfric end ind",ind
6221 ind=iposd_chain(ichain)
6222 innt=chain_border(1,ichain)
6223 inct=chain_border(2,ichain)
6225 if (iabs(itype(i,1)).ne.10.and.mnum.le.2) then
6228 DU1fric(ind)=gamvec(i-nnt+1)/4
6236 ind=iposd_chain(ichain)
6237 innt=chain_border(1,ichain)
6238 inct=chain_border(2,ichain)
6240 if (iabs(itype(i,1)).ne.10.and.mnum.le.2) then
6241 DU2fric(ind)=gamvec(i-nnt+1)/4
6242 DU2fric(ind+1)=0.0d0
6251 write(iout,*)"The upper part of the five-diagonal friction matrix"
6253 write (iout,'(a,i5)') 'Chain',ichain
6254 innt=iposd_chain(ichain)
6255 inct=iposd_chain(ichain)+dimen_chain(ichain)-1
6257 if (i.lt.inct-1) then
6258 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMfric(i),DU1fric(i),&
6260 else if (i.eq.inct-1) then
6261 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMfric(i),DU1fric(i)
6263 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMfric(i)
6272 ! Zeroing out fricmat
6278 ! Load the friction coefficients corresponding to peptide groups
6283 gamvec(ind1)=gamp(mnum)
6286 if (molnum(nct).eq.5) then
6289 gamvec(ind1)=gamp(mnum)
6291 ! Load the friction coefficients corresponding to side chains
6295 gamsc(ntyp1_molec(j),j)=1.0d0
6302 gamvec(ii)=gamsc(iabs(iti),mnum)
6304 if (surfarea) call sdarea(gamvec)
6306 ! write (iout,*) "Matrix A and vector gamma"
6308 ! write (iout,'(i2,$)') i
6310 ! write (iout,'(f4.1,$)') A(i,j)
6312 ! write (iout,'(f8.3)') gamvec(i)
6316 write (iout,*) "Vector gamvec"
6318 write (iout,'(i5,f10.5)') i, gamvec(i)
6322 ! The friction matrix
6327 dtdi=dtdi+A(j,k)*A(j,i)*gamvec(j)
6334 write (iout,'(//a)') "Matrix fricmat"
6335 call matout2(dimen,dimen,nres2,nres2,fricmat)
6337 if (lang.eq.2 .or. lang.eq.3) then
6338 ! Mass-scale the friction matrix if non-direct integration will be performed
6344 Ginvfric(i,j)=Ginvfric(i,j)+ &
6345 Gsqrm(i,k)*Gsqrm(l,j)*fricmat(k,l)
6350 ! Diagonalize the friction matrix
6355 Ghalf(ind)=Ginvfric(i,j)
6358 call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
6361 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
6362 " mass-scaled friction matrix"
6363 call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
6365 ! Precompute matrices for tinker stochastic integrator
6372 mt1(i,j)=mt1(i,j)+fricvec(k,i)*gsqrm(k,j)
6373 mt2(i,j)=mt2(i,j)+fricvec(k,i)*gsqrp(k,j)
6379 else if (lang.eq.4) then
6380 ! Diagonalize the friction matrix
6385 Ghalf(ind)=fricmat(i,j)
6388 call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
6391 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
6393 call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
6395 ! Determine the number of zero eigenvalues of the friction matrix
6396 nzero=max0(dimen-dimen1,0)
6397 ! do while (fricgam(nzero+1).le.1.0d-5 .and. nzero.lt.dimen)
6400 write (iout,*) "Number of zero eigenvalues:",nzero
6405 fricmat(i,j)=fricmat(i,j) &
6406 +fricvec(i,k)*fricvec(j,k)/fricgam(k)
6411 write (iout,'(//a)') "Generalized inverse of fricmat"
6412 call matout(dimen,dimen,nres6,nres6,fricmat)
6417 if (nfgtasks.gt.1) then
6418 if (fg_rank.eq.0) then
6419 ! The matching BROADCAST for fg processors is called in ERGASTULUM
6425 call MPI_Bcast(10,1,MPI_INTEGER,king,FG_COMM,IERROR)
6427 time_Bcast=time_Bcast+MPI_Wtime()-time00
6429 time_Bcast=time_Bcast+tcpu()-time00
6431 ! print *,"Processor",myrank,
6432 ! & " BROADCAST iorder in SETUP_FRICMAT"
6435 write (iout,*) "setup_fricmat licznik"!,licznik !sp
6441 ! Scatter the friction matrix
6442 call MPI_Scatterv(fricmat(1,1),nginv_counts(0),&
6443 nginv_start(0),MPI_DOUBLE_PRECISION,fcopy(1,1),&
6444 myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
6447 time_scatter=time_scatter+MPI_Wtime()-time00
6448 time_scatter_fmat=time_scatter_fmat+MPI_Wtime()-time00
6450 time_scatter=time_scatter+tcpu()-time00
6451 time_scatter_fmat=time_scatter_fmat+tcpu()-time00
6455 do j=1,2*my_ng_count
6456 fricmat(j,i)=fcopy(i,j)
6459 ! write (iout,*) "My chunk of fricmat"
6460 ! call MATOUT2(my_ng_count,dimen,maxres2,maxres2,fcopy)
6465 end subroutine setup_fricmat
6466 !-----------------------------------------------------------------------------
6467 subroutine stochastic_force(stochforcvec)
6470 use random, only:anorm_distr
6471 ! implicit real*8 (a-h,o-z)
6472 ! include 'DIMENSIONS'
6473 use control, only: tcpu
6478 ! include 'COMMON.VAR'
6479 ! include 'COMMON.CHAIN'
6480 ! include 'COMMON.DERIV'
6481 ! include 'COMMON.GEO'
6482 ! include 'COMMON.LOCAL'
6483 ! include 'COMMON.INTERACT'
6484 ! include 'COMMON.MD'
6485 ! include 'COMMON.TIME1'
6487 ! include 'COMMON.LANGEVIN'
6489 ! include 'COMMON.LANGEVIN.lang0'
6491 ! include 'COMMON.IOUNITS'
6493 real(kind=8) :: x,sig,lowb,highb
6494 real(kind=8) :: ff(3),force(3,0:2*nres),zeta2,lowb2
6495 real(kind=8) :: highb2,sig2,forcvec(6*nres),stochforcvec(6*nres)
6496 real(kind=8) :: time00
6497 logical :: lprn = .false.
6498 integer :: i,j,ind,mnum
6500 integer ichain,innt,inct,iposc
6505 stochforc(j,i)=0.0d0
6515 ! Compute the stochastic forces acting on bodies. Store in force.
6521 force(j,i)=anorm_distr(x,sig,lowb,highb)
6529 force(j,i+nres)=anorm_distr(x,sig2,lowb2,highb2)
6533 time_fsample=time_fsample+MPI_Wtime()-time00
6535 time_fsample=time_fsample+tcpu()-time00
6540 innt=chain_border(1,ichain)
6541 inct=chain_border(2,ichain)
6542 iposc=iposd_chain(ichain)
6543 !c for debugging only
6544 !c innt=chain_border(1,1)
6545 !c inct=chain_border(2,1)
6546 !c iposc=iposd_chain(1)
6547 !c write (iout,*)"stochastic_force ichain=",ichain," innt",innt,
6548 !c & " inct",inct," iposc",iposc
6550 stochforcvec(ind+j)=0.5d0*force(j,innt)
6552 if (iabs(itype(innt,molnum(innt))).eq.10.or.molnum(innt).ge.3) then
6554 stochforcvec(ind+j)=stochforcvec(ind+j)+force(j,innt+nres)
6560 stochforcvec(ind+j)=force(j,innt+nres)
6566 stochforcvec(ind+j)=0.5d0*(force(j,i)+force(j,i-1))
6568 if (iabs(itype(i,1).eq.10).or.molnum(i).ge.3) then
6570 stochforcvec(ind+j)=stochforcvec(ind+j)+force(j,i+nres)
6576 stochforcvec(ind+j)=force(j,i+nres)
6582 stochforcvec(ind+j)=0.5d0*force(j,inct-1)
6584 if (iabs(itype(inct,molnum(inct))).eq.10.or.molnum(inct).ge.3) then
6586 stochforcvec(ind+j)=stochforcvec(ind+j)+force(j,inct+nres)
6592 stochforcvec(ind+j)=force(j,inct+nres)
6596 !c write (iout,*) "chain",ichain," ind",ind
6599 write (iout,*) "stochforcvec"
6600 write (iout,'(3f10.5)') (stochforcvec(j),j=1,ind)
6603 ! Compute the stochastic forces acting on virtual-bond vectors.
6609 stochforc(j,i)=ff(j)+0.5d0*force(j,i)
6612 ff(j)=ff(j)+force(j,i)
6614 ! if (itype(i+1,1).ne.ntyp1) then
6616 if (itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
6618 stochforc(j,i)=stochforc(j,i)+force(j,i+nres+1)
6619 ff(j)=ff(j)+force(j,i+nres+1)
6624 stochforc(j,0)=ff(j)+force(j,nnt+nres)
6628 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
6629 .and.(mnum.lt.4)) then
6630 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
6632 stochforc(j,i+nres)=force(j,i+nres)
6638 stochforcvec(j)=stochforc(j,0)
6643 stochforcvec(ind+j)=stochforc(j,i)
6649 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
6650 .and.(mnum.lt.4)) then
6651 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
6653 stochforcvec(ind+j)=stochforc(j,i+nres)
6659 write (iout,*) "stochforcvec"
6661 write(iout,'(i5,e15.5)') i,stochforcvec(i)
6663 write(iout,*) "Stochastic forces backbone"
6665 write(iout,'(i5,3e15.5)') i,(stochforc(j,i),j=1,3)
6667 write(iout,*) "Stochastic forces side chain"
6669 write(iout,'(i5,3e15.5)') &
6670 i,(stochforc(j,i+nres),j=1,3)
6678 write (iout,*) i,ind
6680 forcvec(ind+j)=force(j,i)
6685 write (iout,*) i,ind
6687 forcvec(j+ind)=force(j,i+nres)
6692 write (iout,*) "forcvec"
6696 write (iout,'(2i3,2f10.5)') i,j,force(j,i),&
6703 write (iout,'(2i3,2f10.5)') i,j,force(j,i+nres),&
6712 end subroutine stochastic_force
6713 !-----------------------------------------------------------------------------
6714 subroutine sdarea(gamvec)
6716 ! Scale the friction coefficients according to solvent accessible surface areas
6717 ! Code adapted from TINKER
6721 ! implicit real*8 (a-h,o-z)
6722 ! include 'DIMENSIONS'
6723 ! include 'COMMON.CONTROL'
6724 ! include 'COMMON.VAR'
6725 ! include 'COMMON.MD'
6727 ! include 'COMMON.LANGEVIN'
6729 ! include 'COMMON.LANGEVIN.lang0'
6731 ! include 'COMMON.CHAIN'
6732 ! include 'COMMON.DERIV'
6733 ! include 'COMMON.GEO'
6734 ! include 'COMMON.LOCAL'
6735 ! include 'COMMON.INTERACT'
6736 ! include 'COMMON.IOUNITS'
6737 ! include 'COMMON.NAMES'
6738 real(kind=8),dimension(2*nres) :: radius,gamvec !(maxres2)
6739 real(kind=8),parameter :: twosix = 1.122462048309372981d0
6740 logical :: lprn = .false.
6741 real(kind=8) :: probe,area,ratio
6742 integer :: i,j,ind,iti,mnum
6744 ! determine new friction coefficients every few SD steps
6746 ! set the atomic radii to estimates of sigma values
6748 ! print *,"Entered sdarea"
6754 ! Load peptide group radii
6757 radius(i)=pstok(mnum)
6759 ! Load side chain radii
6763 radius(i+nres)=restok(iti,mnum)
6766 ! write (iout,*) "i",i," radius",radius(i)
6769 radius(i) = radius(i) / twosix
6770 if (radius(i) .ne. 0.0d0) radius(i) = radius(i) + probe
6773 ! scale atomic friction coefficients by accessible area
6775 if (lprn) write (iout,*) &
6776 "Original gammas, surface areas, scaling factors, new gammas, ",&
6777 "std's of stochastic forces"
6780 if (radius(i).gt.0.0d0) then
6781 call surfatom (i,area,radius)
6782 ratio = dmax1(area/(4.0d0*pi*radius(i)**2),1.0d-1)
6783 if (lprn) write (iout,'(i5,3f10.5,$)') &
6784 i,gamvec(ind+1),area,ratio
6787 gamvec(ind) = ratio * gamvec(ind)
6790 stdforcp(i)=stdfp(mnum)*dsqrt(gamvec(ind))
6791 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcp(i)
6795 if (radius(i+nres).gt.0.0d0) then
6796 call surfatom (i+nres,area,radius)
6797 ratio = dmax1(area/(4.0d0*pi*radius(i+nres)**2),1.0d-1)
6798 if (lprn) write (iout,'(i5,3f10.5,$)') &
6799 i,gamvec(ind+1),area,ratio
6802 gamvec(ind) = ratio * gamvec(ind)
6805 stdforcsc(i)=stdfsc(itype(i,mnum),mnum)*dsqrt(gamvec(ind))
6806 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcsc(i)
6811 end subroutine sdarea
6812 !-----------------------------------------------------------------------------
6814 !-----------------------------------------------------------------------------
6817 ! ###################################################
6818 ! ## COPYRIGHT (C) 1996 by Jay William Ponder ##
6819 ! ## All Rights Reserved ##
6820 ! ###################################################
6822 ! ################################################################
6824 ! ## subroutine surfatom -- exposed surface area of an atom ##
6826 ! ################################################################
6829 ! "surfatom" performs an analytical computation of the surface
6830 ! area of a specified atom; a simplified version of "surface"
6832 ! literature references:
6834 ! T. J. Richmond, "Solvent Accessible Surface Area and
6835 ! Excluded Volume in Proteins", Journal of Molecular Biology,
6838 ! L. Wesson and D. Eisenberg, "Atomic Solvation Parameters
6839 ! Applied to Molecular Dynamics of Proteins in Solution",
6840 ! Protein Science, 1, 227-235 (1992)
6842 ! variables and parameters:
6844 ! ir number of atom for which area is desired
6845 ! area accessible surface area of the atom
6846 ! radius radii of each of the individual atoms
6849 subroutine surfatom(ir,area,radius)
6851 ! implicit real*8 (a-h,o-z)
6852 ! include 'DIMENSIONS'
6854 ! include 'COMMON.GEO'
6855 ! include 'COMMON.IOUNITS'
6857 integer :: nsup,nstart_sup
6858 ! double precision c,dc,dc_old,d_c_work,xloc,xrot,dc_norm
6859 ! common /chain/ c(3,maxres2+2),dc(3,0:maxres2),dc_old(3,0:maxres2),
6860 ! & xloc(3,maxres),xrot(3,maxres),dc_norm(3,0:maxres2),
6861 ! & dc_work(MAXRES6),nres,nres0
6862 integer,parameter :: maxarc=300
6866 integer :: mi,ni,narc
6867 integer :: key(maxarc)
6868 integer :: intag(maxarc)
6869 integer :: intag1(maxarc)
6870 real(kind=8) :: area,arcsum
6871 real(kind=8) :: arclen,exang
6872 real(kind=8) :: delta,delta2
6873 real(kind=8) :: eps,rmove
6874 real(kind=8) :: xr,yr,zr
6875 real(kind=8) :: rr,rrsq
6876 real(kind=8) :: rplus,rminus
6877 real(kind=8) :: axx,axy,axz
6878 real(kind=8) :: ayx,ayy
6879 real(kind=8) :: azx,azy,azz
6880 real(kind=8) :: uxj,uyj,uzj
6881 real(kind=8) :: tx,ty,tz
6882 real(kind=8) :: txb,tyb,td
6883 real(kind=8) :: tr2,tr,txr,tyr
6884 real(kind=8) :: tk1,tk2
6885 real(kind=8) :: thec,the,t,tb
6886 real(kind=8) :: txk,tyk,tzk
6887 real(kind=8) :: t1,ti,tf,tt
6888 real(kind=8) :: txj,tyj,tzj
6889 real(kind=8) :: ccsq,cc,xysq
6890 real(kind=8) :: bsqk,bk,cosine
6891 real(kind=8) :: dsqj,gi,pix2
6892 real(kind=8) :: therk,dk,gk
6893 real(kind=8) :: risqk,rik
6894 real(kind=8) :: radius(2*nres) !(maxatm) (maxatm=maxres2)
6895 real(kind=8) :: ri(maxarc),risq(maxarc)
6896 real(kind=8) :: ux(maxarc),uy(maxarc),uz(maxarc)
6897 real(kind=8) :: xc(maxarc),yc(maxarc),zc(maxarc)
6898 real(kind=8) :: xc1(maxarc),yc1(maxarc),zc1(maxarc)
6899 real(kind=8) :: dsq(maxarc),bsq(maxarc)
6900 real(kind=8) :: dsq1(maxarc),bsq1(maxarc)
6901 real(kind=8) :: arci(maxarc),arcf(maxarc)
6902 real(kind=8) :: ex(maxarc),lt(maxarc),gr(maxarc)
6903 real(kind=8) :: b(maxarc),b1(maxarc),bg(maxarc)
6904 real(kind=8) :: kent(maxarc),kout(maxarc)
6905 real(kind=8) :: ther(maxarc)
6906 logical :: moved,top
6907 logical :: omit(maxarc)
6910 ! maxatm = 2*nres !maxres2 maxres2=2*maxres
6911 ! maxlight = 8*maxatm
6914 ! maxtors = 4*maxatm
6916 ! zero out the surface area for the sphere of interest
6919 ! write (2,*) "ir",ir," radius",radius(ir)
6920 if (radius(ir) .eq. 0.0d0) return
6922 ! set the overlap significance and connectivity shift
6926 delta2 = delta * delta
6931 ! store coordinates and radius of the sphere of interest
6939 ! initialize values of some counters and summations
6948 ! test each sphere to see if it overlaps the sphere of interest
6951 if (i.eq.ir .or. radius(i).eq.0.0d0) goto 30
6952 rplus = rr + radius(i)
6954 if (abs(tx) .ge. rplus) goto 30
6956 if (abs(ty) .ge. rplus) goto 30
6958 if (abs(tz) .ge. rplus) goto 30
6960 ! check for sphere overlap by testing distance against radii
6962 xysq = tx*tx + ty*ty
6963 if (xysq .lt. delta2) then
6970 if (rplus-cc .le. delta) goto 30
6971 rminus = rr - radius(i)
6973 ! check to see if sphere of interest is completely buried
6975 if (cc-abs(rminus) .le. delta) then
6976 if (rminus .le. 0.0d0) goto 170
6980 ! check for too many overlaps with sphere of interest
6982 if (io .ge. maxarc) then
6984 20 format (/,' SURFATOM -- Increase the Value of MAXARC')
6988 ! get overlap between current sphere and sphere of interest
6997 gr(io) = (ccsq+rplus*rminus) / (2.0d0*rr*b1(io))
7003 ! case where no other spheres overlap the sphere of interest
7006 area = 4.0d0 * pi * rrsq
7010 ! case where only one sphere overlaps the sphere of interest
7013 area = pix2 * (1.0d0 + gr(1))
7014 area = mod(area,4.0d0*pi) * rrsq
7018 ! case where many spheres intersect the sphere of interest;
7019 ! sort the intersecting spheres by their degree of overlap
7021 call sort2 (io,gr,key)
7024 intag(i) = intag1(k)
7033 ! get radius of each overlap circle on surface of the sphere
7038 risq(i) = rrsq - gi*gi
7039 ri(i) = sqrt(risq(i))
7040 ther(i) = 0.5d0*pi - asin(min(1.0d0,max(-1.0d0,gr(i))))
7043 ! find boundary of inaccessible area on sphere of interest
7046 if (.not. omit(k)) then
7053 ! check to see if J circle is intersecting K circle;
7054 ! get distance between circle centers and sum of radii
7057 if (omit(j)) goto 60
7058 cc = (txk*xc(j)+tyk*yc(j)+tzk*zc(j))/(bk*b(j))
7059 cc = acos(min(1.0d0,max(-1.0d0,cc)))
7060 td = therk + ther(j)
7062 ! check to see if circles enclose separate regions
7064 if (cc .ge. td) goto 60
7066 ! check for circle J completely inside circle K
7068 if (cc+ther(j) .lt. therk) goto 40
7070 ! check for circles that are essentially parallel
7072 if (cc .gt. delta) goto 50
7077 ! check to see if sphere of interest is completely buried
7080 if (pix2-cc .le. td) goto 170
7086 ! find T value of circle intersections
7089 if (omit(k)) goto 110
7104 ! rotation matrix elements
7116 if (.not. omit(j)) then
7121 ! rotate spheres so K vector colinear with z-axis
7123 uxj = txj*axx + tyj*axy - tzj*axz
7124 uyj = tyj*ayy - txj*ayx
7125 uzj = txj*azx + tyj*azy + tzj*azz
7126 cosine = min(1.0d0,max(-1.0d0,uzj/b(j)))
7127 if (acos(cosine) .lt. therk+ther(j)) then
7128 dsqj = uxj*uxj + uyj*uyj
7133 tr2 = risqk*dsqj - tb*tb
7139 ! get T values of intersection for K circle
7142 tb = min(1.0d0,max(-1.0d0,tb))
7144 if (tyb-txr .lt. 0.0d0) tk1 = pix2 - tk1
7146 tb = min(1.0d0,max(-1.0d0,tb))
7148 if (tyb+txr .lt. 0.0d0) tk2 = pix2 - tk2
7149 thec = (rrsq*uzj-gk*bg(j)) / (rik*ri(j)*b(j))
7150 if (abs(thec) .lt. 1.0d0) then
7152 else if (thec .ge. 1.0d0) then
7154 else if (thec .le. -1.0d0) then
7158 ! see if "tk1" is entry or exit point; check t=0 point;
7159 ! "ti" is exit point, "tf" is entry point
7161 cosine = min(1.0d0,max(-1.0d0, &
7162 (uzj*gk-uxj*rik)/(b(j)*rr)))
7163 if ((acos(cosine)-ther(j))*(tk2-tk1) .le. 0.0d0) then
7171 if (narc .ge. maxarc) then
7173 70 format (/,' SURFATOM -- Increase the Value',&
7177 if (tf .le. ti) then
7198 ! special case; K circle without intersections
7200 if (narc .le. 0) goto 90
7202 ! general case; sum up arclength and set connectivity code
7204 call sort2 (narc,arci,key)
7209 if (narc .gt. 1) then
7212 if (t .lt. arci(j)) then
7213 arcsum = arcsum + arci(j) - t
7214 exang = exang + ex(ni)
7216 if (jb .ge. maxarc) then
7218 80 format (/,' SURFATOM -- Increase the Value',&
7223 kent(jb) = maxarc*i + k
7225 kout(jb) = maxarc*k + i
7234 arcsum = arcsum + pix2 - t
7236 exang = exang + ex(ni)
7239 kent(jb) = maxarc*i + k
7241 kout(jb) = maxarc*k + i
7248 arclen = arclen + gr(k)*arcsum
7251 if (arclen .eq. 0.0d0) goto 170
7252 if (jb .eq. 0) goto 150
7254 ! find number of independent boundaries and check connectivity
7258 if (kout(k) .ne. 0) then
7265 if (m .eq. kent(ii)) then
7268 if (j .eq. jb) goto 150
7280 ! attempt to fix connectivity error by moving atom slightly
7284 140 format (/,' SURFATOM -- Connectivity Error at Atom',i6)
7293 ! compute the exposed surface area for the sphere of interest
7296 area = ib*pix2 + exang + arclen
7297 area = mod(area,4.0d0*pi) * rrsq
7299 ! attempt to fix negative area by moving atom slightly
7301 if (area .lt. 0.0d0) then
7304 160 format (/,' SURFATOM -- Negative Area at Atom',i6)
7315 end subroutine surfatom
7316 !----------------------------------------------------------------
7317 !----------------------------------------------------------------
7318 subroutine alloc_MD_arrays
7319 !EL Allocation of arrays used by MD module
7321 integer :: nres2,nres6
7324 !----------------------
7328 allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
7329 allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
7330 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
7331 allocate(fricvec(nres2,nres2))
7332 allocate(pfric_mat(nres2,nres2),vfric_mat(nres2,nres2))
7333 allocate(afric_mat(nres2,nres2),prand_mat(nres2,nres2))
7334 allocate(vrand_mat1(nres2,nres2),vrand_mat2(nres2,nres2)) !(MAXRES2,MAXRES2)
7335 allocate(pfric0_mat(nres2,nres2,0:maxflag_stoch))
7336 allocate(afric0_mat(nres2,nres2,0:maxflag_stoch))
7337 allocate(vfric0_mat(nres2,nres2,0:maxflag_stoch))
7338 allocate(prand0_mat(nres2,nres2,0:maxflag_stoch))
7339 allocate(vrand0_mat1(nres2,nres2,0:maxflag_stoch))
7340 allocate(vrand0_mat2(nres2,nres2,0:maxflag_stoch)) !(MAXRES2,MAXRES2,0:maxflag_stoch)
7341 allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
7343 allocate(mt1(nres2,nres2),mt2(nres2,nres2),mt3(nres2,nres2)) !(maxres2,maxres2)
7344 !----------------------
7346 ! commom.langevin.lang0
7348 allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
7350 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
7351 allocate(fricvec(nres2,nres2)) !(MAXRES2,MAXRES2)
7353 allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
7354 allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
7357 if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2))
7359 !----------------------
7360 ! commom.hairpin in CSA module
7361 !----------------------
7362 ! common.mce in MCM_MD module
7363 !----------------------
7365 ! common /mdgrad/ in module.energy
7366 ! common /back_constr/ in module.energy
7367 ! common /qmeas/ in module.energy
7370 allocate(potEcomp(0:n_ene+4)) !(0:n_ene+4)
7372 allocate(d_t(3,0:nres2),d_a(3,0:nres2),d_t_old(3,0:nres2)) !(3,0:MAXRES2)
7373 allocate(d_a_work(nres6)) !(6*MAXRES)
7375 allocate(DM(nres2),DU1(nres2),DU2(nres2))
7376 allocate(DMorig(nres2),DU1orig(nres2),DU2orig(nres2))
7378 allocate(Gvec(nres2,nres2))
7381 write (iout,*) "Before A Allocation",nfgtasks-1
7383 allocate(Gmat(nres2,nres2),A(nres2,nres2))
7384 if(.not.allocated(Ginv)) allocate(Ginv(nres2,nres2)) !in control: ergastulum
7385 allocate(Gsqrp(nres2,nres2),Gsqrm(nres2,nres2),Gvec(nres2,nres2)) !(maxres2,maxres2)
7387 allocate(Geigen(nres2)) !(maxres2)
7388 if(.not.allocated(vtot)) allocate(vtot(nres2)) !(maxres2)
7389 ! common /inertia/ in io_conf: parmread
7390 ! real(kind=8),dimension(:),allocatable :: ISC,msc !(ntyp+1)
7391 ! common /langevin/in io read_MDpar
7392 ! real(kind=8),dimension(:),allocatable :: gamsc !(ntyp1)
7393 ! real(kind=8),dimension(:),allocatable :: stdfsc !(ntyp)
7394 ! in io_conf: parmread
7395 ! real(kind=8),dimension(:),allocatable :: restok !(ntyp+1)
7396 ! common /mdpmpi/ in control: ergastulum
7397 if(.not.allocated(ng_start)) allocate(ng_start(0:nfgtasks-1))
7398 if(.not.allocated(ng_counts)) allocate(ng_counts(0:nfgtasks-1))
7399 if(.not.allocated(nginv_counts)) allocate(nginv_counts(0:nfgtasks-1)) !(0:MaxProcs-1)
7400 if(.not.allocated(nginv_start)) allocate(nginv_start(0:nfgtasks)) !(0:MaxProcs)
7401 !----------------------
7402 ! common.muca in read_muca
7403 ! common /double_muca/
7404 ! real(kind=8) :: elow,ehigh,factor,hbin,factor_min
7405 ! real(kind=8),dimension(:),allocatable :: emuca,nemuca,&
7406 ! nemuca2,hist !(4*maxres)
7407 ! common /integer_muca/
7408 ! integer :: nmuca,imtime,muca_smooth
7410 ! real(kind=8),dimension(:),allocatable :: elowi,ehighi !(maxprocs)
7411 !----------------------
7413 ! common /mdgrad/ in module.energy
7414 ! common /back_constr/ in module.energy
7415 ! common /qmeas/ in module.energy
7419 allocate(d_t_work(nres6),d_t_work_new(nres6),d_af_work(nres6))
7420 allocate(d_as_work(nres6),kinetic_force(nres6)) !(MAXRES6)
7421 allocate(d_t_new(3,0:nres2),d_a_old(3,0:nres2),d_a_short(3,0:nres2)) !,d_a !(3,0:MAXRES2)
7422 allocate(stdforcp(nres),stdforcsc(nres)) !(MAXRES)
7423 !----------------------
7425 allocate(D_ban(nres6)) !(MAXRES6) maxres6=6*maxres
7426 ! common /stochcalc/ stochforcvec
7427 allocate(stochforcvec(nres6)) !(MAXRES6) maxres6=6*maxres
7428 !----------------------
7430 end subroutine alloc_MD_arrays
7431 !-----------------------------------------------------------------------------
7432 !-----------------------------------------------------------------------------