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
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)=msc(itype(i,mnum),mnum)
680 !c write (iout,*) "Kinetic trp:",i,(v(j),j=1,3)
681 vtot(i)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
682 KEt_p=KEt_p+mp(mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
685 incr(j)=incr(j)+d_t(j,i)
688 !c write(iout,*) 'KEt_p', KEt_p
689 !c The translational part for the side chain virtual bond
690 !c Only now we can initialize incr with zeros. It must be equal
691 !c to the velocities of the first Calpha.
697 iti=iabs(itype(i,mnum))
698 if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
705 v(j)=incr(j)+d_t(j,nres+i)
708 write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3)
709 write (iout,*) "i",i," msc",msc(iti,1)," v",(v(j),j=1,3)
710 KEt_sc=KEt_sc+msc(iti,mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
711 vtot(i+nres)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
713 incr(j)=incr(j)+d_t(j,i)
717 ! write(iout,*) 'KEt_sc', KEt_sc
718 ! The part due to stretching and rotation of the peptide groups
721 if (itype(i,mnum).ne.ntyp1_molec(mnum)&
722 .and.itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
723 if (mnum.eq.5) Ip(mnum)=Isc(itype(i,mnum),mnum)
724 ! write (iout,*) "i",i
725 ! write (iout,*) "i",i," mag1",mag1," mag2",mag2
729 !c write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3)
730 KEr_p=KEr_p+Ip(mnum)*(incr(1)*incr(1)+incr(2)*incr(2) &
735 !c write(iout,*) 'KEr_p', KEr_p
736 !c The rotational part of the side chain virtual bond
739 iti=iabs(itype(i,mnum))
740 if (itype(i,1).ne.10.and.itype(i,mnum).ne.ntyp1_molec(mnum)&
743 incr(j)=d_t(j,nres+i)
745 ! write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3)
746 KEr_sc=KEr_sc+Isc(iti,mnum)*(incr(1)*incr(1)+incr(2)*incr(2)+&
750 !c The total kinetic energy
752 !c write(iout,*) ' KEt_p',KEt_p,' KEt_sc',KEt_sc,' KEr_p',KEr_p,
753 !c & ' KEr_sc', KEr_sc
754 KE_total=0.5d0*(KEt_p+KEt_sc+0.25d0*KEr_p+KEr_sc)
755 !c write (iout,*) "KE_total",KE_total
757 end subroutine kinetic
760 !-----------------------------------------------------------------------------
761 subroutine kinetic(KE_total)
762 !----------------------------------------------------------------
763 ! This subroutine calculates the total kinetic energy of the chain
764 !-----------------------------------------------------------------
766 ! implicit real*8 (a-h,o-z)
767 ! include 'DIMENSIONS'
768 ! include 'COMMON.VAR'
769 ! include 'COMMON.CHAIN'
770 ! include 'COMMON.DERIV'
771 ! include 'COMMON.GEO'
772 ! include 'COMMON.LOCAL'
773 ! include 'COMMON.INTERACT'
774 ! include 'COMMON.MD'
775 ! include 'COMMON.IOUNITS'
776 real(kind=8) :: KE_total,mscab
778 integer :: i,j,k,iti,mnum,term
779 real(kind=8) :: KEt_p,KEt_sc,KEr_p,KEr_sc,incr(3),&
782 write (iout,*) "Velocities, kietic"
784 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
785 (d_t(j,i+nres),j=1,3)
790 ! write (iout,*) "ISC",(isc(itype(i,1)),i=1,nres)
791 ! The translational part for peptide virtual bonds
796 ! if (molnum(nct).gt.3) term=nct
799 if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
800 ! write (iout,*) "Kinetic trp:",i,(incr(j),j=1,3),mp(mnum)
803 v(j)=incr(j)+0.5d0*d_t(j,i)
807 v(j)=incr(j)+0.5d0*d_t(j,i)
810 vtot(i)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
811 KEt_p=KEt_p+mp(mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
813 incr(j)=incr(j)+d_t(j,i)
816 ! write(iout,*) 'KEt_p', KEt_p
817 ! The translational part for the side chain virtual bond
818 ! Only now we can initialize incr with zeros. It must be equal
819 ! to the velocities of the first Calpha.
825 iti=iabs(itype(i,mnum))
826 ! if (mnum.ge.4) then
831 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
832 if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
833 .or.(mnum.ge.4)) then
839 v(j)=incr(j)+d_t(j,nres+i)
842 ! write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3)
843 ! write (iout,*) "i",i," msc",msc(iti,mnum)," v",(v(j),j=1,3)
844 KEt_sc=KEt_sc+mscab*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
845 vtot(i+nres)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
847 incr(j)=incr(j)+d_t(j,i)
851 ! write(iout,*) 'KEt_sc', KEt_sc
852 ! The part due to stretching and rotation of the peptide groups
856 ! write (iout,*) "i",i
857 ! write (iout,*) "i",i," mag1",mag1," mag2",mag2
861 ! write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3)
862 KEr_p=KEr_p+Ip(mnum)*(incr(1)*incr(1)+incr(2)*incr(2) &
866 ! write(iout,*) 'KEr_p', KEr_p
867 ! The rotational part of the side chain virtual bond
871 iti=iabs(itype(i,mnum))
872 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
873 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
874 .and.(mnum.lt.4)) then
876 incr(j)=d_t(j,nres+i)
878 ! write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3)
879 KEr_sc=KEr_sc+Isc(iti,mnum)*(incr(1)*incr(1)+incr(2)*incr(2)+ &
883 ! The total kinetic energy
885 ! write(iout,*) 'KEr_sc', KEr_sc
886 KE_total=0.5d0*(KEt_p+KEt_sc+0.25d0*KEr_p+KEr_sc)
887 ! write (iout,*) "KE_total",KE_total
889 end subroutine kinetic
890 !-----------------------------------------------------------------------------
892 subroutine kinetic_CASC(KE_total)
893 !c----------------------------------------------------------------
894 !c Compute the kinetic energy of the system using the Calpha-SC
896 !c-----------------------------------------------------------------
898 double precision KE_total
900 integer i,j,k,iti,ichain,innt,inct,mnum
901 double precision KEt_p,KEt_sc,KEr_p,KEr_sc,incr(3),&
908 !c write (iout,*) "ISC",(isc(itype(i)),i=1,nres)
909 !c The translational part for peptide virtual bonds
912 innt=chain_border(1,ichain)
913 inct=chain_border(2,ichain)
914 !c write (iout,*) "Kinetic_CASC chain",ichain," innt",innt,
919 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
920 !c write (iout,*) i,(d_t(j,i),j=1,3),(d_t(j,i+1),j=1,3)
922 v(j)=0.5d0*(d_t(j,i)+d_t(j,i+1))
924 !c write (iout,*) "Kinetic trp i",i," v",(v(j),j=1,3)
925 KEt_p=KEt_p+mp(mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
927 !c write(iout,*) 'KEt_p', KEt_p
928 !c The translational part for the side chain virtual bond
929 !c Only now we can initialize incr with zeros. It must be equal
930 !c to the velocities of the first Calpha.
933 iti=iabs(itype(i,mnum))
934 if (itype(i,1).eq.10.or.mnum.ge.3.or. itype(i,mnum).eq.ntyp1_molec(mnum)) then
935 !c write (iout,*) i,iti,(d_t(j,i),j=1,3)
940 !c write (iout,*) i,iti,(d_t(j,nres+i),j=1,3)
945 !c write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3)
946 !c write (iout,*) "i",i," msc",msc(iti)," v",(v(j),j=1,3)
947 KEt_sc=KEt_sc+msc(iti,mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
950 !c write(iout,*) 'KEt_sc', KEt_sc
951 !c The part due to stretching and rotation of the peptide groups
955 incr(j)=d_t(j,i+1)-d_t(j,i)
957 if (mnum.eq.5) Ip(mnum)=Isc(itype(i,mnum),mnum)
958 !c write (iout,*) i,(incr(j),j=1,3)
959 !c write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3)
960 KEr_p=KEr_p+Ip(mnum)*(incr(1)*incr(1)+incr(2)*incr(2)&
964 !c write(iout,*) 'KEr_p', KEr_p
965 !c The rotational part of the side chain virtual bond
968 iti=iabs(itype(i,mnum))
969 ! if (iti.ne.10.and.mnum.lt.3) then
970 if (itype(i,1).ne.10.and.mnum.lt.3.and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
972 incr(j)=d_t(j,nres+i)-d_t(j,i)
974 !c write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3)
975 KEr_sc=KEr_sc+Isc(iti,mnum)*(incr(1)*incr(1)+incr(2)*incr(2)+&
981 !c The total kinetic energy
983 !c write(iout,*) ' KEt_p',KEt_p,' KEt_sc',KEt_sc,' KEr_p',KEr_p,
984 !c & ' KEr_sc', KEr_sc
985 KE_total=0.5d0*(KEt_p+KEt_sc+0.25d0*KEr_p+KEr_sc)
986 !c write (iout,*) "KE_total",KE_tota
988 write (iout,*) "Need to compile with -DFIVEDIAG to use this sub!"
992 end subroutine kinetic_CASC
995 !-----------------------------------------------------------------------------
997 !------------------------------------------------
998 ! The driver for molecular dynamics subroutines
999 !------------------------------------------------
1002 use control, only:tcpu,ovrtim
1003 ! use io_comm, only:ilen
1005 use compare, only:secondary2,hairpin
1006 use io, only:cartout,statout
1007 ! implicit real*8 (a-h,o-z)
1008 ! include 'DIMENSIONS'
1011 integer :: IERROR,ERRCODE
1013 ! include 'COMMON.SETUP'
1014 ! include 'COMMON.CONTROL'
1015 ! include 'COMMON.VAR'
1016 ! include 'COMMON.MD'
1018 ! include 'COMMON.LANGEVIN'
1020 ! include 'COMMON.LANGEVIN.lang0'
1022 ! include 'COMMON.CHAIN'
1023 ! include 'COMMON.DERIV'
1024 ! include 'COMMON.GEO'
1025 ! include 'COMMON.LOCAL'
1026 ! include 'COMMON.INTERACT'
1027 ! include 'COMMON.IOUNITS'
1028 ! include 'COMMON.NAMES'
1029 ! include 'COMMON.TIME1'
1030 ! include 'COMMON.HAIRPIN'
1031 real(kind=8),dimension(3) :: L,vcm
1033 real(kind=8),dimension(6*nres) :: v_work,v_transf !(maxres6) maxres6=6*maxres
1035 integer :: rstcount !ilen,
1037 character(len=50) :: tytul
1038 !el common /gucio/ cm
1039 integer :: i,j,nharp
1040 integer,dimension(4,nres) :: iharp !(4,nres/3)(4,maxres/3)
1042 real(kind=8) :: tt0,scalfac
1043 integer :: nres2,itime
1048 print *,"MY tmpdir",tmpdir,ilen(tmpdir)
1049 if (ilen(tmpdir).gt.0) &
1050 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_" &
1051 //liczba(:ilen(liczba))//'.rst')
1053 if (ilen(tmpdir).gt.0) &
1054 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"//'.rst')
1061 write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
1067 print *,"just befor setup matix",nres
1068 ! Determine the inverse of the inertia matrix.
1069 call setup_MD_matrices
1071 print *,"AFTER SETUP MATRICES"
1073 print *,"AFTER INIT MD"
1076 t_MDsetup = MPI_Wtime()-tt0
1078 t_MDsetup = tcpu()-tt0
1081 ! Entering the MD loop
1087 if (lang.eq.2 .or. lang.eq.3) then
1091 call sd_verlet_p_setup
1093 call sd_verlet_ciccotti_setup
1097 pfric0_mat(i,j,0)=pfric_mat(i,j)
1098 afric0_mat(i,j,0)=afric_mat(i,j)
1099 vfric0_mat(i,j,0)=vfric_mat(i,j)
1100 prand0_mat(i,j,0)=prand_mat(i,j)
1101 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
1102 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
1105 flag_stoch(0)=.true.
1106 do i=1,maxflag_stoch
1107 flag_stoch(i)=.false.
1111 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
1113 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1117 else if (lang.eq.1 .or. lang.eq.4) then
1118 print *,"before setup_fricmat"
1120 print *,"after setup_fricmat"
1123 t_langsetup=MPI_Wtime()-tt0
1126 t_langsetup=tcpu()-tt0
1129 do itime=1,n_timestep
1130 if (large) print *,itime,ntwe
1132 if (large.and. mod(itime,ntwe).eq.0) &
1133 write (iout,*) "itime",itime
1135 if (lang.gt.0 .and. surfarea .and. &
1136 mod(itime,reset_fricmat).eq.0) then
1137 if (lang.eq.2 .or. lang.eq.3) then
1141 call sd_verlet_p_setup
1143 call sd_verlet_ciccotti_setup
1147 pfric0_mat(i,j,0)=pfric_mat(i,j)
1148 afric0_mat(i,j,0)=afric_mat(i,j)
1149 vfric0_mat(i,j,0)=vfric_mat(i,j)
1150 prand0_mat(i,j,0)=prand_mat(i,j)
1151 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
1152 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
1155 flag_stoch(0)=.true.
1156 do i=1,maxflag_stoch
1157 flag_stoch(i)=.false.
1160 else if (lang.eq.1 .or. lang.eq.4) then
1161 print *,"before setup_fricmat"
1163 print *,"after setup_fricmat"
1165 write (iout,'(a,i10)') &
1166 "Friction matrix reset based on surface area, itime",itime
1168 if (reset_vel .and. tbf .and. lang.eq.0 &
1169 .and. mod(itime,count_reset_vel).eq.0) then
1171 write(iout,'(a,f20.2)') &
1172 "Velocities reset to random values, time",totT
1175 d_t_old(j,i)=d_t(j,i)
1179 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
1183 d_t(j,0)=d_t(j,0)-vcm(j)
1186 kinetic_T=2.0d0/(dimen3*Rb)*EK
1187 scalfac=dsqrt(T_bath/kinetic_T)
1188 write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
1191 d_t_old(j,i)=scalfac*d_t(j,i)
1197 ! Time-reversible RESPA algorithm
1198 ! (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
1199 call RESPA_step(itime)
1201 ! Variable time step algorithm.
1202 if (large) print *,"before verlet_step"
1203 call velverlet_step(itime)
1204 if (large) print *,"after verlet_step"
1208 call brown_step(itime)
1210 print *,"Brown dynamics not here!"
1212 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1219 if (mod(itime,ntwe).eq.0) then
1223 ! call check_ecartint
1233 v_work(ind)=d_t(j,i)
1238 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum).and.mnum.lt.4) then
1241 v_work(ind)=d_t(j,i+nres)
1246 write (66,'(80f10.5)') &
1247 ((d_t(j,i),j=1,3),i=0,nres-1),((d_t(j,i+nres),j=1,3),i=1,nres)
1251 v_transf(i)=v_transf(i)+gvec(j,i)*v_work(j)
1253 v_transf(i)= v_transf(i)*dsqrt(geigen(i))
1255 write (67,'(80f10.5)') (v_transf(i),i=1,ind)
1258 if (mod(itime,ntwx).eq.0) then
1260 write (tytul,'("time",f8.2)') totT
1262 call hairpin(.true.,nharp,iharp)
1263 call secondary2(.true.)
1264 call pdbout(potE,tytul,ipdb)
1265 call enerprint(potEcomp)
1270 write(iout,*) "starting fodstep"
1271 call fodstep(nfodstep)
1272 write(iout,*) "after fodstep"
1275 call hairpin(.true.,nharp,iharp)
1276 call secondary2(.true.)
1277 call pdbout(potE,tytul,ipdb)
1284 if (rstcount.eq.1000.or.itime.eq.n_timestep) then
1285 open(irest2,file=rest2name,status='unknown')
1286 write(irest2,*) totT,EK,potE,totE,t_bath
1288 ! AL 4/17/17: Now writing d_t(0,:) too
1290 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
1292 ! AL 4/17/17: Now writing d_c(0,:) too
1294 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
1302 t_MD=MPI_Wtime()-tt0
1306 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))') &
1308 'MD calculations setup:',t_MDsetup,&
1309 'Energy & gradient evaluation:',t_enegrad,&
1310 'Stochastic MD setup:',t_langsetup,&
1311 'Stochastic MD step setup:',t_sdsetup,&
1313 write (iout,'(/28(1h=),a25,27(1h=))') &
1314 ' End of MD calculation '
1316 write (iout,*) "time for etotal",t_etotal," elong",t_elong,&
1318 write (iout,*) "time_fric",time_fric," time_stoch",time_stoch,&
1319 " time_fricmatmult",time_fricmatmult," time_fsample ",&
1324 !-----------------------------------------------------------------------------
1325 subroutine velverlet_step(itime)
1326 !-------------------------------------------------------------------------------
1327 ! Perform a single velocity Verlet step; the time step can be rescaled if
1328 ! increments in accelerations exceed the threshold
1329 !-------------------------------------------------------------------------------
1330 ! implicit real*8 (a-h,o-z)
1331 ! include 'DIMENSIONS'
1333 use control, only:tcpu
1335 use minimm, only:minim_dc
1338 integer :: ierror,ierrcode
1339 real(kind=8) :: errcode
1341 ! include 'COMMON.SETUP'
1342 ! include 'COMMON.VAR'
1343 ! include 'COMMON.MD'
1345 ! include 'COMMON.LANGEVIN'
1347 ! include 'COMMON.LANGEVIN.lang0'
1349 ! include 'COMMON.CHAIN'
1350 ! include 'COMMON.DERIV'
1351 ! include 'COMMON.GEO'
1352 ! include 'COMMON.LOCAL'
1353 ! include 'COMMON.INTERACT'
1354 ! include 'COMMON.IOUNITS'
1355 ! include 'COMMON.NAMES'
1356 ! include 'COMMON.TIME1'
1357 ! include 'COMMON.MUCA'
1358 real(kind=8),dimension(3) :: vcm,incr
1359 real(kind=8),dimension(3) :: L
1360 integer :: count,rstcount !ilen,
1362 character(len=50) :: tytul
1363 integer :: maxcount_scale = 30
1364 !el common /gucio/ cm
1365 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1366 !el common /stochcalc/ stochforcvec
1367 integer :: icount_scale,itime_scal,i,j,ifac_time,iretcode,itime
1369 real(kind=8) :: epdrift,tt0,fac_time
1371 if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres)) !(MAXRES6) maxres6=6*maxres
1377 if (large) print *,"after sddir_precalc"
1378 else if (lang.eq.2 .or. lang.eq.3) then
1380 call stochastic_force(stochforcvec)
1383 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
1385 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1392 icount_scale=icount_scale+1
1393 ! write(iout,*) "icount_scale",icount_scale,scalel
1394 if (icount_scale.gt.maxcount_scale) then
1396 "ERROR: too many attempts at scaling down the time step. ",&
1397 "amax=",amax,"epdrift=",epdrift,&
1398 "damax=",damax,"edriftmax=",edriftmax,&
1402 call MPI_Abort(MPI_COMM_WORLD,IERROR,IERRCODE)
1406 ! First step of the velocity Verlet algorithm
1411 else if (lang.eq.3) then
1413 call sd_verlet1_ciccotti
1415 else if (lang.eq.1) then
1420 ! Build the chain from the newly calculated coordinates
1421 call chainbuild_cart
1422 if (rattle) call rattle1
1424 if (large) then !.and. mod(itime,ntwe).eq.0) then
1425 write (iout,*) "Cartesian and internal coordinates: step 1"
1430 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),&
1431 (dc(j,i+nres),j=1,3)
1433 write (iout,*) "Accelerations"
1435 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1436 (d_a(j,i+nres),j=1,3)
1438 write (iout,*) "Velocities, step 1"
1440 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1441 (d_t(j,i+nres),j=1,3)
1450 ! Calculate energy and forces
1452 call etotal(potEcomp)
1453 ! AL 4/17/17: Reduce the steps if NaNs occurred.
1454 if (potEcomp(0).gt.0.99e18 .or. isnan(potEcomp(0)).gt.0) then
1455 call enerprint(potEcomp)
1457 if (icount_scale.gt.15) then
1458 write (iout,*) "Tu jest problem",potEcomp(0),d_time
1459 ! call gen_rand_conf(1,*335)
1460 ! call minim_dc(potEcomp(0),iretcode,100)
1463 ! call etotal(potEcomp)
1464 ! write(iout,*) "needed to repara,",potEcomp
1467 ! 335 write(iout,*) "Failed genrand"
1471 if (large.and. mod(itime,ntwe).eq.0) &
1472 call enerprint(potEcomp)
1475 t_etotal=t_etotal+MPI_Wtime()-tt0
1477 t_etotal=t_etotal+tcpu()-tt0
1480 potE=potEcomp(0)-potEcomp(51)
1482 ! Get the new accelerations
1485 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1487 t_enegrad=t_enegrad+tcpu()-tt0
1489 ! Determine maximum acceleration and scale down the timestep if needed
1491 amax=amax/(itime_scal+1)**2
1492 call predict_edrift(epdrift)
1493 ! write(iout,*) "amax=",amax,damax,epdrift,edriftmax,amax/(itime_scal+1)
1495 ! write (iout,*) "before enter if",scalel,icount_scale
1496 if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then
1497 ! write(iout,*) "I enter if"
1498 ! Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
1500 ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax)) &
1502 itime_scal=itime_scal+ifac_time
1503 ! fac_time=dmin1(damax/amax,0.5d0)
1504 fac_time=0.5d0**ifac_time
1505 d_time=d_time*fac_time
1506 if (lang.eq.2 .or. lang.eq.3) then
1508 ! write (iout,*) "Calling sd_verlet_setup: 1"
1509 ! Rescale the stochastic forces and recalculate or restore
1510 ! the matrices of tinker integrator
1511 if (itime_scal.gt.maxflag_stoch) then
1512 if (large) write (iout,'(a,i5,a)') &
1513 "Calculate matrices for stochastic step;",&
1514 " itime_scal ",itime_scal
1516 call sd_verlet_p_setup
1518 call sd_verlet_ciccotti_setup
1520 write (iout,'(2a,i3,a,i3,1h.)') &
1521 "Warning: cannot store matrices for stochastic",&
1522 " integration because the index",itime_scal,&
1523 " is greater than",maxflag_stoch
1524 write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",&
1525 " integration Langevin algorithm for better efficiency."
1526 else if (flag_stoch(itime_scal)) then
1527 if (large) write (iout,'(a,i5,a,l1)') &
1528 "Restore matrices for stochastic step; itime_scal ",&
1529 itime_scal," flag ",flag_stoch(itime_scal)
1532 pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
1533 afric_mat(i,j)=afric0_mat(i,j,itime_scal)
1534 vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
1535 prand_mat(i,j)=prand0_mat(i,j,itime_scal)
1536 vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
1537 vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
1541 if (large) write (iout,'(2a,i5,a,l1)') &
1542 "Calculate & store matrices for stochastic step;",&
1543 " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
1545 call sd_verlet_p_setup
1547 call sd_verlet_ciccotti_setup
1549 flag_stoch(ifac_time)=.true.
1552 pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
1553 afric0_mat(i,j,itime_scal)=afric_mat(i,j)
1554 vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
1555 prand0_mat(i,j,itime_scal)=prand_mat(i,j)
1556 vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
1557 vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
1561 fac_time=1.0d0/dsqrt(fac_time)
1563 stochforcvec(i)=fac_time*stochforcvec(i)
1566 else if (lang.eq.1) then
1567 ! Rescale the accelerations due to stochastic forces
1568 fac_time=1.0d0/dsqrt(fac_time)
1570 d_as_work(i)=d_as_work(i)*fac_time
1573 if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)') &
1574 "itime",itime," Timestep scaled down to ",&
1575 d_time," ifac_time",ifac_time," itime_scal",itime_scal
1577 ! Second step of the velocity Verlet algorithm
1582 else if (lang.eq.3) then
1584 call sd_verlet2_ciccotti
1586 else if (lang.eq.1) then
1591 if (rattle) call rattle2
1594 if (d_time.ne.d_time0) then
1597 if (lang.eq.2 .or. lang.eq.3) then
1598 if (large) write (iout,'(a)') &
1599 "Restore original matrices for stochastic step"
1600 ! write (iout,*) "Calling sd_verlet_setup: 2"
1601 ! Restore the matrices of tinker integrator if the time step has been restored
1604 pfric_mat(i,j)=pfric0_mat(i,j,0)
1605 afric_mat(i,j)=afric0_mat(i,j,0)
1606 vfric_mat(i,j)=vfric0_mat(i,j,0)
1607 prand_mat(i,j)=prand0_mat(i,j,0)
1608 vrand_mat1(i,j)=vrand0_mat1(i,j,0)
1609 vrand_mat2(i,j)=vrand0_mat2(i,j,0)
1617 ! Calculate the kinetic and the total energy and the kinetic temperature
1621 ! call kinetic1(EK1)
1622 ! write (iout,*) "step",itime," EK",EK," EK1",EK1
1624 ! Couple the system to Berendsen bath if needed
1625 if (tbf .and. lang.eq.0) then
1628 kinetic_T=2.0d0/(dimen3*Rb)*EK
1629 ! Backup the coordinates, velocities, and accelerations
1633 d_t_old(j,i)=d_t(j,i)
1634 d_a_old(j,i)=d_a(j,i)
1638 if (mod(itime,ntwe).eq.0 .and. large) then
1639 write (iout,*) "Velocities, step 2"
1641 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1642 (d_t(j,i+nres),j=1,3)
1647 end subroutine velverlet_step
1648 !-----------------------------------------------------------------------------
1649 subroutine RESPA_step(itime)
1650 !-------------------------------------------------------------------------------
1651 ! Perform a single RESPA step.
1652 !-------------------------------------------------------------------------------
1653 ! implicit real*8 (a-h,o-z)
1654 ! include 'DIMENSIONS'
1658 use control, only:tcpu
1660 ! use io_conf, only:cartprint
1663 integer :: IERROR,ERRCODE
1665 ! include 'COMMON.SETUP'
1666 ! include 'COMMON.CONTROL'
1667 ! include 'COMMON.VAR'
1668 ! include 'COMMON.MD'
1670 ! include 'COMMON.LANGEVIN'
1672 ! include 'COMMON.LANGEVIN.lang0'
1674 ! include 'COMMON.CHAIN'
1675 ! include 'COMMON.DERIV'
1676 ! include 'COMMON.GEO'
1677 ! include 'COMMON.LOCAL'
1678 ! include 'COMMON.INTERACT'
1679 ! include 'COMMON.IOUNITS'
1680 ! include 'COMMON.NAMES'
1681 ! include 'COMMON.TIME1'
1682 real(kind=8),dimension(0:n_ene) :: energia_short,energia_long
1683 real(kind=8),dimension(3) :: L,vcm,incr
1684 real(kind=8),dimension(3,0:2*nres) :: dc_old0,d_t_old0,d_a_old0 !(3,0:maxres2) maxres2=2*maxres
1685 logical :: PRINT_AMTS_MSG = .false.
1686 integer :: count,rstcount !ilen,
1688 character(len=50) :: tytul
1689 integer :: maxcount_scale = 10
1690 !el common /gucio/ cm
1691 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1692 !el common /stochcalc/ stochforcvec
1693 integer :: itt,i,j,itsplit,itime
1695 !el common /cipiszcze/ itt
1697 real(kind=8) :: epdrift,tt0,epdriftmax
1700 if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres)) !(MAXRES6) maxres6=6*maxres
1704 if (large.and. mod(itime,ntwe).eq.0) then
1705 write (iout,*) "***************** RESPA itime",itime
1706 write (iout,*) "Cartesian and internal coordinates: step 0"
1708 call pdbout(0.0d0,"cipiszcze",iout)
1710 write (iout,*) "Accelerations from long-range forces"
1712 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1713 (d_a(j,i+nres),j=1,3)
1715 write (iout,*) "Velocities, step 0"
1717 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1718 (d_t(j,i+nres),j=1,3)
1723 ! Perform the initial RESPA step (increment velocities)
1724 ! write (iout,*) "*********************** RESPA ini"
1727 if (mod(itime,ntwe).eq.0 .and. large) then
1728 write (iout,*) "Velocities, end"
1730 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1731 (d_t(j,i+nres),j=1,3)
1735 ! Compute the short-range forces
1741 ! 7/2/2009 commented out
1743 ! call etotal_short(energia_short)
1746 ! 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
1749 d_a(j,i)=d_a_short(j,i)
1753 if (large.and. mod(itime,ntwe).eq.0) then
1754 write (iout,*) "energia_short",energia_short(0)
1755 write (iout,*) "Accelerations from short-range forces"
1757 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1758 (d_a(j,i+nres),j=1,3)
1763 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1765 t_enegrad=t_enegrad+tcpu()-tt0
1770 d_t_old(j,i)=d_t(j,i)
1771 d_a_old(j,i)=d_a(j,i)
1774 ! 6/30/08 A-MTS: attempt at increasing the split number
1777 dc_old0(j,i)=dc_old(j,i)
1778 d_t_old0(j,i)=d_t_old(j,i)
1779 d_a_old0(j,i)=d_a_old(j,i)
1782 if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
1783 if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
1790 ! write (iout,*) "itime",itime," ntime_split",ntime_split
1791 ! Split the time step
1792 d_time=d_time0/ntime_split
1793 ! Perform the short-range RESPA steps (velocity Verlet increments of
1794 ! positions and velocities using short-range forces)
1795 ! write (iout,*) "*********************** RESPA split"
1796 do itsplit=1,ntime_split
1799 else if (lang.eq.2 .or. lang.eq.3) then
1801 call stochastic_force(stochforcvec)
1804 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
1806 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1811 ! First step of the velocity Verlet algorithm
1816 else if (lang.eq.3) then
1818 call sd_verlet1_ciccotti
1820 else if (lang.eq.1) then
1825 ! Build the chain from the newly calculated coordinates
1826 call chainbuild_cart
1827 if (rattle) call rattle1
1829 if (large.and. mod(itime,ntwe).eq.0) then
1830 write (iout,*) "***** ITSPLIT",itsplit
1831 write (iout,*) "Cartesian and internal coordinates: step 1"
1832 call pdbout(0.0d0,"cipiszcze",iout)
1835 write (iout,*) "Velocities, step 1"
1837 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1838 (d_t(j,i+nres),j=1,3)
1847 ! Calculate energy and forces
1849 call etotal_short(energia_short)
1850 ! AL 4/17/17: Exit itime_split loop when energy goes infinite
1851 if (energia_short(0).gt.0.99e20 .or. isnan(energia_short(0)) ) then
1852 if (PRINT_AMTS_MSG) &
1853 write (iout,*) "Infinities/NaNs in energia_short",energia_short(0),"; increasing ntime_split to",ntime_split
1854 ntime_split=ntime_split*2
1855 if (ntime_split.gt.maxtime_split) then
1858 "Cannot rescue the run; aborting job. Retry with a smaller time step"
1860 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1863 "Cannot rescue the run; terminating. Retry with a smaller time step"
1869 if (large.and. mod(itime,ntwe).eq.0) &
1870 call enerprint(energia_short)
1873 t_eshort=t_eshort+MPI_Wtime()-tt0
1875 t_eshort=t_eshort+tcpu()-tt0
1879 ! Get the new accelerations
1881 ! 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
1884 d_a_short(j,i)=d_a(j,i)
1888 if (large.and. mod(itime,ntwe).eq.0) then
1889 write (iout,*)"energia_short",energia_short(0)
1890 write (iout,*) "Accelerations from short-range forces"
1892 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1893 (d_a(j,i+nres),j=1,3)
1898 ! Determine maximum acceleration and scale down the timestep if needed
1900 amax=amax/ntime_split**2
1901 call predict_edrift(epdrift)
1902 if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0) &
1903 write (iout,*) "amax",amax," damax",damax,&
1904 " epdrift",epdrift," epdriftmax",epdriftmax
1905 ! Exit loop and try with increased split number if the change of
1906 ! acceleration is too big
1907 if (amax.gt.damax .or. epdrift.gt.edriftmax) then
1908 if (ntime_split.lt.maxtime_split) then
1910 ntime_split=ntime_split*2
1911 ! AL 4/17/17: We should exit the itime_split loop when acceleration change is too big
1915 dc_old(j,i)=dc_old0(j,i)
1916 d_t_old(j,i)=d_t_old0(j,i)
1917 d_a_old(j,i)=d_a_old0(j,i)
1920 if (PRINT_AMTS_MSG) then
1921 write (iout,*) "acceleration/energy drift too large",amax,&
1922 epdrift," split increased to ",ntime_split," itime",itime,&
1928 "Uh-hu. Bumpy landscape. Maximum splitting number",&
1930 " already reached!!! Trying to carry on!"
1934 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1936 t_enegrad=t_enegrad+tcpu()-tt0
1938 ! Second step of the velocity Verlet algorithm
1943 else if (lang.eq.3) then
1945 call sd_verlet2_ciccotti
1947 else if (lang.eq.1) then
1952 if (rattle) call rattle2
1953 ! Backup the coordinates, velocities, and accelerations
1957 d_t_old(j,i)=d_t(j,i)
1958 d_a_old(j,i)=d_a(j,i)
1965 ! Restore the time step
1967 ! Compute long-range forces
1974 call etotal_long(energia_long)
1975 if (energia_long(0).gt.0.99e20 .or. isnan(energia_long(0))) then
1978 "Infinitied/NaNs in energia_long, Aborting MPI job."
1980 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1982 write (iout,*) "Infinitied/NaNs in energia_long, terminating."
1986 if (large.and. mod(itime,ntwe).eq.0) &
1987 call enerprint(energia_long)
1990 t_elong=t_elong+MPI_Wtime()-tt0
1992 t_elong=t_elong+tcpu()-tt0
1995 potE=potEcomp(0)-potEcomp(51)
1999 t_enegrad=t_enegrad+MPI_Wtime()-tt0
2001 t_enegrad=t_enegrad+tcpu()-tt0
2003 ! Compute accelerations from long-range forces
2005 if (large.and. mod(itime,ntwe).eq.0) then
2006 write (iout,*) "energia_long",energia_long(0)
2007 write (iout,*) "Cartesian and internal coordinates: step 2"
2009 call pdbout(0.0d0,"cipiszcze",iout)
2011 write (iout,*) "Accelerations from long-range forces"
2013 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
2014 (d_a(j,i+nres),j=1,3)
2016 write (iout,*) "Velocities, step 2"
2018 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2019 (d_t(j,i+nres),j=1,3)
2023 ! Compute the final RESPA step (increment velocities)
2024 ! write (iout,*) "*********************** RESPA fin"
2026 ! Compute the complete potential energy
2028 potEcomp(i)=energia_short(i)+energia_long(i)
2030 potE=potEcomp(0)-potEcomp(51)
2031 ! potE=energia_short(0)+energia_long(0)
2034 ! Calculate the kinetic and the total energy and the kinetic temperature
2037 ! Couple the system to Berendsen bath if needed
2038 if (tbf .and. lang.eq.0) then
2041 kinetic_T=2.0d0/(dimen3*Rb)*EK
2042 ! Backup the coordinates, velocities, and accelerations
2044 if (mod(itime,ntwe).eq.0 .and. large) then
2045 write (iout,*) "Velocities, end"
2047 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2048 (d_t(j,i+nres),j=1,3)
2053 end subroutine RESPA_step
2054 !-----------------------------------------------------------------------------
2055 subroutine RESPA_vel
2056 ! First and last RESPA step (incrementing velocities using long-range
2059 ! implicit real*8 (a-h,o-z)
2060 ! include 'DIMENSIONS'
2061 ! include 'COMMON.CONTROL'
2062 ! include 'COMMON.VAR'
2063 ! include 'COMMON.MD'
2064 ! include 'COMMON.CHAIN'
2065 ! include 'COMMON.DERIV'
2066 ! include 'COMMON.GEO'
2067 ! include 'COMMON.LOCAL'
2068 ! include 'COMMON.INTERACT'
2069 ! include 'COMMON.IOUNITS'
2070 ! include 'COMMON.NAMES'
2071 integer :: i,j,inres,mnum
2074 d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
2078 d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
2083 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2084 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2085 .and.(mnum.lt.4)) then
2088 d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
2093 end subroutine RESPA_vel
2094 !-----------------------------------------------------------------------------
2096 ! Applying velocity Verlet algorithm - step 1 to coordinates
2098 ! implicit real*8 (a-h,o-z)
2099 ! include 'DIMENSIONS'
2100 ! include 'COMMON.CONTROL'
2101 ! include 'COMMON.VAR'
2102 ! include 'COMMON.MD'
2103 ! include 'COMMON.CHAIN'
2104 ! include 'COMMON.DERIV'
2105 ! include 'COMMON.GEO'
2106 ! include 'COMMON.LOCAL'
2107 ! include 'COMMON.INTERACT'
2108 ! include 'COMMON.IOUNITS'
2109 ! include 'COMMON.NAMES'
2110 real(kind=8) :: adt,adt2
2111 integer :: i,j,inres,mnum
2114 write (iout,*) "VELVERLET1 START: DC"
2116 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),&
2117 (dc(j,i+nres),j=1,3)
2121 adt=d_a_old(j,0)*d_time
2123 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
2124 d_t_new(j,0)=d_t_old(j,0)+adt2
2125 d_t(j,0)=d_t_old(j,0)+adt
2129 adt=d_a_old(j,i)*d_time
2131 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
2132 d_t_new(j,i)=d_t_old(j,i)+adt2
2133 d_t(j,i)=d_t_old(j,i)+adt
2138 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2139 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2140 .and.(mnum.lt.4)) then
2143 adt=d_a_old(j,inres)*d_time
2145 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
2146 d_t_new(j,inres)=d_t_old(j,inres)+adt2
2147 d_t(j,inres)=d_t_old(j,inres)+adt
2152 write (iout,*) "VELVERLET1 END: DC"
2154 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),&
2155 (dc(j,i+nres),j=1,3)
2159 end subroutine verlet1
2160 !-----------------------------------------------------------------------------
2162 ! Step 2 of the velocity Verlet algorithm: update velocities
2164 ! implicit real*8 (a-h,o-z)
2165 ! include 'DIMENSIONS'
2166 ! include 'COMMON.CONTROL'
2167 ! include 'COMMON.VAR'
2168 ! include 'COMMON.MD'
2169 ! include 'COMMON.CHAIN'
2170 ! include 'COMMON.DERIV'
2171 ! include 'COMMON.GEO'
2172 ! include 'COMMON.LOCAL'
2173 ! include 'COMMON.INTERACT'
2174 ! include 'COMMON.IOUNITS'
2175 ! include 'COMMON.NAMES'
2176 integer :: i,j,inres,mnum
2179 d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
2183 d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
2188 ! iti=iabs(itype(i,mnum))
2189 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2190 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2191 .and.(mnum.lt.4)) then
2194 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
2199 end subroutine verlet2
2200 !-----------------------------------------------------------------------------
2201 subroutine sddir_precalc
2202 ! Applying velocity Verlet algorithm - step 1 to coordinates
2203 ! implicit real*8 (a-h,o-z)
2204 ! include 'DIMENSIONS'
2210 ! include 'COMMON.CONTROL'
2211 ! include 'COMMON.VAR'
2212 ! include 'COMMON.MD'
2214 ! include 'COMMON.LANGEVIN'
2216 ! include 'COMMON.LANGEVIN.lang0'
2218 ! include 'COMMON.CHAIN'
2219 ! include 'COMMON.DERIV'
2220 ! include 'COMMON.GEO'
2221 ! include 'COMMON.LOCAL'
2222 ! include 'COMMON.INTERACT'
2223 ! include 'COMMON.IOUNITS'
2224 ! include 'COMMON.NAMES'
2225 ! include 'COMMON.TIME1'
2226 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
2227 !el common /stochcalc/ stochforcvec
2228 real(kind=8) :: time00
2231 ! Compute friction and stochastic forces
2235 if (large) print *,"before friction_force"
2237 if (large) print *,"after friction_force"
2238 time_fric=time_fric+MPI_Wtime()-time00
2240 call stochastic_force(stochforcvec)
2241 time_stoch=time_stoch+MPI_Wtime()-time00
2244 ! Compute the acceleration due to friction forces (d_af_work) and stochastic
2245 ! forces (d_as_work)
2247 ! call ginv_mult(fric_work, d_af_work)
2248 ! call ginv_mult(stochforcvec, d_as_work)
2250 call fivediaginv_mult(dimen,fric_work, d_af_work)
2251 call fivediaginv_mult(dimen,stochforcvec, d_as_work)
2253 write(iout,*),"dimen",dimen
2255 write(iout,*) "fricwork",fric_work(i), d_af_work(i)
2256 write(iout,*) "stochforcevec", stochforcvec(i), d_as_work(i)
2260 call ginv_mult(fric_work, d_af_work)
2261 call ginv_mult(stochforcvec, d_as_work)
2265 end subroutine sddir_precalc
2266 !-----------------------------------------------------------------------------
2267 subroutine sddir_verlet1
2268 ! Applying velocity Verlet algorithm - step 1 to velocities
2271 ! implicit real*8 (a-h,o-z)
2272 ! include 'DIMENSIONS'
2273 ! include 'COMMON.CONTROL'
2274 ! include 'COMMON.VAR'
2275 ! include 'COMMON.MD'
2277 ! include 'COMMON.LANGEVIN'
2279 ! include 'COMMON.LANGEVIN.lang0'
2281 ! include 'COMMON.CHAIN'
2282 ! include 'COMMON.DERIV'
2283 ! include 'COMMON.GEO'
2284 ! include 'COMMON.LOCAL'
2285 ! include 'COMMON.INTERACT'
2286 ! include 'COMMON.IOUNITS'
2287 ! include 'COMMON.NAMES'
2288 ! Revised 3/31/05 AL: correlation between random contributions to
2289 ! position and velocity increments included.
2290 real(kind=8) :: sqrt13 = 0.57735026918962576451d0 ! 1/sqrt(3)
2291 real(kind=8) :: adt,adt2
2292 integer :: i,j,ind,inres,mnum
2294 ! Add the contribution from BOTH friction and stochastic force to the
2295 ! coordinates, but ONLY the contribution from the friction forces to velocities
2298 adt=(d_a_old(j,0)+d_af_work(j))*d_time
2299 adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
2300 ! write(iout,*) i,"adt",adt,"ads",adt2,d_a_old(j,0),d_af_work(j),d_time
2301 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
2302 d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
2303 d_t(j,0)=d_t_old(j,0)+adt
2308 adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
2309 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
2310 ! write(iout,*) i,"adt",adt,"ads",adt2,d_a_old(j,i),d_af_work(ind+j)
2311 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
2312 d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
2313 d_t(j,i)=d_t_old(j,i)+adt
2319 ! iti=iabs(itype(i,mnum))
2320 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2321 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2322 .and.(mnum.lt.4)) then
2325 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
2326 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
2327 ! write(iout,*) i,"adt",adt,"ads",adt2,d_a_old(j,inres),d_af_work(ind+j)
2328 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
2329 d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
2330 d_t(j,inres)=d_t_old(j,inres)+adt
2337 end subroutine sddir_verlet1
2338 !-----------------------------------------------------------------------------
2339 subroutine sddir_verlet2
2340 ! Calculating the adjusted velocities for accelerations
2343 ! implicit real*8 (a-h,o-z)
2344 ! include 'DIMENSIONS'
2345 ! include 'COMMON.CONTROL'
2346 ! include 'COMMON.VAR'
2347 ! include 'COMMON.MD'
2349 ! include 'COMMON.LANGEVIN'
2351 ! include 'COMMON.LANGEVIN.lang0'
2353 ! include 'COMMON.CHAIN'
2354 ! include 'COMMON.DERIV'
2355 ! include 'COMMON.GEO'
2356 ! include 'COMMON.LOCAL'
2357 ! include 'COMMON.INTERACT'
2358 ! include 'COMMON.IOUNITS'
2359 ! include 'COMMON.NAMES'
2360 real(kind=8),dimension(6*nres) :: stochforcvec,d_as_work1 !(MAXRES6) maxres6=6*maxres
2361 real(kind=8) :: cos60 = 0.5d0, sin60 = 0.86602540378443864676d0
2362 integer :: i,j,ind,inres,mnum
2363 ! Revised 3/31/05 AL: correlation between random contributions to
2364 ! position and velocity increments included.
2365 ! The correlation coefficients are calculated at low-friction limit.
2366 ! Also, friction forces are now not calculated with new velocities.
2368 ! call friction_force
2369 call stochastic_force(stochforcvec)
2371 ! Compute the acceleration due to friction forces (d_af_work) and stochastic
2372 ! forces (d_as_work)
2375 call fivediaginv_mult(6*nres,stochforcvec, d_as_work1)
2377 call ginv_mult(stochforcvec, d_as_work1)
2384 d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j)) &
2385 +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
2390 d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j)) &
2391 +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
2397 ! iti=iabs(itype(i,mnum))
2398 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2399 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2400 .and.(mnum.lt.4)) then
2403 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres) &
2404 +d_af_work(ind+j))+sin60*d_as_work(ind+j) &
2405 +cos60*d_as_work1(ind+j))*d_time
2411 end subroutine sddir_verlet2
2412 !-----------------------------------------------------------------------------
2413 subroutine max_accel
2415 ! Find the maximum difference in the accelerations of the the sites
2416 ! at the beginning and the end of the time step.
2419 ! implicit real*8 (a-h,o-z)
2420 ! include 'DIMENSIONS'
2421 ! include 'COMMON.CONTROL'
2422 ! include 'COMMON.VAR'
2423 ! include 'COMMON.MD'
2424 ! include 'COMMON.CHAIN'
2425 ! include 'COMMON.DERIV'
2426 ! include 'COMMON.GEO'
2427 ! include 'COMMON.LOCAL'
2428 ! include 'COMMON.INTERACT'
2429 ! include 'COMMON.IOUNITS'
2430 real(kind=8),dimension(3) :: aux,accel,accel_old
2431 real(kind=8) :: dacc
2435 ! aux(j)=d_a(j,0)-d_a_old(j,0)
2436 accel_old(j)=d_a_old(j,0)
2443 ! 7/3/08 changed to asymmetric difference
2445 ! accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
2446 accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
2447 accel(j)=accel(j)+0.5d0*d_a(j,i)
2448 ! if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
2449 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
2450 dacc=dabs(accel(j)-accel_old(j))
2451 ! write (iout,*) i,dacc
2452 if (dacc.gt.amax) amax=dacc
2460 accel_old(j)=d_a_old(j,0)
2465 accel_old(j)=accel_old(j)+d_a_old(j,1)
2466 accel(j)=accel(j)+d_a(j,1)
2471 ! iti=iabs(itype(i,mnum))
2472 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2473 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2474 .and.(mnum.lt.4)) then
2476 ! accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
2477 accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
2478 accel(j)=accel(j)+d_a(j,i+nres)
2482 ! if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
2483 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
2484 dacc=dabs(accel(j)-accel_old(j))
2485 ! write (iout,*) "side-chain",i,dacc
2486 if (dacc.gt.amax) amax=dacc
2490 accel_old(j)=accel_old(j)+d_a_old(j,i)
2491 accel(j)=accel(j)+d_a(j,i)
2492 ! aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
2496 end subroutine max_accel
2497 !-----------------------------------------------------------------------------
2498 subroutine predict_edrift(epdrift)
2500 ! Predict the drift of the potential energy
2503 use control_data, only: lmuca
2504 ! implicit real*8 (a-h,o-z)
2505 ! include 'DIMENSIONS'
2506 ! include 'COMMON.CONTROL'
2507 ! include 'COMMON.VAR'
2508 ! include 'COMMON.MD'
2509 ! include 'COMMON.CHAIN'
2510 ! include 'COMMON.DERIV'
2511 ! include 'COMMON.GEO'
2512 ! include 'COMMON.LOCAL'
2513 ! include 'COMMON.INTERACT'
2514 ! include 'COMMON.IOUNITS'
2515 ! include 'COMMON.MUCA'
2516 real(kind=8) :: epdrift,epdriftij
2518 ! Drift of the potential energy
2524 epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
2525 if (lmuca) epdriftij=epdriftij*factor
2526 ! write (iout,*) "back",i,j,epdriftij
2527 if (epdriftij.gt.epdrift) epdrift=epdriftij
2531 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1.and.&
2532 molnum(i).lt.4) then
2535 dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
2536 if (lmuca) epdriftij=epdriftij*factor
2537 ! write (iout,*) "side",i,j,epdriftij
2538 if (epdriftij.gt.epdrift) epdrift=epdriftij
2542 epdrift=0.5d0*epdrift*d_time*d_time
2543 ! write (iout,*) "epdrift",epdrift
2545 end subroutine predict_edrift
2546 !-----------------------------------------------------------------------------
2547 subroutine verlet_bath
2549 ! Coupling to the thermostat by using the Berendsen algorithm
2552 ! implicit real*8 (a-h,o-z)
2553 ! include 'DIMENSIONS'
2554 ! include 'COMMON.CONTROL'
2555 ! include 'COMMON.VAR'
2556 ! include 'COMMON.MD'
2557 ! include 'COMMON.CHAIN'
2558 ! include 'COMMON.DERIV'
2559 ! include 'COMMON.GEO'
2560 ! include 'COMMON.LOCAL'
2561 ! include 'COMMON.INTERACT'
2562 ! include 'COMMON.IOUNITS'
2563 ! include 'COMMON.NAMES'
2564 real(kind=8) :: T_half,fact
2565 integer :: i,j,inres,mnum
2567 T_half=2.0d0/(dimen3*Rb)*EK
2568 fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
2569 ! write(iout,*) "T_half", T_half
2570 ! write(iout,*) "EK", EK
2571 ! write(iout,*) "fact", fact
2573 d_t(j,0)=fact*d_t(j,0)
2577 d_t(j,i)=fact*d_t(j,i)
2582 ! iti=iabs(itype(i,mnum))
2583 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2584 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2585 .and.(mnum.lt.4)) then
2588 d_t(j,inres)=fact*d_t(j,inres)
2593 end subroutine verlet_bath
2594 !-----------------------------------------------------------------------------
2596 ! Set up the initial conditions of a MD simulation
2599 use control, only:tcpu
2600 !el use io_basic, only:ilen
2603 use minimm, only:minim_dc,minimize,sc_move
2604 use io_config, only:readrst
2605 use io, only:statout
2606 use random, only: iran_num
2607 ! implicit real*8 (a-h,o-z)
2608 ! include 'DIMENSIONS'
2611 character(len=16) :: form
2612 integer :: IERROR,ERRCODE
2614 integer :: iranmin,itrial,itmp,n_model_try,k, &
2616 integer, dimension(:),allocatable :: list_model_try
2617 integer, dimension(0:nodes-1) :: i_start_models
2618 ! include 'COMMON.SETUP'
2619 ! include 'COMMON.CONTROL'
2620 ! include 'COMMON.VAR'
2621 ! include 'COMMON.MD'
2623 ! include 'COMMON.LANGEVIN'
2625 ! include 'COMMON.LANGEVIN.lang0'
2627 ! include 'COMMON.CHAIN'
2628 ! include 'COMMON.DERIV'
2629 ! include 'COMMON.GEO'
2630 ! include 'COMMON.LOCAL'
2631 ! include 'COMMON.INTERACT'
2632 ! include 'COMMON.IOUNITS'
2633 ! include 'COMMON.NAMES'
2634 ! include 'COMMON.REMD'
2635 real(kind=8),dimension(0:n_ene) :: energia_long,energia_short,energia
2636 real(kind=8),dimension(3) :: vcm,incr,L
2637 real(kind=8) :: xv,sigv,lowb,highb
2638 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
2639 character(len=256) :: qstr
2642 character(len=50) :: tytul
2643 logical :: file_exist
2644 !el common /gucio/ cm
2645 integer :: i,j,ipos,iq,iw,nft_sc,iretcode,nfun,ierr,mnum,itime
2646 real(kind=8) :: etot,tt0
2650 ! write(iout,*) "d_time", d_time
2651 ! Compute the standard deviations of stochastic forces for Langevin dynamics
2652 ! if the friction coefficients do not depend on surface area
2653 if (lang.gt.0 .and. .not.surfarea) then
2656 stdforcp(i)=stdfp(mnum)*dsqrt(gamp(mnum))
2660 stdforcsc(i)=stdfsc(iabs(itype(i,mnum)),mnum) &
2661 *dsqrt(gamsc(iabs(itype(i,mnum)),mnum))
2665 ! Open the pdb file for snapshotshots
2668 if (ilen(tmpdir).gt.0) &
2669 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2670 liczba(:ilen(liczba))//".pdb")
2672 file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2676 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) &
2677 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2678 liczba(:ilen(liczba))//".x")
2679 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2682 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) &
2683 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2684 liczba(:ilen(liczba))//".cx")
2685 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2691 if (ilen(tmpdir).gt.0) &
2692 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
2693 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
2695 if (ilen(tmpdir).gt.0) &
2696 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
2697 cartname=prefix(:ilen(prefix))//"_MD.cx"
2701 write (qstr,'(256(1h ))')
2704 iq = qinfrag(i,iset)*10
2705 iw = wfrag(i,iset)/100
2707 if(me.eq.king.or..not.out1file) &
2708 write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
2709 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
2714 iq = qinpair(i,iset)*10
2715 iw = wpair(i,iset)/100
2717 if(me.eq.king.or..not.out1file) &
2718 write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
2719 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
2723 ! pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
2725 ! cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
2727 ! cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
2729 ! statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
2733 if (restart1file) then
2735 inquire(file=mremd_rst_name,exist=file_exist)
2736 write (*,*) me," Before broadcast: file_exist",file_exist
2738 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,&
2741 write (*,*) me," After broadcast: file_exist",file_exist
2742 ! inquire(file=mremd_rst_name,exist=file_exist)
2743 if(me.eq.king.or..not.out1file) &
2744 write(iout,*) "Initial state read by master and distributed"
2746 if (ilen(tmpdir).gt.0) &
2747 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_' &
2748 //liczba(:ilen(liczba))//'.rst')
2749 inquire(file=rest2name,exist=file_exist)
2752 if(.not.restart1file) then
2753 if(me.eq.king.or..not.out1file) &
2754 write(iout,*) "Initial state will be read from file ",&
2755 rest2name(:ilen(rest2name))
2758 call rescale_weights(t_bath)
2760 if(me.eq.king.or..not.out1file)then
2761 if (restart1file) then
2762 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),&
2765 write(iout,*) "File ",rest2name(:ilen(rest2name)),&
2768 write(iout,*) "Initial velocities randomly generated"
2775 ! Generate initial velocities
2776 if(me.eq.king.or..not.out1file) &
2777 write(iout,*) "Initial velocities randomly generated"
2782 ! rest2name = prefix(:ilen(prefix))//'.rst'
2783 if(me.eq.king.or..not.out1file)then
2784 write (iout,*) "Initial velocities"
2786 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2787 (d_t(j,i+nres),j=1,3)
2789 ! Zeroing the total angular momentum of the system
2790 write(iout,*) "Calling the zero-angular momentum subroutine"
2793 ! Getting the potential energy and forces and velocities and accelerations
2795 ! write (iout,*) "velocity of the center of the mass:"
2796 ! write (iout,*) (vcm(j),j=1,3)
2798 d_t(j,0)=d_t(j,0)-vcm(j)
2800 ! Removing the velocity of the center of mass
2802 if(me.eq.king.or..not.out1file)then
2803 write (iout,*) "vcm right after adjustment:"
2804 write (iout,*) (vcm(j),j=1,3)
2811 if ((.not.rest).or.(forceminim)) then
2812 if (forceminim) call chainbuild_cart
2814 if(iranconf.ne.0 .or.indpdb.gt.0.and..not.unres_pdb .or.preminim) then
2816 print *, 'Calling OVERLAP_SC'
2817 call overlap_sc(fail)
2818 print *,'after OVERLAP'
2821 print *,'call SC_MOVE'
2822 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
2823 print *,'SC_move',nft_sc,etot
2824 if(me.eq.king.or..not.out1file) &
2825 write(iout,*) 'SC_move',nft_sc,etot
2829 print *, 'Calling MINIM_DC'
2830 call minim_dc(etot,iretcode,nfun)
2832 call geom_to_var(nvar,varia)
2833 print *,'Calling MINIMIZE.'
2834 call minimize(etot,varia,iretcode,nfun)
2835 call var_to_geom(nvar,varia)
2837 write(iout,*) "just before minimin"
2839 if(me.eq.king.or..not.out1file) &
2840 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2842 write(iout,*) "just after minimin"
2844 if(iranconf.ne.0) then
2845 !c 8/22/17 AL Loop to produce a low-energy random conformation
2848 if(me.eq.king.or..not.out1file) &
2849 write (iout,*) 'Calling OVERLAP_SC'
2850 call overlap_sc(fail)
2851 endif !endif overlap
2854 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
2855 print *,'SC_move',nft_sc,etot
2856 if(me.eq.king.or..not.out1file) &
2857 write(iout,*) 'SC_move',nft_sc,etot
2861 print *, 'Calling MINIM_DC'
2862 call minim_dc(etot,iretcode,nfun)
2863 call int_from_cart1(.false.)
2865 call geom_to_var(nvar,varia)
2866 print *,'Calling MINIMIZE.'
2867 call minimize(etot,varia,iretcode,nfun)
2868 call var_to_geom(nvar,varia)
2870 if(me.eq.king.or..not.out1file) &
2871 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2872 write(iout,*) "just after minimin"
2874 if (isnan(etot) .or. etot.gt.4.0d6) then
2875 write (iout,*) "Energy too large",etot, &
2876 " trying another random conformation"
2879 call gen_rand_conf(itmp,*30)
2881 30 write (iout,*) 'Failed to generate random conformation', &
2883 write (*,*) 'Processor:',me, &
2884 ' Failed to generate random conformation',&
2893 write (iout,'(a,i3,a)') 'Processor:',me, &
2894 ' error in generating random conformation.'
2895 write (*,'(a,i3,a)') 'Processor:',me, &
2896 ' error in generating random conformation.'
2899 ! call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
2900 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
2910 write (iout,'(a,i3,a)') 'Processor:',me, &
2911 ' failed to generate a low-energy random conformation.'
2912 write (*,'(a,i3,a,f10.3)') 'Processor:',me, &
2913 ' failed to generate a low-energy random conformation.',etot
2917 ! call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
2918 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
2923 else if (preminim) then
2924 if (start_from_model) then
2928 do while (fail .and. n_model_try.lt.nmodel_start)
2929 write (iout,*) "n_model_try",n_model_try
2931 i_model=iran_num(1,nmodel_start)
2933 if (i_model.eq.list_model_try(k)) exit
2935 if (k.gt.n_model_try) exit
2937 n_model_try=n_model_try+1
2938 list_model_try(n_model_try)=i_model
2939 if (me.eq.king .or. .not. out1file) &
2940 write (iout,*) 'Trying to start from model ',&
2941 pdbfiles_chomo(i_model)(:ilen(pdbfiles_chomo(i_model)))
2944 c(j,i)=chomo(j,i,i_model)
2947 call int_from_cart(.true.,.false.)
2948 call sc_loc_geom(.false.)
2952 dc(j,i)=c(j,i+1)-c(j,i)
2953 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
2958 dc(j,i+nres)=c(j,i+nres)-c(j,i)
2959 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
2962 if (me.eq.king.or..not.out1file) then
2963 write (iout,*) "Energies before removing overlaps"
2964 call etotal(energia(0))
2965 call enerprint(energia(0))
2967 ! Remove SC overlaps if requested
2969 write (iout,*) 'Calling OVERLAP_SC'
2970 call overlap_sc(fail)
2973 "Failed to remove overlap from model",i_model
2977 if (me.eq.king.or..not.out1file) then
2978 write (iout,*) "Energies after removing overlaps"
2979 call etotal(energia(0))
2980 call enerprint(energia(0))
2983 ! Search for better SC rotamers if requested
2985 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
2986 print *,'SC_move',nft_sc,etot
2987 if (me.eq.king.or..not.out1file)&
2988 write(iout,*) 'SC_move',nft_sc,etot
2990 call etotal(energia(0))
2993 call MPI_Gather(i_model,1,MPI_INTEGER,i_start_models(0),&
2994 1,MPI_INTEGER,king,CG_COMM,IERROR)
2995 if (n_model_try.gt.nmodel_start .and.&
2996 (me.eq.king .or. out1file)) then
2998 "All models have irreparable overlaps. Trying randoms starts."
3000 i_model=nmodel_start+1
3004 ! Remove SC overlaps if requested
3006 write (iout,*) 'Calling OVERLAP_SC'
3007 call overlap_sc(fail)
3010 "Failed to remove overlap"
3013 if (me.eq.king.or..not.out1file) then
3014 write (iout,*) "Energies after removing overlaps"
3015 call etotal(energia(0))
3016 call enerprint(energia(0))
3019 ! 8/22/17 AL Minimize initial structure
3021 if (me.eq.king.or..not.out1file) write(iout,*)&
3022 'Minimizing initial PDB structure: Calling MINIM_DC'
3023 call minim_dc(etot,iretcode,nfun)
3025 call geom_to_var(nvar,varia)
3026 if(me.eq.king.or..not.out1file) write (iout,*)&
3027 'Minimizing initial PDB structure: Calling MINIMIZE.'
3028 call minimize(etot,varia,iretcode,nfun)
3029 call var_to_geom(nvar,varia)
3031 if (me.eq.king.or..not.out1file)&
3032 write(iout,*) 'LBFGS return code is ',status,' eval ',nfun
3033 if(me.eq.king.or..not.out1file)&
3034 write(iout,*) 'LBFGS return code is ',status,' eval ',nfun
3036 if (me.eq.king.or..not.out1file)&
3037 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
3038 if(me.eq.king.or..not.out1file)&
3039 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
3043 if (nmodel_start.gt.0 .and. me.eq.king) then
3044 write (iout,'(a)') "Task Starting model"
3046 if (i_start_models(i).gt.nmodel_start) then
3047 write (iout,'(i4,2x,a)') i,"RANDOM STRUCTURE"
3049 write(iout,'(i4,2x,a)')i,pdbfiles_chomo(i_start_models(i)) &
3050 (:ilen(pdbfiles_chomo(i_start_models(i))))
3055 call chainbuild_cart
3057 write(iout,*) "just after kinetic"
3062 kinetic_T=2.0d0/(dimen3*Rb)*EK
3063 if(me.eq.king.or..not.out1file)then
3064 write(iout,*) "just after verlet_bath"
3074 write(iout,*) "before ETOTAL"
3075 call etotal(potEcomp)
3076 if (large) call enerprint(potEcomp)
3079 t_etotal=t_etotal+MPI_Wtime()-tt0
3081 t_etotal=t_etotal+tcpu()-tt0
3086 write(iout,*) "before lagrangian"
3088 write(iout,*) "before max_accel"
3090 if (amax*d_time .gt. dvmax) then
3091 d_time=d_time*dvmax/amax
3092 if(me.eq.king.or..not.out1file) write (iout,*) &
3093 "Time step reduced to",d_time,&
3094 " because of too large initial acceleration."
3096 if(me.eq.king.or..not.out1file)then
3097 write(iout,*) "Potential energy and its components"
3098 call enerprint(potEcomp)
3099 ! write(iout,*) (potEcomp(i),i=0,n_ene)
3101 potE=potEcomp(0)-potEcomp(51)
3105 if (ntwe.ne.0) call statout(itime)
3106 if(me.eq.king.or..not.out1file) &
3107 write (iout,'(/a/3(a25,1pe14.5/))') "Initial:", &
3108 " Kinetic energy",EK," Potential energy",potE, &
3109 " Total energy",totE," Maximum acceleration ", &
3112 write (iout,*) "Initial coordinates"
3114 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),&
3117 write (iout,*) "Initial dC"
3119 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),&
3120 (dc(j,i+nres),j=1,3)
3122 write (iout,*) "Initial velocities"
3123 write (iout,"(13x,' backbone ',23x,' side chain')")
3125 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
3126 (d_t(j,i+nres),j=1,3)
3128 write (iout,*) "Initial accelerations"
3130 ! write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
3131 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),&
3132 (d_a(j,i+nres),j=1,3)
3138 d_t_old(j,i)=d_t(j,i)
3139 d_a_old(j,i)=d_a(j,i)
3141 ! write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
3150 call etotal_short(energia_short)
3151 if (large) call enerprint(potEcomp)
3154 t_eshort=t_eshort+MPI_Wtime()-tt0
3156 t_eshort=t_eshort+tcpu()-tt0
3161 if(.not.out1file .and. large) then
3162 write (iout,*) "energia_long",energia_long(0),&
3163 " energia_short",energia_short(0),&
3164 " total",energia_long(0)+energia_short(0)
3165 write (iout,*) "Initial fast-force accelerations"
3167 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
3168 (d_a(j,i+nres),j=1,3)
3171 ! 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
3174 d_a_short(j,i)=d_a(j,i)
3183 call etotal_long(energia_long)
3184 if (large) call enerprint(potEcomp)
3187 t_elong=t_elong+MPI_Wtime()-tt0
3189 t_elong=t_elong+tcpu()-tt0
3194 if(.not.out1file .and. large) then
3195 write (iout,*) "energia_long",energia_long(0)
3196 write (iout,*) "Initial slow-force accelerations"
3198 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
3199 (d_a(j,i+nres),j=1,3)
3203 t_enegrad=t_enegrad+MPI_Wtime()-tt0
3205 t_enegrad=t_enegrad+tcpu()-tt0
3209 end subroutine init_MD
3210 !-----------------------------------------------------------------------------
3211 subroutine random_vel
3213 ! implicit real*8 (a-h,o-z)
3215 use random, only:anorm_distr
3217 ! include 'DIMENSIONS'
3218 ! include 'COMMON.CONTROL'
3219 ! include 'COMMON.VAR'
3220 ! include 'COMMON.MD'
3222 ! include 'COMMON.LANGEVIN'
3224 ! include 'COMMON.LANGEVIN.lang0'
3226 ! include 'COMMON.CHAIN'
3227 ! include 'COMMON.DERIV'
3228 ! include 'COMMON.GEO'
3229 ! include 'COMMON.LOCAL'
3230 ! include 'COMMON.INTERACT'
3231 ! include 'COMMON.IOUNITS'
3232 ! include 'COMMON.NAMES'
3233 ! include 'COMMON.TIME1'
3234 real(kind=8) :: xv,sigv,lowb,highb ,Ek1
3236 integer ichain,n,innt,inct,ibeg,ierr
3237 double precision work(48*nres)
3238 integer iwork(6*nres)
3239 ! double precision Ghalf(mmaxres2_chain),Geigen(maxres2_chain),&
3240 ! Gvec(maxres2_chain,maxres2_chain)
3241 ! common /przechowalnia/Ghalf,Geigen,Gvec
3243 double precision inertia(maxres2_chain,maxres2_chain)
3248 real(kind=8) ,allocatable, dimension(:) :: DDU1,DDU2,DL2,DL1,xsolv,DML,rs
3249 real(kind=8) :: sumx,Ek2,Ek3,aux
3251 real(kind=8) ,allocatable, dimension(:) :: rsold
3252 real (kind=8),allocatable,dimension(:,:) :: matold,inertia
3255 integer :: i,j,ii,k,ind,mark,imark,mnum
3256 ! Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
3257 ! First generate velocities in the eigenspace of the G matrix
3258 ! write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
3262 write (iout,*) "Random_vel, fivediag"
3263 allocate(inertia(2*nres,2*nres))
3270 write(iout,*), nchain
3275 n=dimen_chain(ichain)
3276 innt=iposd_chain(ichain)
3279 write (iout,*) "Chain",ichain," n",n," start",innt
3281 if (i.lt.inct-1) then
3282 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i),DU1orig(i),&
3284 else if (i.eq.inct-1) then
3285 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i),DU1orig(i)
3287 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i)
3292 ghalf(ind+1)=dmorig(innt)
3293 ghalf(ind+2)=du1orig(innt)
3294 ghalf(ind+3)=dmorig(innt+1)
3298 write (iout,*) "i",i," ind",ind," indu2",innt+i-2,&
3299 " indu1",innt+i-1," indm",innt+i
3300 ghalf(ind+1)=du2orig(innt-1+i-2)
3301 ghalf(ind+2)=du1orig(innt-1+i-1)
3302 ghalf(ind+3)=dmorig(innt-1+i)
3303 !c write (iout,'(3(a,i2,1x))') "DU2",innt-1+i-2,
3304 !c & "DU1",innt-1+i-1,"DM ",innt-1+i
3312 inertia(i,j)=ghalf(ind)
3313 inertia(j,i)=ghalf(ind)
3318 write (iout,*) "Chain ",ichain," ind",ind," dim",n*(n+1)/2
3319 write (iout,*) "Five-diagonal inertia matrix, lower triangle"
3320 ! call matoutr(n,ghalf)
3322 call gldiag(nres*2,n,n,Ghalf,work,Geigen,Gvec,ierr,iwork)
3324 write (iout,'(//a,i3)')&
3325 "Eigenvectors and eigenvalues of the G matrix chain",ichain
3326 call eigout(n,n,nres*2,nres*2,Gvec,Geigen)
3329 !c check diagonalization
3335 aux=aux+gvec(k,i)*gvec(l,j)*inertia(k,l)
3339 write (iout,*) i,j,aux,geigen(i)
3341 write (iout,*) i,j,aux
3351 sigv=dsqrt((Rb*t_bath)/geigen(i))
3354 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
3355 EK=EK+0.5d0*geigen(i)*d_t_work_new(ii)**2
3356 !c write (iout,*) "i",i," ii",ii," geigen",geigen(i),
3357 !c & " d_t_work_new",d_t_work_new(ii)
3365 d_t_work(ind)=d_t_work(ind)&
3366 +Gvec(i,j)*d_t_work_new((j-1)*3+k)
3368 !c write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
3377 aux=aux+inertia(i,j)*d_t_work(3*(i-1)+k)*d_t_work(3*(j-1)+k)
3383 !c Transfer to the d_t vector
3384 innt=chain_border(1,ichain)
3385 inct=chain_border(2,ichain)
3387 !c write (iout,*) "ichain",ichain," innt",innt," inct",inct
3391 d_t(j,i)=d_t_work(ind)
3394 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum).and.mnum.le.2) then
3397 d_t(j,i+nres)=d_t_work(ind)
3404 write (iout,*) "Random velocities in the Calpha,SC space"
3407 write (iout,'(a3,1h(,i5,1h),3f10.5,3x,3f10.5)')&
3408 restyp(itype(i,mnum),mnum),i,(d_t(j,i),j=1,3),(d_t(j,i+nres),j=1,3)
3411 call kinetic_CASC(Ek1)
3413 ! Transform the velocities to virtual-bond space
3422 if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum).or.mnum.ge.3) then
3424 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
3428 d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
3429 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
3440 !c d_a(:,0)=d_a(:,1)
3442 !c write (iout,*) "Shifting accelerations"
3444 !c write (iout,*) "ichain",chain_border1(1,ichain)-1,
3445 !c & chain_border1(1,ichain)
3446 d_t(:,chain_border1(1,ichain)-1)=d_t(:,chain_border1(1,ichain))
3447 d_t(:,chain_border1(1,ichain))=0.0d0
3449 !c write (iout,*) "Adding accelerations"
3451 !c write (iout,*) "chain",ichain,chain_border1(1,ichain)-1,
3452 !c & chain_border(2,ichain-1)
3453 d_t(:,chain_border1(1,ichain)-1)=&
3454 d_t(:,chain_border1(1,ichain)-1)+d_t(:,chain_border(2,ichain-1))
3455 d_t(:,chain_border(2,ichain-1))=0.0d0
3458 write (iout,*) "chain",ichain,chain_border1(1,ichain)-1,&
3459 chain_border(2,ichain-1)
3460 d_t(:,chain_border1(1,ichain)-1)=&
3461 d_t(:,chain_border1(1,ichain)-1)+d_t(:,chain_border(2,ichain-1))
3462 d_t(:,chain_border(2,ichain-1))=0.0d0
3467 !c d_t(j,0)=d_t(j,nnt)
3470 innt=chain_border(1,ichain)
3471 inct=chain_border(2,ichain)
3472 !c write (iout,*) "ichain",ichain," innt",innt," inct",inct
3473 !c write (iout,*) "ibeg",ibeg
3475 d_t(j,ibeg)=d_t(j,innt)
3480 if (iabs(itype(i,1).eq.10).or.mnum.ge.3) then
3481 !c write (iout,*) "i",i,(d_t(j,i),j=1,3),(d_t(j,i+1),j=1,3)
3483 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
3487 d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
3488 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
3497 "Random velocities in the virtual-bond-vector space"
3498 write (iout,'(3hORG,1h(,i5,1h),3f10.5)') 0,(d_t(j,0),j=1,3)
3500 write (iout,'(a3,1h(,i5,1h),3f10.5,3x,3f10.5)')&
3501 restyp(itype(i,mnum),mnum),i,(d_t(j,i),j=1,3),(d_t(j,i+nres),j=1,3)
3504 write (iout,*) "Kinetic energy from inertia matrix eigenvalues",&
3507 "Kinetic temperatures from inertia matrix eigenvalues",&
3510 write (iout,*) "Kinetic energy from inertia matrix",Ek3
3511 write (iout,*) "Kinetic temperatures from inertia",&
3514 write (iout,*) "Kinetic energy from velocities in CA-SC space",&
3517 "Kinetic temperatures from velovities in CA-SC space",&
3521 "Kinetic energy from virtual-bond-vector velocities",Ek1
3523 "Kinetic temperature from virtual-bond-vector velocities ",&
3532 sigv=dsqrt((Rb*t_bath)/geigen(i))
3535 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
3537 write (iout,*) "i",i," ii",ii," geigen",geigen(i),&
3538 " d_t_work_new",d_t_work_new(ii)
3549 Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
3552 write (iout,*) "Ek from eigenvectors",Ek1
3553 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
3562 d_t_work(ind)=d_t_work(ind) &
3563 +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
3565 ! write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
3569 ! Transfer to the d_t vector
3571 d_t(j,0)=d_t_work(j)
3577 d_t(j,i)=d_t_work(ind)
3582 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3583 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3584 .and.(mnum.lt.4)) then
3587 d_t(j,i+nres)=d_t_work(ind)
3593 ! write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",&
3594 ! 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
3596 ! write(iout,*) "end init MD"
3599 end subroutine random_vel
3600 !-----------------------------------------------------------------------------
3602 subroutine sd_verlet_p_setup
3603 ! Sets up the parameters of stochastic Verlet algorithm
3604 ! implicit real*8 (a-h,o-z)
3605 ! include 'DIMENSIONS'
3606 use control, only: tcpu
3611 ! include 'COMMON.CONTROL'
3612 ! include 'COMMON.VAR'
3613 ! include 'COMMON.MD'
3615 ! include 'COMMON.LANGEVIN'
3617 ! include 'COMMON.LANGEVIN.lang0'
3619 ! include 'COMMON.CHAIN'
3620 ! include 'COMMON.DERIV'
3621 ! include 'COMMON.GEO'
3622 ! include 'COMMON.LOCAL'
3623 ! include 'COMMON.INTERACT'
3624 ! include 'COMMON.IOUNITS'
3625 ! include 'COMMON.NAMES'
3626 ! include 'COMMON.TIME1'
3627 real(kind=8),dimension(6*nres) :: emgdt !(MAXRES6) maxres6=6*maxres
3628 real(kind=8) :: pterm,vterm,rho,rhoc,vsig
3629 real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
3630 prand_vec,vrand_vec1,vrand_vec2 !(MAXRES6) maxres6=6*maxres
3631 logical :: lprn = .false.
3632 real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
3633 real(kind=8) :: ktm,gdt,egdt,gdt2,gdt3,gdt4,gdt5,gdt6,gdt7,gdt8,&
3635 integer :: i,maxres2
3642 ! AL 8/17/04 Code adapted from tinker
3644 ! Get the frictional and random terms for stochastic dynamics in the
3645 ! eigenspace of mass-scaled UNRES friction matrix
3649 gdt = fricgam(i) * d_time
3651 ! Stochastic dynamics reduces to simple MD for zero friction
3653 if (gdt .le. zero) then
3654 pfric_vec(i) = 1.0d0
3655 vfric_vec(i) = d_time
3656 afric_vec(i) = 0.5d0 * d_time * d_time
3657 prand_vec(i) = 0.0d0
3658 vrand_vec1(i) = 0.0d0
3659 vrand_vec2(i) = 0.0d0
3661 ! Analytical expressions when friction coefficient is large
3664 if (gdt .ge. gdt_radius) then
3667 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
3668 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
3669 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
3670 vterm = 1.0d0 - egdt**2
3671 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
3673 ! Use series expansions when friction coefficient is small
3684 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0 &
3685 - gdt5/120.0d0 + gdt6/720.0d0 &
3686 - gdt7/5040.0d0 + gdt8/40320.0d0 &
3687 - gdt9/362880.0d0) / fricgam(i)**2
3688 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
3689 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
3690 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0 &
3691 + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0 &
3692 + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0 &
3693 + 127.0d0*gdt9/90720.0d0
3694 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0 &
3695 - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0 &
3696 - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0 &
3697 - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
3698 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0 &
3699 - 17.0d0*gdt2/1280.0d0 &
3700 + 17.0d0*gdt3/6144.0d0 &
3701 + 40967.0d0*gdt4/34406400.0d0 &
3702 - 57203.0d0*gdt5/275251200.0d0 &
3703 - 1429487.0d0*gdt6/13212057600.0d0)
3706 ! Compute the scaling factors of random terms for the nonzero friction case
3708 ktm = 0.5d0*d_time/fricgam(i)
3709 psig = dsqrt(ktm*pterm) / fricgam(i)
3710 vsig = dsqrt(ktm*vterm)
3711 rhoc = dsqrt(1.0d0 - rho*rho)
3713 vrand_vec1(i) = vsig * rho
3714 vrand_vec2(i) = vsig * rhoc
3719 "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",&
3722 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),&
3723 afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
3727 ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
3730 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
3731 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
3732 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
3733 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
3734 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
3735 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
3738 t_sdsetup=t_sdsetup+MPI_Wtime()
3740 t_sdsetup=t_sdsetup+tcpu()-tt0
3743 end subroutine sd_verlet_p_setup
3744 !-----------------------------------------------------------------------------
3745 subroutine eigtransf1(n,ndim,ab,d,c)
3749 real(kind=8) :: ab(ndim,ndim,n),c(ndim,n),d(ndim)
3755 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
3760 end subroutine eigtransf1
3761 !-----------------------------------------------------------------------------
3762 subroutine eigtransf(n,ndim,a,b,d,c)
3766 real(kind=8) :: a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
3772 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
3777 end subroutine eigtransf
3778 !-----------------------------------------------------------------------------
3779 subroutine sd_verlet1
3781 ! Applying stochastic velocity Verlet algorithm - step 1 to velocities
3783 ! implicit real*8 (a-h,o-z)
3784 ! include 'DIMENSIONS'
3785 ! include 'COMMON.CONTROL'
3786 ! include 'COMMON.VAR'
3787 ! include 'COMMON.MD'
3789 ! include 'COMMON.LANGEVIN'
3791 ! include 'COMMON.LANGEVIN.lang0'
3793 ! include 'COMMON.CHAIN'
3794 ! include 'COMMON.DERIV'
3795 ! include 'COMMON.GEO'
3796 ! include 'COMMON.LOCAL'
3797 ! include 'COMMON.INTERACT'
3798 ! include 'COMMON.IOUNITS'
3799 ! include 'COMMON.NAMES'
3800 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
3801 !el common /stochcalc/ stochforcvec
3802 logical :: lprn = .false.
3803 real(kind=8) :: ddt1,ddt2
3804 integer :: i,j,ind,inres
3806 ! write (iout,*) "dc_old"
3808 ! write (iout,'(i5,3f10.5,5x,3f10.5)')
3809 ! & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3812 dc_work(j)=dc_old(j,0)
3813 d_t_work(j)=d_t_old(j,0)
3814 d_a_work(j)=d_a_old(j,0)
3819 dc_work(ind+j)=dc_old(j,i)
3820 d_t_work(ind+j)=d_t_old(j,i)
3821 d_a_work(ind+j)=d_a_old(j,i)
3827 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3828 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3829 .and.(mnum.lt.4)) then
3831 dc_work(ind+j)=dc_old(j,i+nres)
3832 d_t_work(ind+j)=d_t_old(j,i+nres)
3833 d_a_work(ind+j)=d_a_old(j,i+nres)
3841 "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",&
3845 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
3846 vfric_mat(i,j),afric_mat(i,j),&
3847 prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3855 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) &
3856 +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3857 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3858 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3860 d_t_work_new(i)=ddt1+0.5d0*ddt2
3861 d_t_work(i)=ddt1+ddt2
3866 d_t(j,0)=d_t_work(j)
3871 dc(j,i)=dc_work(ind+j)
3872 d_t(j,i)=d_t_work(ind+j)
3878 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3879 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3880 .and.(mnum.lt.4)) then
3883 dc(j,inres)=dc_work(ind+j)
3884 d_t(j,inres)=d_t_work(ind+j)
3890 end subroutine sd_verlet1
3891 !-----------------------------------------------------------------------------
3892 subroutine sd_verlet2
3894 ! Calculating the adjusted velocities for accelerations
3896 ! implicit real*8 (a-h,o-z)
3897 ! include 'DIMENSIONS'
3898 ! include 'COMMON.CONTROL'
3899 ! include 'COMMON.VAR'
3900 ! include 'COMMON.MD'
3902 ! include 'COMMON.LANGEVIN'
3904 ! include 'COMMON.LANGEVIN.lang0'
3906 ! include 'COMMON.CHAIN'
3907 ! include 'COMMON.DERIV'
3908 ! include 'COMMON.GEO'
3909 ! include 'COMMON.LOCAL'
3910 ! include 'COMMON.INTERACT'
3911 ! include 'COMMON.IOUNITS'
3912 ! include 'COMMON.NAMES'
3913 !el real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV !(MAXRES6) maxres6=6*maxres
3914 real(kind=8),dimension(6*nres) :: stochforcvecV !(MAXRES6) maxres6=6*maxres
3915 !el common /stochcalc/ stochforcvec
3917 real(kind=8) :: ddt1,ddt2
3918 integer :: i,j,ind,inres
3919 ! Compute the stochastic forces which contribute to velocity change
3921 call stochastic_force(stochforcvecV)
3928 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3929 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+ &
3930 vrand_mat2(i,j)*stochforcvecV(j)
3932 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
3936 d_t(j,0)=d_t_work(j)
3941 d_t(j,i)=d_t_work(ind+j)
3947 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3948 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3949 .and.(mnum.lt.4)) then
3952 d_t(j,inres)=d_t_work(ind+j)
3958 end subroutine sd_verlet2
3959 !-----------------------------------------------------------------------------
3960 subroutine sd_verlet_ciccotti_setup
3962 ! Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
3964 ! implicit real*8 (a-h,o-z)
3965 ! include 'DIMENSIONS'
3966 use control, only: tcpu
3971 ! include 'COMMON.CONTROL'
3972 ! include 'COMMON.VAR'
3973 ! include 'COMMON.MD'
3975 ! include 'COMMON.LANGEVIN'
3977 ! include 'COMMON.LANGEVIN.lang0'
3979 ! include 'COMMON.CHAIN'
3980 ! include 'COMMON.DERIV'
3981 ! include 'COMMON.GEO'
3982 ! include 'COMMON.LOCAL'
3983 ! include 'COMMON.INTERACT'
3984 ! include 'COMMON.IOUNITS'
3985 ! include 'COMMON.NAMES'
3986 ! include 'COMMON.TIME1'
3987 real(kind=8),dimension(6*nres) :: emgdt !(MAXRES6) maxres6=6*maxres
3988 real(kind=8) :: pterm,vterm,rho,rhoc,vsig
3989 real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
3990 prand_vec,vrand_vec1,vrand_vec2 !(MAXRES6) maxres6=6*maxres
3991 logical :: lprn = .false.
3992 real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
3993 real(kind=8) :: ktm,gdt,egdt,tt0
3994 integer :: i,maxres2
4001 ! AL 8/17/04 Code adapted from tinker
4003 ! Get the frictional and random terms for stochastic dynamics in the
4004 ! eigenspace of mass-scaled UNRES friction matrix
4008 write (iout,*) "i",i," fricgam",fricgam(i)
4009 gdt = fricgam(i) * d_time
4011 ! Stochastic dynamics reduces to simple MD for zero friction
4013 if (gdt .le. zero) then
4014 pfric_vec(i) = 1.0d0
4015 vfric_vec(i) = d_time
4016 afric_vec(i) = 0.5d0*d_time*d_time
4017 prand_vec(i) = afric_vec(i)
4018 vrand_vec2(i) = vfric_vec(i)
4020 ! Analytical expressions when friction coefficient is large
4025 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
4026 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
4027 prand_vec(i) = afric_vec(i)
4028 vrand_vec2(i) = vfric_vec(i)
4030 ! Compute the scaling factors of random terms for the nonzero friction case
4032 ! ktm = 0.5d0*d_time/fricgam(i)
4033 ! psig = dsqrt(ktm*pterm) / fricgam(i)
4034 ! vsig = dsqrt(ktm*vterm)
4035 ! prand_vec(i) = psig*afric_vec(i)
4036 ! vrand_vec2(i) = vsig*vfric_vec(i)
4041 "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",&
4044 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),&
4045 afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
4049 ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
4051 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
4052 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
4053 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
4054 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
4055 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
4057 t_sdsetup=t_sdsetup+MPI_Wtime()
4059 t_sdsetup=t_sdsetup+tcpu()-tt0
4062 end subroutine sd_verlet_ciccotti_setup
4063 !-----------------------------------------------------------------------------
4064 subroutine sd_verlet1_ciccotti
4066 ! Applying stochastic velocity Verlet algorithm - step 1 to velocities
4067 ! implicit real*8 (a-h,o-z)
4069 ! include 'DIMENSIONS'
4073 ! include 'COMMON.CONTROL'
4074 ! include 'COMMON.VAR'
4075 ! include 'COMMON.MD'
4077 ! include 'COMMON.LANGEVIN'
4079 ! include 'COMMON.LANGEVIN.lang0'
4081 ! include 'COMMON.CHAIN'
4082 ! include 'COMMON.DERIV'
4083 ! include 'COMMON.GEO'
4084 ! include 'COMMON.LOCAL'
4085 ! include 'COMMON.INTERACT'
4086 ! include 'COMMON.IOUNITS'
4087 ! include 'COMMON.NAMES'
4088 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
4089 !el common /stochcalc/ stochforcvec
4090 logical :: lprn = .false.
4091 real(kind=8) :: ddt1,ddt2
4092 integer :: i,j,ind,inres
4093 ! write (iout,*) "dc_old"
4095 ! write (iout,'(i5,3f10.5,5x,3f10.5)')
4096 ! & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
4099 dc_work(j)=dc_old(j,0)
4100 d_t_work(j)=d_t_old(j,0)
4101 d_a_work(j)=d_a_old(j,0)
4106 dc_work(ind+j)=dc_old(j,i)
4107 d_t_work(ind+j)=d_t_old(j,i)
4108 d_a_work(ind+j)=d_a_old(j,i)
4113 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
4115 dc_work(ind+j)=dc_old(j,i+nres)
4116 d_t_work(ind+j)=d_t_old(j,i+nres)
4117 d_a_work(ind+j)=d_a_old(j,i+nres)
4126 "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",&
4130 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
4131 vfric_mat(i,j),afric_mat(i,j),&
4132 prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
4140 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) &
4141 +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
4142 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
4143 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
4145 d_t_work_new(i)=ddt1+0.5d0*ddt2
4146 d_t_work(i)=ddt1+ddt2
4151 d_t(j,0)=d_t_work(j)
4156 dc(j,i)=dc_work(ind+j)
4157 d_t(j,i)=d_t_work(ind+j)
4163 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
4164 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4165 .and.(mnum.lt.4)) then
4168 dc(j,inres)=dc_work(ind+j)
4169 d_t(j,inres)=d_t_work(ind+j)
4175 end subroutine sd_verlet1_ciccotti
4176 !-----------------------------------------------------------------------------
4177 subroutine sd_verlet2_ciccotti
4179 ! Calculating the adjusted velocities for accelerations
4181 ! implicit real*8 (a-h,o-z)
4182 ! include 'DIMENSIONS'
4183 ! include 'COMMON.CONTROL'
4184 ! include 'COMMON.VAR'
4185 ! include 'COMMON.MD'
4187 ! include 'COMMON.LANGEVIN'
4189 ! include 'COMMON.LANGEVIN.lang0'
4191 ! include 'COMMON.CHAIN'
4192 ! include 'COMMON.DERIV'
4193 ! include 'COMMON.GEO'
4194 ! include 'COMMON.LOCAL'
4195 ! include 'COMMON.INTERACT'
4196 ! include 'COMMON.IOUNITS'
4197 ! include 'COMMON.NAMES'
4198 !el real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV !(MAXRES6) maxres6=6*maxres
4199 real(kind=8),dimension(6*nres) :: stochforcvecV !(MAXRES6) maxres6=6*maxres
4200 !el common /stochcalc/ stochforcvec
4201 real(kind=8) :: ddt1,ddt2
4202 integer :: i,j,ind,inres
4204 ! Compute the stochastic forces which contribute to velocity change
4206 call stochastic_force(stochforcvecV)
4213 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
4214 ! ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
4215 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
4217 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
4221 d_t(j,0)=d_t_work(j)
4226 d_t(j,i)=d_t_work(ind+j)
4232 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4234 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
4237 d_t(j,inres)=d_t_work(ind+j)
4243 end subroutine sd_verlet2_ciccotti
4245 !-----------------------------------------------------------------------------
4247 !-----------------------------------------------------------------------------
4249 subroutine inertia_tensor
4252 real(kind=8) Im(3,3),Imcp(3,3),pr(3),M_SC,&
4253 eigvec(3,3),Id(3,3),eigval(3),L(3),vp(3),vrot(3),&
4254 vpp(3,0:MAXRES),vs_p(3),pr1(3,3),&
4255 pr2(3,3),pp(3),incr(3),v(3),mag,mag2,M_PEP
4256 integer iti,inres,i,j,k,mnum,mnum1
4269 !c caulating the center of the mass of the protein
4273 if (itype(i,mnum).eq.ntyp1_molec(mnum)&
4274 .or. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) cycle
4275 if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
4276 M_PEP=M_PEP+mp(mnum)
4279 cm(j)=cm(j)+(c(j,i)+0.5d0*dc(j,i))*mp(mnum)
4289 iti=iabs(itype(i,mnum))
4290 if (iti.eq.ntyp1_molec(mnum)) cycle
4291 M_SC=M_SC+msc(iabs(iti),mnum)
4294 cm(j)=cm(j)+msc(iabs(iti),mnum)*c(j,inres)
4298 cm(j)=cm(j)/(M_SC+M_PEP)
4303 if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
4304 if (itype(i,mnum).eq.ntyp1_molec(mnum)&
4305 .or. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) cycle
4307 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
4309 Im(1,1)=Im(1,1)+mp(mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
4310 Im(1,2)=Im(1,2)-mp(mnum)*pr(1)*pr(2)
4311 Im(1,3)=Im(1,3)-mp(mnum)*pr(1)*pr(3)
4312 Im(2,3)=Im(2,3)-mp(mnum)*pr(2)*pr(3)
4313 Im(2,2)=Im(2,2)+mp(mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
4314 Im(3,3)=Im(3,3)+mp(mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
4319 iti=iabs(itype(i,mnum))
4320 if (iti.eq.ntyp1_molec(mnum)) cycle
4323 pr(j)=c(j,inres)-cm(j)
4325 Im(1,1)=Im(1,1)+msc(iabs(iti),mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
4326 Im(1,2)=Im(1,2)-msc(iabs(iti),mnum)*pr(1)*pr(2)
4327 Im(1,3)=Im(1,3)-msc(iabs(iti),mnum)*pr(1)*pr(3)
4328 Im(2,3)=Im(2,3)-msc(iabs(iti),mnum)*pr(2)*pr(3)
4329 Im(2,2)=Im(2,2)+msc(iabs(iti),mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
4330 Im(3,3)=Im(3,3)+msc(iabs(iti),mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
4335 if (itype(i,mnum).eq.ntyp1_molec(mnum)&
4336 .or. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) cycle
4337 Im(1,1)=Im(1,1)+Ip(mnum)*(1-dc_norm(1,i)*dc_norm(1,i))*&
4338 vbld(i+1)*vbld(i+1)*0.25d0
4339 Im(1,2)=Im(1,2)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(2,i))*&
4340 vbld(i+1)*vbld(i+1)*0.25d0
4341 Im(1,3)=Im(1,3)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(3,i))*&
4342 vbld(i+1)*vbld(i+1)*0.25d0
4343 Im(2,3)=Im(2,3)+Ip(mnum)*(-dc_norm(2,i)*dc_norm(3,i))*&
4344 vbld(i+1)*vbld(i+1)*0.25d0
4345 Im(2,2)=Im(2,2)+Ip(mnum)*(1-dc_norm(2,i)*dc_norm(2,i))*&
4346 vbld(i+1)*vbld(i+1)*0.25d0
4347 Im(3,3)=Im(3,3)+Ip(mnum)*(1-dc_norm(3,i)*dc_norm(3,i))*&
4348 vbld(i+1)*vbld(i+1)*0.25d0
4353 iti=iabs(itype(i,mnum))
4354 if (iti.ne.10 .and. iti.ne.ntyp1_molec(mnum).and.mnum.le.2) then
4356 Im(1,1)=Im(1,1)+Isc(iti,mnum)*(1-dc_norm(1,inres)*&
4357 dc_norm(1,inres))*vbld(inres)*vbld(inres)
4358 Im(1,2)=Im(1,2)-Isc(iti,mnum)*(dc_norm(1,inres)*&
4359 dc_norm(2,inres))*vbld(inres)*vbld(inres)
4360 Im(1,3)=Im(1,3)-Isc(iti,mnum)*(dc_norm(1,inres)*&
4361 dc_norm(3,inres))*vbld(inres)*vbld(inres)
4362 Im(2,3)=Im(2,3)-Isc(iti,mnum)*(dc_norm(2,inres)*&
4363 dc_norm(3,inres))*vbld(inres)*vbld(inres)
4364 Im(2,2)=Im(2,2)+Isc(iti,mnum)*(1-dc_norm(2,inres)*&
4365 dc_norm(2,inres))*vbld(inres)*vbld(inres)
4366 Im(3,3)=Im(3,3)+Isc(iti,mnum)*(1-dc_norm(3,inres)*&
4367 dc_norm(3,inres))*vbld(inres)*vbld(inres)
4376 !c Copng the Im matrix for the djacob subroutine
4383 !c Finding the eigenvectors and eignvalues of the inertia tensor
4384 call djacob(3,3,10000,1.0d-10,Imcp,eigvec,eigval)
4386 if (dabs(eigval(i)).gt.1.0d-15) then
4387 Id(i,i)=1.0d0/eigval(i)
4394 Imcp(i,j)=eigvec(j,i)
4400 pr1(i,j)=pr1(i,j)+Id(i,k)*Imcp(k,j)
4407 pr2(i,j)=pr2(i,j)+eigvec(i,k)*pr1(k,j)
4411 !c Calculating the total rotational velocity of the molecule
4414 vrot(i)=vrot(i)+pr2(i,j)*L(j)
4420 ! write (iout,*) itype(i+1,mnum1),itype(i,mnum)
4421 if (itype(i+1,mnum1).ne.ntyp1_molec(mnum1) &
4422 .and. itype(i,mnum).eq.ntyp1_molec(mnum) .or.&
4423 itype(i,mnum).ne.ntyp1_molec(mnum) &
4424 .and. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) cycle
4425 call vecpr(vrot(1),dc(1,i),vp)
4427 d_t(j,i)=d_t(j,i)-vp(j)
4432 if(itype(i,1).ne.10 .and.itype(i,mnum).ne.ntyp1_molec(mnum)&
4433 .and.mnum.le.2) then
4435 call vecpr(vrot(1),dc(1,inres),vp)
4437 d_t(j,inres)=d_t(j,inres)-vp(j)
4443 end subroutine inertia_tensor
4444 !------------------------------------------------------------
4445 subroutine angmom(cm,L)
4447 double precision L(3),cm(3),pr(3),vp(3),vrot(3),incr(3),v(3),&
4449 integer iti,inres,i,j,mnum,mnum1
4450 !c Calculate the angular momentum
4460 if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
4461 if (itype(i,mnum).eq.ntyp1_molec(mnum)&
4462 .or. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) cycle
4464 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
4467 v(j)=incr(j)+0.5d0*d_t(j,i)
4470 incr(j)=incr(j)+d_t(j,i)
4472 call vecpr(pr(1),v(1),vp)
4474 L(j)=L(j)+mp(mnum)*vp(j)
4478 pp(j)=0.5d0*d_t(j,i)
4480 call vecpr(pr(1),pp(1),vp)
4482 L(j)=L(j)+Ip(mnum)*vp(j)
4490 iti=iabs(itype(i,mnum))
4493 pr(j)=c(j,inres)-cm(j)
4495 if (itype(i,1).ne.10 .and.itype(i,mnum).ne.ntyp1_molec(mnum)&
4496 .and.mnum.le.2) then
4498 v(j)=incr(j)+d_t(j,inres)
4505 call vecpr(pr(1),v(1),vp)
4506 !c write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),
4507 !c & " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
4514 L(j)=L(j)+mscab*vp(j)
4516 !c write (iout,*) "L",(l(j),j=1,3)
4517 if (itype(i,1).ne.10 .and.itype(i,mnum).ne.ntyp1_molec(mnum)&
4518 .and.mnum.le.2) then
4520 v(j)=incr(j)+d_t(j,inres)
4522 call vecpr(dc(1,inres),d_t(1,inres),vp)
4524 L(j)=L(j)+Isc(iti,mnum)*vp(j)
4528 incr(j)=incr(j)+d_t(j,i)
4532 end subroutine angmom
4533 !---------------------------------------------------------------
4534 subroutine vcm_vel(vcm)
4535 double precision vcm(3),vv(3),summas,amas
4536 integer i,j,mnum,mnum1
4544 if ((mnum.ge.5).or.(mnum.eq.3))&
4545 mp(mnum)=msc(itype(i,mnum),mnum)
4547 summas=summas+mp(mnum)
4549 vcm(j)=vcm(j)+mp(mnum)*(vv(j)+0.5d0*d_t(j,i))
4553 amas=msc(iabs(itype(i,mnum)),mnum)
4557 ! amas=msc(iabs(itype(i)))
4559 ! if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
4560 if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4561 .and.(mnum.lt.4)) then
4563 vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres))
4567 vcm(j)=vcm(j)+amas*vv(j)
4571 vv(j)=vv(j)+d_t(j,i)
4574 !c write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas
4576 vcm(j)=vcm(j)/summas
4579 end subroutine vcm_vel
4581 subroutine inertia_tensor
4583 ! Calculating the intertia tensor for the entire protein in order to
4584 ! remove the perpendicular components of velocity matrix which cause
4585 ! the molecule to rotate.
4588 ! implicit real*8 (a-h,o-z)
4589 ! include 'DIMENSIONS'
4590 ! include 'COMMON.CONTROL'
4591 ! include 'COMMON.VAR'
4592 ! include 'COMMON.MD'
4593 ! include 'COMMON.CHAIN'
4594 ! include 'COMMON.DERIV'
4595 ! include 'COMMON.GEO'
4596 ! include 'COMMON.LOCAL'
4597 ! include 'COMMON.INTERACT'
4598 ! include 'COMMON.IOUNITS'
4599 ! include 'COMMON.NAMES'
4601 real(kind=8),dimension(3,3) :: Im,Imcp,eigvec,Id
4602 real(kind=8),dimension(3) :: pr,eigval,L,vp,vrot
4603 real(kind=8) :: M_SC,mag,mag2,M_PEP
4604 real(kind=8),dimension(3,0:nres) :: vpp !(3,0:MAXRES)
4605 real(kind=8),dimension(3) :: vs_p,pp,incr,v
4606 real(kind=8),dimension(3,3) :: pr1,pr2
4608 !el common /gucio/ cm
4609 integer :: iti,inres,i,j,k,mnum
4620 ! calculating the center of the mass of the protein
4624 if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
4625 write(iout,*) "WTF",itype(i,mnum),i,mnum,mp(mnum)
4626 ! if (itype(i,mnum).eq.ntyp1_molec(mnum)) cycle
4627 M_PEP=M_PEP+mp(mnum)
4630 cm(j)=cm(j)+(c(j,i)+0.5d0*dc(j,i))*mp(mnum)
4639 if (mnum.ge.5) cycle
4640 iti=iabs(itype(i,mnum))
4641 M_SC=M_SC+msc(iabs(iti),mnum)
4644 cm(j)=cm(j)+msc(iabs(iti),mnum)*c(j,inres)
4648 cm(j)=cm(j)/(M_SC+M_PEP)
4650 ! write(iout,*) "Center of mass:",cm
4653 if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
4655 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
4657 Im(1,1)=Im(1,1)+mp(mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
4658 Im(1,2)=Im(1,2)-mp(mnum)*pr(1)*pr(2)
4659 Im(1,3)=Im(1,3)-mp(mnum)*pr(1)*pr(3)
4660 Im(2,3)=Im(2,3)-mp(mnum)*pr(2)*pr(3)
4661 Im(2,2)=Im(2,2)+mp(mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
4662 Im(3,3)=Im(3,3)+mp(mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
4665 ! write(iout,*) "The angular momentum before msc add"
4667 ! write (iout,*) (Im(i,j),j=1,3)
4672 iti=iabs(itype(i,mnum))
4673 if (mnum.ge.5) cycle
4676 pr(j)=c(j,inres)-cm(j)
4678 Im(1,1)=Im(1,1)+msc(iabs(iti),mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
4679 Im(1,2)=Im(1,2)-msc(iabs(iti),mnum)*pr(1)*pr(2)
4680 Im(1,3)=Im(1,3)-msc(iabs(iti),mnum)*pr(1)*pr(3)
4681 Im(2,3)=Im(2,3)-msc(iabs(iti),mnum)*pr(2)*pr(3)
4682 Im(2,2)=Im(2,2)+msc(iabs(iti),mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
4683 Im(3,3)=Im(3,3)+msc(iabs(iti),mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
4685 ! write(iout,*) "The angular momentum before Ip add"
4687 ! write (iout,*) (Im(i,j),j=1,3)
4692 Im(1,1)=Im(1,1)+Ip(mnum)*(1-dc_norm(1,i)*dc_norm(1,i))* &
4693 vbld(i+1)*vbld(i+1)*0.25d0
4694 Im(1,2)=Im(1,2)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(2,i))* &
4695 vbld(i+1)*vbld(i+1)*0.25d0
4696 Im(1,3)=Im(1,3)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(3,i))* &
4697 vbld(i+1)*vbld(i+1)*0.25d0
4698 Im(2,3)=Im(2,3)+Ip(mnum)*(-dc_norm(2,i)*dc_norm(3,i))* &
4699 vbld(i+1)*vbld(i+1)*0.25d0
4700 Im(2,2)=Im(2,2)+Ip(mnum)*(1-dc_norm(2,i)*dc_norm(2,i))* &
4701 vbld(i+1)*vbld(i+1)*0.25d0
4702 Im(3,3)=Im(3,3)+Ip(mnum)*(1-dc_norm(3,i)*dc_norm(3,i))* &
4703 vbld(i+1)*vbld(i+1)*0.25d0
4705 ! write(iout,*) "The angular momentum before Isc add"
4707 ! write (iout,*) (Im(i,j),j=1,3)
4713 ! if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
4714 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4715 .and.(mnum.lt.4)) then
4716 iti=iabs(itype(i,mnum))
4718 Im(1,1)=Im(1,1)+Isc(iti,mnum)*(1-dc_norm(1,inres)* &
4719 dc_norm(1,inres))*vbld(inres)*vbld(inres)
4720 Im(1,2)=Im(1,2)-Isc(iti,mnum)*(dc_norm(1,inres)* &
4721 dc_norm(2,inres))*vbld(inres)*vbld(inres)
4722 Im(1,3)=Im(1,3)-Isc(iti,mnum)*(dc_norm(1,inres)* &
4723 dc_norm(3,inres))*vbld(inres)*vbld(inres)
4724 Im(2,3)=Im(2,3)-Isc(iti,mnum)*(dc_norm(2,inres)* &
4725 dc_norm(3,inres))*vbld(inres)*vbld(inres)
4726 Im(2,2)=Im(2,2)+Isc(iti,mnum)*(1-dc_norm(2,inres)* &
4727 dc_norm(2,inres))*vbld(inres)*vbld(inres)
4728 Im(3,3)=Im(3,3)+Isc(iti,mnum)*(1-dc_norm(3,inres)* &
4729 dc_norm(3,inres))*vbld(inres)*vbld(inres)
4733 ! write(iout,*) "The angular momentum before agnom:"
4735 ! write (iout,*) (Im(i,j),j=1,3)
4739 ! write(iout,*) "The angular momentum before adjustment:"
4740 ! write(iout,*) (L(j),j=1,3)
4742 ! write (iout,*) (Im(i,j),j=1,3)
4748 ! Copying the Im matrix for the djacob subroutine
4756 ! Finding the eigenvectors and eignvalues of the inertia tensor
4757 call djacob(3,3,10000,1.0d-10,Imcp,eigvec,eigval)
4758 ! write (iout,*) "Eigenvalues & Eigenvectors"
4759 ! write (iout,'(5x,3f10.5)') (eigval(i),i=1,3)
4762 ! write (iout,'(i5,3f10.5)') i,(eigvec(i,j),j=1,3)
4764 ! Constructing the diagonalized matrix
4766 if (dabs(eigval(i)).gt.1.0d-15) then
4767 Id(i,i)=1.0d0/eigval(i)
4774 Imcp(i,j)=eigvec(j,i)
4780 pr1(i,j)=pr1(i,j)+Id(i,k)*Imcp(k,j)
4787 pr2(i,j)=pr2(i,j)+eigvec(i,k)*pr1(k,j)
4791 ! Calculating the total rotational velocity of the molecule
4794 vrot(i)=vrot(i)+pr2(i,j)*L(j)
4797 ! Resetting the velocities
4799 call vecpr(vrot(1),dc(1,i),vp)
4801 ! print *,"HERE2",d_t(j,i),vp(j)
4802 d_t(j,i)=d_t(j,i)-vp(j)
4803 ! print *,"HERE2",d_t(j,i),vp(j)
4808 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4809 .and.(mnum.lt.4)) then
4810 ! if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
4811 ! if(itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
4813 call vecpr(vrot(1),dc(1,inres),vp)
4815 d_t(j,inres)=d_t(j,inres)-vp(j)
4820 ! write(iout,*) "The angular momentum after adjustment:"
4821 ! write(iout,*) (L(j),j=1,3)
4824 end subroutine inertia_tensor
4825 !-----------------------------------------------------------------------------
4826 subroutine angmom(cm,L)
4829 ! implicit real*8 (a-h,o-z)
4830 ! include 'DIMENSIONS'
4831 ! include 'COMMON.CONTROL'
4832 ! include 'COMMON.VAR'
4833 ! include 'COMMON.MD'
4834 ! include 'COMMON.CHAIN'
4835 ! include 'COMMON.DERIV'
4836 ! include 'COMMON.GEO'
4837 ! include 'COMMON.LOCAL'
4838 ! include 'COMMON.INTERACT'
4839 ! include 'COMMON.IOUNITS'
4840 ! include 'COMMON.NAMES'
4841 real(kind=8) :: mscab
4842 real(kind=8),dimension(3) :: L,cm,pr,vp,vrot,incr,v,pp
4843 integer :: iti,inres,i,j,mnum
4844 ! Calculate the angular momentum
4853 if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
4855 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
4858 v(j)=incr(j)+0.5d0*d_t(j,i)
4861 incr(j)=incr(j)+d_t(j,i)
4863 call vecpr(pr(1),v(1),vp)
4865 L(j)=L(j)+mp(mnum)*vp(j)
4866 ! print *,"HERE3",J,i,L(j),mp(mnum),Ip(mnum),mnum
4870 pp(j)=0.5d0*d_t(j,i)
4872 call vecpr(pr(1),pp(1),vp)
4875 L(j)=L(j)+Ip(mnum)*vp(j)
4883 iti=iabs(itype(i,mnum))
4891 pr(j)=c(j,inres)-cm(j)
4894 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4895 .and.(mnum.lt.4)) then
4897 v(j)=incr(j)+d_t(j,inres)
4904 ! print *,i,pr,"pr",v
4905 call vecpr(pr(1),v(1),vp)
4906 ! write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),&
4907 ! " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
4909 L(j)=L(j)+mscab*vp(j)
4911 ! write (iout,*) "L",(l(j),j=1,3)
4912 ! print *,"L",(l(j),j=1,3),i,vp(1)
4914 if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4915 .and.(mnum.lt.4)) then
4917 v(j)=incr(j)+d_t(j,inres)
4919 call vecpr(dc(1,inres),d_t(1,inres),vp)
4921 L(j)=L(j)+Isc(iti,mnum)*vp(j)
4925 incr(j)=incr(j)+d_t(j,i)
4929 end subroutine angmom
4930 !-----------------------------------------------------------------------------
4931 subroutine vcm_vel(vcm)
4934 ! implicit real*8 (a-h,o-z)
4935 ! include 'DIMENSIONS'
4936 ! include 'COMMON.VAR'
4937 ! include 'COMMON.MD'
4938 ! include 'COMMON.CHAIN'
4939 ! include 'COMMON.DERIV'
4940 ! include 'COMMON.GEO'
4941 ! include 'COMMON.LOCAL'
4942 ! include 'COMMON.INTERACT'
4943 ! include 'COMMON.IOUNITS'
4944 real(kind=8),dimension(3) :: vcm,vv
4945 real(kind=8) :: summas,amas
4955 if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
4957 summas=summas+mp(mnum)
4959 vcm(j)=vcm(j)+mp(mnum)*(vv(j)+0.5d0*d_t(j,i))
4960 ! print *,i,j,vv(j),d_t(j,i)
4964 amas=msc(iabs(itype(i,mnum)),mnum)
4969 if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4970 .and.(mnum.lt.4)) then
4972 vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres))
4976 vcm(j)=vcm(j)+amas*vv(j)
4980 vv(j)=vv(j)+d_t(j,i)
4983 ! write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas
4985 vcm(j)=vcm(j)/summas
4988 end subroutine vcm_vel
4990 !-----------------------------------------------------------------------------
4992 !-----------------------------------------------------------------------------
4994 ! RATTLE algorithm for velocity Verlet - step 1, UNRES
4998 ! implicit real*8 (a-h,o-z)
4999 ! include 'DIMENSIONS'
5001 ! include 'COMMON.CONTROL'
5002 ! include 'COMMON.VAR'
5003 ! include 'COMMON.MD'
5005 ! include 'COMMON.LANGEVIN'
5007 ! include 'COMMON.LANGEVIN.lang0'
5009 ! include 'COMMON.CHAIN'
5010 ! include 'COMMON.DERIV'
5011 ! include 'COMMON.GEO'
5012 ! include 'COMMON.LOCAL'
5013 ! include 'COMMON.INTERACT'
5014 ! include 'COMMON.IOUNITS'
5015 ! include 'COMMON.NAMES'
5016 ! include 'COMMON.TIME1'
5017 !el real(kind=8) :: gginv(2*nres,2*nres),&
5018 !el gdc(3,2*nres,2*nres)
5019 real(kind=8) :: dC_uncor(3,2*nres) !,&
5020 !el real(kind=8) :: Cmat(2*nres,2*nres)
5021 real(kind=8) :: x(2*nres),xcorr(3,2*nres) !maxres2=2*maxres
5022 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
5023 !el common /przechowalnia/ nbond
5024 integer :: max_rattle = 5
5025 logical :: lprn = .false., lprn1 = .false., not_done
5026 real(kind=8) :: tol_rattle = 1.0d-5
5028 integer :: ii,i,j,jj,l,ind,ind1,nres2
5031 !el /common/ przechowalnia
5033 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
5034 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
5035 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
5037 if (lprn) write (iout,*) "RATTLE1"
5041 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5042 .and.(mnum.lt.4)) nbond=nbond+1
5044 ! Make a folded form of the Ginv-matrix
5057 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
5062 if (itype(k,1).ne.10 .and. itype(k,mnum).ne.ntyp1_molec(mnum)&
5063 .and.(mnum.lt.4)) then
5067 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
5075 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5086 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
5090 if (itype(k,1).ne.10) then
5094 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
5102 write (iout,*) "Matrix GGinv"
5103 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
5109 if (iter.gt.max_rattle) then
5110 write (iout,*) "Error - too many iterations in RATTLE."
5113 ! Calculate the matrix C = GG**(-1) dC_old o dC
5118 dC_uncor(j,ind1)=dC(j,i)
5122 if (itype(i,1).ne.10) then
5125 dC_uncor(j,ind1)=dC(j,i+nres)
5134 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
5138 if (itype(k,1).ne.10) then
5141 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
5146 ! Calculate deviations from standard virtual-bond lengths
5150 x(ind)=vbld(i+1)**2-vbl**2
5153 if (itype(i,1).ne.10) then
5155 x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
5159 write (iout,*) "Coordinates and violations"
5161 write(iout,'(i5,3f10.5,5x,e15.5)') &
5162 i,(dC_uncor(j,i),j=1,3),x(i)
5164 write (iout,*) "Velocities and violations"
5168 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5169 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
5173 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5174 .and.(mnum.lt.4)) then
5177 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5178 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
5179 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
5182 ! write (iout,*) "gdc"
5184 ! write (iout,*) "i",i
5186 ! write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
5192 if (dabs(x(i)).gt.xmax) then
5196 if (xmax.lt.tol_rattle) then
5200 ! Calculate the matrix of the system of equations
5205 Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
5210 write (iout,*) "Matrix Cmat"
5211 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
5213 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
5214 ! Add constraint term to positions
5221 xx = xx+x(ii)*gdc(j,ind,ii)
5225 d_t_new(j,i)=d_t_new(j,i)-xx/d_time
5229 if (itype(i,1).ne.10) then
5234 xx = xx+x(ii)*gdc(j,ind,ii)
5237 dC(j,i+nres)=dC(j,i+nres)-xx
5238 d_t_new(j,i+nres)=d_t_new(j,i+nres)-xx/d_time
5242 ! Rebuild the chain using the new coordinates
5243 call chainbuild_cart
5245 write (iout,*) "New coordinates, Lagrange multipliers,",&
5246 " and differences between actual and standard bond lengths"
5250 xx=vbld(i+1)**2-vbl**2
5251 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
5252 i,(dC(j,i),j=1,3),x(ind),xx
5256 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5259 xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
5260 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
5261 i,(dC(j,i+nres),j=1,3),x(ind),xx
5264 write (iout,*) "Velocities and violations"
5268 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5269 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
5272 if (itype(i,1).ne.10) then
5274 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5275 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
5276 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
5283 10 write (iout,*) "Error - singularity in solving the system",&
5284 " of equations for Lagrange multipliers."
5288 "RATTLE inactive; use -DRATTLE switch at compile time."
5291 end subroutine rattle1
5292 !-----------------------------------------------------------------------------
5294 ! RATTLE algorithm for velocity Verlet - step 2, UNRES
5298 ! implicit real*8 (a-h,o-z)
5299 ! include 'DIMENSIONS'
5301 ! include 'COMMON.CONTROL'
5302 ! include 'COMMON.VAR'
5303 ! include 'COMMON.MD'
5305 ! include 'COMMON.LANGEVIN'
5307 ! include 'COMMON.LANGEVIN.lang0'
5309 ! include 'COMMON.CHAIN'
5310 ! include 'COMMON.DERIV'
5311 ! include 'COMMON.GEO'
5312 ! include 'COMMON.LOCAL'
5313 ! include 'COMMON.INTERACT'
5314 ! include 'COMMON.IOUNITS'
5315 ! include 'COMMON.NAMES'
5316 ! include 'COMMON.TIME1'
5317 !el real(kind=8) :: gginv(2*nres,2*nres),&
5318 !el gdc(3,2*nres,2*nres)
5319 real(kind=8) :: dC_uncor(3,2*nres) !,&
5320 !el Cmat(2*nres,2*nres)
5321 real(kind=8) :: x(2*nres) !maxres2=2*maxres
5322 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
5323 !el common /przechowalnia/ nbond
5324 integer :: max_rattle = 5
5325 logical :: lprn = .false., lprn1 = .false., not_done
5326 real(kind=8) :: tol_rattle = 1.0d-5
5330 !el /common/ przechowalnia
5331 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
5332 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
5333 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
5335 if (lprn) write (iout,*) "RATTLE2"
5336 if (lprn) write (iout,*) "Velocity correction"
5337 ! Calculate the matrix G dC
5343 gdc(j,i,ind)=GGinv(i,ind)*dC(j,k)
5348 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5349 .and.(mnum.lt.4)) then
5352 gdc(j,i,ind)=GGinv(i,ind)*dC(j,k+nres)
5358 ! write (iout,*) "gdc"
5360 ! write (iout,*) "i",i
5362 ! write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
5366 ! Calculate the matrix of the system of equations
5373 Cmat(ind,j)=Cmat(ind,j)+dC(k,i)*gdc(k,ind,j)
5379 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5380 .and.(mnum.lt.4)) then
5385 Cmat(ind,j)=Cmat(ind,j)+dC(k,i+nres)*gdc(k,ind,j)
5390 ! Calculate the scalar product dC o d_t_new
5394 x(ind)=scalar(d_t(1,i),dC(1,i))
5398 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5399 .and.(mnum.lt.4)) then
5401 x(ind)=scalar(d_t(1,i+nres),dC(1,i+nres))
5405 write (iout,*) "Velocities and violations"
5409 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5410 i,ind,(d_t(j,i),j=1,3),x(ind)
5414 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5415 .and.(mnum.lt.4)) then
5417 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5418 i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind)
5424 if (dabs(x(i)).gt.xmax) then
5428 if (xmax.lt.tol_rattle) then
5433 write (iout,*) "Matrix Cmat"
5434 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
5436 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
5437 ! Add constraint term to velocities
5444 xx = xx+x(ii)*gdc(j,ind,ii)
5446 d_t(j,i)=d_t(j,i)-xx
5451 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5452 .and.(mnum.lt.4)) then
5457 xx = xx+x(ii)*gdc(j,ind,ii)
5459 d_t(j,i+nres)=d_t(j,i+nres)-xx
5465 "New velocities, Lagrange multipliers violations"
5469 if (lprn) write (iout,'(2i5,3f10.5,5x,2e15.5)') &
5470 i,ind,(d_t(j,i),j=1,3),x(ind),scalar(d_t(1,i),dC(1,i))
5474 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5477 write (iout,'(2i5,3f10.5,5x,2e15.5)') &
5478 i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind),&
5479 scalar(d_t(1,i+nres),dC(1,i+nres))
5485 10 write (iout,*) "Error - singularity in solving the system",&
5486 " of equations for Lagrange multipliers."
5490 "RATTLE inactive; use -DRATTLE option at compile time."
5493 end subroutine rattle2
5494 !-----------------------------------------------------------------------------
5495 subroutine rattle_brown
5496 ! RATTLE/LINCS algorithm for Brownian dynamics, UNRES
5500 ! implicit real*8 (a-h,o-z)
5501 ! include 'DIMENSIONS'
5503 ! include 'COMMON.CONTROL'
5504 ! include 'COMMON.VAR'
5505 ! include 'COMMON.MD'
5507 ! include 'COMMON.LANGEVIN'
5509 ! include 'COMMON.LANGEVIN.lang0'
5511 ! include 'COMMON.CHAIN'
5512 ! include 'COMMON.DERIV'
5513 ! include 'COMMON.GEO'
5514 ! include 'COMMON.LOCAL'
5515 ! include 'COMMON.INTERACT'
5516 ! include 'COMMON.IOUNITS'
5517 ! include 'COMMON.NAMES'
5518 ! include 'COMMON.TIME1'
5519 !el real(kind=8) :: gginv(2*nres,2*nres),&
5520 !el gdc(3,2*nres,2*nres)
5521 real(kind=8) :: dC_uncor(3,2*nres) !,&
5522 !el real(kind=8) :: Cmat(2*nres,2*nres)
5523 real(kind=8) :: x(2*nres) !maxres2=2*maxres
5524 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
5525 !el common /przechowalnia/ nbond
5526 integer :: max_rattle = 5
5527 logical :: lprn = .false., lprn1 = .false., not_done
5528 real(kind=8) :: tol_rattle = 1.0d-5
5532 !el /common/ przechowalnia
5533 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
5534 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
5535 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
5538 if (lprn) write (iout,*) "RATTLE_BROWN"
5541 if (itype(i,1).ne.10) nbond=nbond+1
5543 ! Make a folded form of the Ginv-matrix
5556 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
5560 if (itype(k,1).ne.10) then
5564 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
5571 if (itype(i,1).ne.10) then
5581 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
5585 if (itype(k,1).ne.10) then
5589 if (j.eq.1 .and. l.eq.1)GGinv(ii,jj)=fricmat(ind,ind1)
5597 write (iout,*) "Matrix GGinv"
5598 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
5604 if (iter.gt.max_rattle) then
5605 write (iout,*) "Error - too many iterations in RATTLE."
5608 ! Calculate the matrix C = GG**(-1) dC_old o dC
5613 dC_uncor(j,ind1)=dC(j,i)
5617 if (itype(i,1).ne.10) then
5620 dC_uncor(j,ind1)=dC(j,i+nres)
5629 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
5633 if (itype(k,1).ne.10) then
5636 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
5641 ! Calculate deviations from standard virtual-bond lengths
5645 x(ind)=vbld(i+1)**2-vbl**2
5648 if (itype(i,1).ne.10) then
5650 x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
5654 write (iout,*) "Coordinates and violations"
5656 write(iout,'(i5,3f10.5,5x,e15.5)') &
5657 i,(dC_uncor(j,i),j=1,3),x(i)
5659 write (iout,*) "Velocities and violations"
5663 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5664 i,ind,(d_t(j,i),j=1,3),scalar(d_t(1,i),dC_old(1,i))
5667 if (itype(i,1).ne.10) then
5669 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5670 i+nres,ind,(d_t(j,i+nres),j=1,3),&
5671 scalar(d_t(1,i+nres),dC_old(1,i+nres))
5674 write (iout,*) "gdc"
5676 write (iout,*) "i",i
5678 write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
5684 if (dabs(x(i)).gt.xmax) then
5688 if (xmax.lt.tol_rattle) then
5692 ! Calculate the matrix of the system of equations
5697 Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
5702 write (iout,*) "Matrix Cmat"
5703 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
5705 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
5706 ! Add constraint term to positions
5713 xx = xx+x(ii)*gdc(j,ind,ii)
5716 d_t(j,i)=d_t(j,i)+xx/d_time
5721 if (itype(i,1).ne.10) then
5726 xx = xx+x(ii)*gdc(j,ind,ii)
5729 d_t(j,i+nres)=d_t(j,i+nres)+xx/d_time
5730 dC(j,i+nres)=dC(j,i+nres)+xx
5734 ! Rebuild the chain using the new coordinates
5735 call chainbuild_cart
5737 write (iout,*) "New coordinates, Lagrange multipliers,",&
5738 " and differences between actual and standard bond lengths"
5742 xx=vbld(i+1)**2-vbl**2
5743 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
5744 i,(dC(j,i),j=1,3),x(ind),xx
5747 if (itype(i,1).ne.10) then
5749 xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
5750 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
5751 i,(dC(j,i+nres),j=1,3),x(ind),xx
5754 write (iout,*) "Velocities and violations"
5758 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5759 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
5762 if (itype(i,1).ne.10) then
5764 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5765 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
5766 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
5773 10 write (iout,*) "Error - singularity in solving the system",&
5774 " of equations for Lagrange multipliers."
5778 "RATTLE inactive; use -DRATTLE option at compile time"
5781 end subroutine rattle_brown
5782 !-----------------------------------------------------------------------------
5784 !-----------------------------------------------------------------------------
5785 subroutine friction_force
5790 ! implicit real*8 (a-h,o-z)
5791 ! include 'DIMENSIONS'
5792 ! include 'COMMON.VAR'
5793 ! include 'COMMON.CHAIN'
5794 ! include 'COMMON.DERIV'
5795 ! include 'COMMON.GEO'
5796 ! include 'COMMON.LOCAL'
5797 ! include 'COMMON.INTERACT'
5798 ! include 'COMMON.MD'
5800 ! include 'COMMON.LANGEVIN'
5802 ! include 'COMMON.LANGEVIN.lang0'
5804 ! include 'COMMON.IOUNITS'
5805 !el real(kind=8),dimension(6*nres) :: gamvec !(MAXRES6) maxres6=6*maxres
5806 !el common /syfek/ gamvec
5808 integer iposc,ichain,n,innt,inct
5809 double precision rs(nres*2)
5810 real(kind=8) ::v_work(3,6*nres),vvec(2*nres)
5812 real(kind=8) :: v_work(6*nres)
5815 real(kind=8) :: vv(3),vvtot(3,nres)!,v_work(6*nres) !,&
5816 !el ginvfric(2*nres,2*nres) !maxres2=2*maxres
5817 !el common /przechowalnia/ ginvfric
5819 logical :: lprn, checkmode
5820 integer :: i,j,ind,k,nres2,nres6,mnum
5825 ! if (large) lprn=.true.
5826 ! if (large) checkmode=.true.
5828 !c Here accelerations due to friction forces are computed right after forces.
5829 d_t_work(:6*nres)=0.0d0
5831 v_work(j,1)=d_t(j,0)
5832 v_work(j,nnt)=d_t(j,0)
5836 v_work(j,i)=v_work(j,i-1)+d_t(j,i-1)
5841 if (iabs(itype(i,1)).ne.10 .and. iabs(itype(i,mnum)).ne.ntyp1_molec(mnum).and.mnum.lt.3) then
5843 v_work(j,i+nres)=v_work(j,i)+d_t(j,i+nres)
5848 write (iout,*) "v_work"
5850 write (iout,'(i5,3f10.5)') i,(v_work(j,i),j=1,3)
5856 n=dimen_chain(ichain)
5857 iposc=iposd_chain(ichain)
5858 !c write (iout,*) "friction_force j",j," ichain",ichain,
5859 !c & " n",n," iposc",iposc,iposc+n-1
5860 innt=chain_border(1,ichain)
5861 inct=chain_border(2,ichain)
5863 !c innt=chain_border(1,1)
5864 !c inct=chain_border(2,1)
5866 vvec(ind+1)=v_work(j,i)
5868 ! if (iabs(itype(i)).ne.10) then
5869 if (iabs(itype(i,1)).ne.10 .and. iabs(itype(i,mnum)).ne.ntyp1_molec(mnum).and.mnum.lt.3) then
5870 vvec(ind+1)=v_work(j,i+nres)
5875 write (iout,*) "vvec ind",ind," n",n
5876 write (iout,'(f10.5)') (vvec(i),i=iposc,ind)
5878 !c write (iout,*) "chain",i," ind",ind," n",n
5886 ! if (large) print *,"before fivediagmult"
5887 call fivediagmult(n,DMfric(iposc),DU1fric(iposc),&
5888 DU2fric(iposc),vvec(iposc),rs)
5889 ! if (large) print *,"after fivediagmult"
5893 time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
5895 time_fricmatmult=time_fricmatmult+tcpu()-time01
5900 write (iout,'(f10.5)') (rs(i),i=1,n)
5902 do i=iposc,iposc+n-1
5903 ! write (iout,*) "ichain",ichain," i",i," j",j,&
5904 ! "index",3*(i-1)+j,"rs",rs(i-iposc+1)
5905 fric_work(3*(i-1)+j)=-rs(i-iposc+1)
5910 write (iout,*) "Vector fric_work dimen3",dimen3
5911 write (iout,'(3f10.5)') (fric_work(j),j=1,dimen3)
5914 if(.not.allocated(gamvec)) allocate(gamvec(nres6)) !(MAXRES6)
5915 if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
5923 d_t_work(j)=d_t(j,0)
5928 d_t_work(ind+j)=d_t(j,i)
5934 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5935 .and.(mnum.lt.4)) then
5937 d_t_work(ind+j)=d_t(j,i+nres)
5943 call fricmat_mult(d_t_work,fric_work)
5945 if (.not.checkmode) return
5948 write (iout,*) "d_t_work and fric_work"
5950 write (iout,'(i3,2e15.5)') i,d_t_work(i),fric_work(i)
5954 friction(j,0)=fric_work(j)
5959 friction(j,i)=fric_work(ind+j)
5965 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5966 .and.(mnum.lt.4)) then
5967 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
5969 friction(j,i+nres)=fric_work(ind+j)
5975 write(iout,*) "Friction backbone"
5977 write(iout,'(i5,3e15.5,5x,3e15.5)') &
5978 i,(friction(j,i),j=1,3),(d_t(j,i),j=1,3)
5980 write(iout,*) "Friction side chain"
5982 write(iout,'(i5,3e15.5,5x,3e15.5)') &
5983 i,(friction(j,i+nres),j=1,3),(d_t(j,i+nres),j=1,3)
5992 vvtot(j,i)=vv(j)+0.5d0*d_t(j,i)
5993 vvtot(j,i+nres)=vv(j)+d_t(j,i+nres)
5994 vv(j)=vv(j)+d_t(j,i)
5997 write (iout,*) "vvtot backbone and sidechain"
5999 write (iout,'(i5,3e15.5,5x,3e15.5)') i,(vvtot(j,i),j=1,3),&
6000 (vvtot(j,i+nres),j=1,3)
6005 v_work(ind+j)=vvtot(j,i)
6011 v_work(ind+j)=vvtot(j,i+nres)
6015 write (iout,*) "v_work gamvec and site-based friction forces"
6017 write (iout,'(i5,3e15.5)') i,v_work(i),gamvec(i),&
6021 ! fric_work1(i)=0.0d0
6023 ! fric_work1(i)=fric_work1(i)-A(j,i)*gamvec(j)*v_work(j)
6026 ! write (iout,*) "fric_work and fric_work1"
6028 ! write (iout,'(i5,2e15.5)') i,fric_work(i),fric_work1(i)
6034 ginvfric(i,j)=ginvfric(i,j)+ginv(i,k)*fricmat(k,j)
6038 write (iout,*) "ginvfric"
6040 write (iout,'(i5,100f8.3)') i,(ginvfric(i,j),j=1,dimen)
6042 write (iout,*) "symmetry check"
6045 write (iout,*) i,j,ginvfric(i,j)-ginvfric(j,i)
6051 end subroutine friction_force
6052 !-----------------------------------------------------------------------------
6053 subroutine setup_fricmat
6057 use control_data, only:time_Bcast
6058 use control, only:tcpu
6060 ! implicit real*8 (a-h,o-z)
6064 real(kind=8) :: time00
6066 ! include 'DIMENSIONS'
6067 ! include 'COMMON.VAR'
6068 ! include 'COMMON.CHAIN'
6069 ! include 'COMMON.DERIV'
6070 ! include 'COMMON.GEO'
6071 ! include 'COMMON.LOCAL'
6072 ! include 'COMMON.INTERACT'
6073 ! include 'COMMON.MD'
6074 ! include 'COMMON.SETUP'
6075 ! include 'COMMON.TIME1'
6076 ! integer licznik /0/
6079 ! include 'COMMON.LANGEVIN'
6081 ! include 'COMMON.LANGEVIN.lang0'
6083 ! include 'COMMON.IOUNITS'
6085 integer :: i,j,ind,ind1,m,ichain,innt,inct
6086 logical :: lprn = .false.
6087 real(kind=8) :: dtdi !el ,gamvec(2*nres)
6088 !el real(kind=8),dimension(2*nres,2*nres) :: ginvfric,fcopy
6089 ! real(kind=8),allocatable,dimension(:,:) :: fcopy
6090 !el real(kind=8),dimension(2*nres*(2*nres+1)/2) :: Ghalf !(mmaxres2) (mmaxres2=(maxres2*(maxres2+1)/2))
6091 !el common /syfek/ gamvec
6092 real(kind=8) :: work(8*2*nres)
6093 integer :: iwork(2*nres)
6094 !el common /przechowalnia/ ginvfric,Ghalf,fcopy
6095 integer :: ii,iti,k,l,nzero,nres2,nres6,ierr,mnum
6100 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
6101 if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
6102 if (fg_rank.ne.king) goto 10
6108 if(.not.allocated(gamvec)) allocate(gamvec(nres2)) !(MAXRES2)
6110 if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
6111 if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
6112 !el allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
6113 if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !maxres2=2*maxres
6115 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
6118 if (.not.allocated(DMfric)) allocate(DMfric(nres2))
6119 if (.not.allocated(DU1fric)) allocate(DU1fric(nres2))
6120 if (.not.allocated(DU2fric)) allocate(DU2fric(nres2))
6121 ! Load the friction coefficients corresponding to peptide groups
6126 gamvec(ind1)=gamp(mnum)
6129 if (molnum(nct).eq.5) then
6132 gamvec(ind1)=gamp(mnum)
6134 ! Load the friction coefficients corresponding to side chains
6138 gamsc(ntyp1_molec(j),j)=1.0d0
6145 gamvec(ii)=gamsc(iabs(iti),mnum)
6147 if (surfarea) call sdarea(gamvec)
6153 innt=chain_border(1,ichain)
6154 inct=chain_border(2,ichain)
6155 !c write (iout,*) "ichain",ichain," innt",innt," inct",inct
6158 DMfric(ind)=gamvec(innt-nnt+1)/4
6159 if (iabs(itype(innt,1)).eq.10.or.mnum.gt.2) then
6160 DMfric(ind)=DMfric(ind)+gamvec(m+innt-nnt+1)
6163 DMfric(ind+1)=gamvec(m+innt-nnt+1)
6166 !c write (iout,*) "DMfric init ind",ind
6170 DMfric(ind)=gamvec(i-nnt+1)/2
6171 if (iabs(itype(i,1)).eq.10.or.mnum.gt.2) then
6172 DMfric(ind)=DMfric(ind)+gamvec(m+i-nnt+1)
6175 DMfric(ind+1)=gamvec(m+i-nnt+1)
6179 !c write (iout,*) "DMfric endloop ind",ind
6180 if (inct.gt.innt) then
6181 DMfric(ind)=gamvec(inct-1-nnt+1)/4
6183 if (iabs(itype(inct,1)).eq.10.or.mnum.gt.2) then
6184 DMfric(ind)=DMfric(ind)+gamvec(inct+m-nnt+1)
6187 DMfric(ind+1)=gamvec(inct+m-nnt+1)
6191 !c write (iout,*) "DMfric end ind",ind
6197 ind=iposd_chain(ichain)
6198 innt=chain_border(1,ichain)
6199 inct=chain_border(2,ichain)
6201 if (iabs(itype(i,1)).ne.10.and.mnum.le.2) then
6204 DU1fric(ind)=gamvec(i-nnt+1)/4
6212 ind=iposd_chain(ichain)
6213 innt=chain_border(1,ichain)
6214 inct=chain_border(2,ichain)
6216 if (iabs(itype(i,1)).ne.10.and.mnum.le.2) then
6217 DU2fric(ind)=gamvec(i-nnt+1)/4
6218 DU2fric(ind+1)=0.0d0
6227 write(iout,*)"The upper part of the five-diagonal friction matrix"
6229 write (iout,'(a,i5)') 'Chain',ichain
6230 innt=iposd_chain(ichain)
6231 inct=iposd_chain(ichain)+dimen_chain(ichain)-1
6233 if (i.lt.inct-1) then
6234 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMfric(i),DU1fric(i),&
6236 else if (i.eq.inct-1) then
6237 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMfric(i),DU1fric(i)
6239 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMfric(i)
6248 ! Zeroing out fricmat
6254 ! Load the friction coefficients corresponding to peptide groups
6259 gamvec(ind1)=gamp(mnum)
6262 if (molnum(nct).eq.5) then
6265 gamvec(ind1)=gamp(mnum)
6267 ! Load the friction coefficients corresponding to side chains
6271 gamsc(ntyp1_molec(j),j)=1.0d0
6278 gamvec(ii)=gamsc(iabs(iti),mnum)
6280 if (surfarea) call sdarea(gamvec)
6282 ! write (iout,*) "Matrix A and vector gamma"
6284 ! write (iout,'(i2,$)') i
6286 ! write (iout,'(f4.1,$)') A(i,j)
6288 ! write (iout,'(f8.3)') gamvec(i)
6292 write (iout,*) "Vector gamvec"
6294 write (iout,'(i5,f10.5)') i, gamvec(i)
6298 ! The friction matrix
6303 dtdi=dtdi+A(j,k)*A(j,i)*gamvec(j)
6310 write (iout,'(//a)') "Matrix fricmat"
6311 call matout2(dimen,dimen,nres2,nres2,fricmat)
6313 if (lang.eq.2 .or. lang.eq.3) then
6314 ! Mass-scale the friction matrix if non-direct integration will be performed
6320 Ginvfric(i,j)=Ginvfric(i,j)+ &
6321 Gsqrm(i,k)*Gsqrm(l,j)*fricmat(k,l)
6326 ! Diagonalize the friction matrix
6331 Ghalf(ind)=Ginvfric(i,j)
6334 call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
6337 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
6338 " mass-scaled friction matrix"
6339 call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
6341 ! Precompute matrices for tinker stochastic integrator
6348 mt1(i,j)=mt1(i,j)+fricvec(k,i)*gsqrm(k,j)
6349 mt2(i,j)=mt2(i,j)+fricvec(k,i)*gsqrp(k,j)
6355 else if (lang.eq.4) then
6356 ! Diagonalize the friction matrix
6361 Ghalf(ind)=fricmat(i,j)
6364 call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
6367 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
6369 call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
6371 ! Determine the number of zero eigenvalues of the friction matrix
6372 nzero=max0(dimen-dimen1,0)
6373 ! do while (fricgam(nzero+1).le.1.0d-5 .and. nzero.lt.dimen)
6376 write (iout,*) "Number of zero eigenvalues:",nzero
6381 fricmat(i,j)=fricmat(i,j) &
6382 +fricvec(i,k)*fricvec(j,k)/fricgam(k)
6387 write (iout,'(//a)') "Generalized inverse of fricmat"
6388 call matout(dimen,dimen,nres6,nres6,fricmat)
6393 if (nfgtasks.gt.1) then
6394 if (fg_rank.eq.0) then
6395 ! The matching BROADCAST for fg processors is called in ERGASTULUM
6401 call MPI_Bcast(10,1,MPI_INTEGER,king,FG_COMM,IERROR)
6403 time_Bcast=time_Bcast+MPI_Wtime()-time00
6405 time_Bcast=time_Bcast+tcpu()-time00
6407 ! print *,"Processor",myrank,
6408 ! & " BROADCAST iorder in SETUP_FRICMAT"
6411 write (iout,*) "setup_fricmat licznik"!,licznik !sp
6417 ! Scatter the friction matrix
6418 call MPI_Scatterv(fricmat(1,1),nginv_counts(0),&
6419 nginv_start(0),MPI_DOUBLE_PRECISION,fcopy(1,1),&
6420 myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
6423 time_scatter=time_scatter+MPI_Wtime()-time00
6424 time_scatter_fmat=time_scatter_fmat+MPI_Wtime()-time00
6426 time_scatter=time_scatter+tcpu()-time00
6427 time_scatter_fmat=time_scatter_fmat+tcpu()-time00
6431 do j=1,2*my_ng_count
6432 fricmat(j,i)=fcopy(i,j)
6435 ! write (iout,*) "My chunk of fricmat"
6436 ! call MATOUT2(my_ng_count,dimen,maxres2,maxres2,fcopy)
6441 end subroutine setup_fricmat
6442 !-----------------------------------------------------------------------------
6443 subroutine stochastic_force(stochforcvec)
6446 use random, only:anorm_distr
6447 ! implicit real*8 (a-h,o-z)
6448 ! include 'DIMENSIONS'
6449 use control, only: tcpu
6454 ! include 'COMMON.VAR'
6455 ! include 'COMMON.CHAIN'
6456 ! include 'COMMON.DERIV'
6457 ! include 'COMMON.GEO'
6458 ! include 'COMMON.LOCAL'
6459 ! include 'COMMON.INTERACT'
6460 ! include 'COMMON.MD'
6461 ! include 'COMMON.TIME1'
6463 ! include 'COMMON.LANGEVIN'
6465 ! include 'COMMON.LANGEVIN.lang0'
6467 ! include 'COMMON.IOUNITS'
6469 real(kind=8) :: x,sig,lowb,highb
6470 real(kind=8) :: ff(3),force(3,0:2*nres),zeta2,lowb2
6471 real(kind=8) :: highb2,sig2,forcvec(6*nres),stochforcvec(6*nres)
6472 real(kind=8) :: time00
6473 logical :: lprn = .false.
6474 integer :: i,j,ind,mnum
6476 integer ichain,innt,inct,iposc
6481 stochforc(j,i)=0.0d0
6491 ! Compute the stochastic forces acting on bodies. Store in force.
6497 force(j,i)=anorm_distr(x,sig,lowb,highb)
6505 force(j,i+nres)=anorm_distr(x,sig2,lowb2,highb2)
6509 time_fsample=time_fsample+MPI_Wtime()-time00
6511 time_fsample=time_fsample+tcpu()-time00
6516 innt=chain_border(1,ichain)
6517 inct=chain_border(2,ichain)
6518 iposc=iposd_chain(ichain)
6519 !c for debugging only
6520 !c innt=chain_border(1,1)
6521 !c inct=chain_border(2,1)
6522 !c iposc=iposd_chain(1)
6523 !c write (iout,*)"stochastic_force ichain=",ichain," innt",innt,
6524 !c & " inct",inct," iposc",iposc
6526 stochforcvec(ind+j)=0.5d0*force(j,innt)
6528 if (iabs(itype(innt,molnum(innt))).eq.10.or.molnum(innt).ge.3) then
6530 stochforcvec(ind+j)=stochforcvec(ind+j)+force(j,innt+nres)
6536 stochforcvec(ind+j)=force(j,innt+nres)
6542 stochforcvec(ind+j)=0.5d0*(force(j,i)+force(j,i-1))
6544 if (iabs(itype(i,1).eq.10).or.molnum(i).ge.3) then
6546 stochforcvec(ind+j)=stochforcvec(ind+j)+force(j,i+nres)
6552 stochforcvec(ind+j)=force(j,i+nres)
6558 stochforcvec(ind+j)=0.5d0*force(j,inct-1)
6560 if (iabs(itype(inct,molnum(inct))).eq.10.or.molnum(inct).ge.3) then
6562 stochforcvec(ind+j)=stochforcvec(ind+j)+force(j,inct+nres)
6568 stochforcvec(ind+j)=force(j,inct+nres)
6572 !c write (iout,*) "chain",ichain," ind",ind
6575 write (iout,*) "stochforcvec"
6576 write (iout,'(3f10.5)') (stochforcvec(j),j=1,ind)
6579 ! Compute the stochastic forces acting on virtual-bond vectors.
6585 stochforc(j,i)=ff(j)+0.5d0*force(j,i)
6588 ff(j)=ff(j)+force(j,i)
6590 ! if (itype(i+1,1).ne.ntyp1) then
6592 if (itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
6594 stochforc(j,i)=stochforc(j,i)+force(j,i+nres+1)
6595 ff(j)=ff(j)+force(j,i+nres+1)
6600 stochforc(j,0)=ff(j)+force(j,nnt+nres)
6604 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
6605 .and.(mnum.lt.4)) then
6606 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
6608 stochforc(j,i+nres)=force(j,i+nres)
6614 stochforcvec(j)=stochforc(j,0)
6619 stochforcvec(ind+j)=stochforc(j,i)
6625 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
6626 .and.(mnum.lt.4)) then
6627 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
6629 stochforcvec(ind+j)=stochforc(j,i+nres)
6635 write (iout,*) "stochforcvec"
6637 write(iout,'(i5,e15.5)') i,stochforcvec(i)
6639 write(iout,*) "Stochastic forces backbone"
6641 write(iout,'(i5,3e15.5)') i,(stochforc(j,i),j=1,3)
6643 write(iout,*) "Stochastic forces side chain"
6645 write(iout,'(i5,3e15.5)') &
6646 i,(stochforc(j,i+nres),j=1,3)
6654 write (iout,*) i,ind
6656 forcvec(ind+j)=force(j,i)
6661 write (iout,*) i,ind
6663 forcvec(j+ind)=force(j,i+nres)
6668 write (iout,*) "forcvec"
6672 write (iout,'(2i3,2f10.5)') i,j,force(j,i),&
6679 write (iout,'(2i3,2f10.5)') i,j,force(j,i+nres),&
6688 end subroutine stochastic_force
6689 !-----------------------------------------------------------------------------
6690 subroutine sdarea(gamvec)
6692 ! Scale the friction coefficients according to solvent accessible surface areas
6693 ! Code adapted from TINKER
6697 ! implicit real*8 (a-h,o-z)
6698 ! include 'DIMENSIONS'
6699 ! include 'COMMON.CONTROL'
6700 ! include 'COMMON.VAR'
6701 ! include 'COMMON.MD'
6703 ! include 'COMMON.LANGEVIN'
6705 ! include 'COMMON.LANGEVIN.lang0'
6707 ! include 'COMMON.CHAIN'
6708 ! include 'COMMON.DERIV'
6709 ! include 'COMMON.GEO'
6710 ! include 'COMMON.LOCAL'
6711 ! include 'COMMON.INTERACT'
6712 ! include 'COMMON.IOUNITS'
6713 ! include 'COMMON.NAMES'
6714 real(kind=8),dimension(2*nres) :: radius,gamvec !(maxres2)
6715 real(kind=8),parameter :: twosix = 1.122462048309372981d0
6716 logical :: lprn = .false.
6717 real(kind=8) :: probe,area,ratio
6718 integer :: i,j,ind,iti,mnum
6720 ! determine new friction coefficients every few SD steps
6722 ! set the atomic radii to estimates of sigma values
6724 ! print *,"Entered sdarea"
6730 ! Load peptide group radii
6733 radius(i)=pstok(mnum)
6735 ! Load side chain radii
6739 radius(i+nres)=restok(iti,mnum)
6742 ! write (iout,*) "i",i," radius",radius(i)
6745 radius(i) = radius(i) / twosix
6746 if (radius(i) .ne. 0.0d0) radius(i) = radius(i) + probe
6749 ! scale atomic friction coefficients by accessible area
6751 if (lprn) write (iout,*) &
6752 "Original gammas, surface areas, scaling factors, new gammas, ",&
6753 "std's of stochastic forces"
6756 if (radius(i).gt.0.0d0) then
6757 call surfatom (i,area,radius)
6758 ratio = dmax1(area/(4.0d0*pi*radius(i)**2),1.0d-1)
6759 if (lprn) write (iout,'(i5,3f10.5,$)') &
6760 i,gamvec(ind+1),area,ratio
6763 gamvec(ind) = ratio * gamvec(ind)
6766 stdforcp(i)=stdfp(mnum)*dsqrt(gamvec(ind))
6767 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcp(i)
6771 if (radius(i+nres).gt.0.0d0) then
6772 call surfatom (i+nres,area,radius)
6773 ratio = dmax1(area/(4.0d0*pi*radius(i+nres)**2),1.0d-1)
6774 if (lprn) write (iout,'(i5,3f10.5,$)') &
6775 i,gamvec(ind+1),area,ratio
6778 gamvec(ind) = ratio * gamvec(ind)
6781 stdforcsc(i)=stdfsc(itype(i,mnum),mnum)*dsqrt(gamvec(ind))
6782 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcsc(i)
6787 end subroutine sdarea
6788 !-----------------------------------------------------------------------------
6790 !-----------------------------------------------------------------------------
6793 ! ###################################################
6794 ! ## COPYRIGHT (C) 1996 by Jay William Ponder ##
6795 ! ## All Rights Reserved ##
6796 ! ###################################################
6798 ! ################################################################
6800 ! ## subroutine surfatom -- exposed surface area of an atom ##
6802 ! ################################################################
6805 ! "surfatom" performs an analytical computation of the surface
6806 ! area of a specified atom; a simplified version of "surface"
6808 ! literature references:
6810 ! T. J. Richmond, "Solvent Accessible Surface Area and
6811 ! Excluded Volume in Proteins", Journal of Molecular Biology,
6814 ! L. Wesson and D. Eisenberg, "Atomic Solvation Parameters
6815 ! Applied to Molecular Dynamics of Proteins in Solution",
6816 ! Protein Science, 1, 227-235 (1992)
6818 ! variables and parameters:
6820 ! ir number of atom for which area is desired
6821 ! area accessible surface area of the atom
6822 ! radius radii of each of the individual atoms
6825 subroutine surfatom(ir,area,radius)
6827 ! implicit real*8 (a-h,o-z)
6828 ! include 'DIMENSIONS'
6830 ! include 'COMMON.GEO'
6831 ! include 'COMMON.IOUNITS'
6833 integer :: nsup,nstart_sup
6834 ! double precision c,dc,dc_old,d_c_work,xloc,xrot,dc_norm
6835 ! common /chain/ c(3,maxres2+2),dc(3,0:maxres2),dc_old(3,0:maxres2),
6836 ! & xloc(3,maxres),xrot(3,maxres),dc_norm(3,0:maxres2),
6837 ! & dc_work(MAXRES6),nres,nres0
6838 integer,parameter :: maxarc=300
6842 integer :: mi,ni,narc
6843 integer :: key(maxarc)
6844 integer :: intag(maxarc)
6845 integer :: intag1(maxarc)
6846 real(kind=8) :: area,arcsum
6847 real(kind=8) :: arclen,exang
6848 real(kind=8) :: delta,delta2
6849 real(kind=8) :: eps,rmove
6850 real(kind=8) :: xr,yr,zr
6851 real(kind=8) :: rr,rrsq
6852 real(kind=8) :: rplus,rminus
6853 real(kind=8) :: axx,axy,axz
6854 real(kind=8) :: ayx,ayy
6855 real(kind=8) :: azx,azy,azz
6856 real(kind=8) :: uxj,uyj,uzj
6857 real(kind=8) :: tx,ty,tz
6858 real(kind=8) :: txb,tyb,td
6859 real(kind=8) :: tr2,tr,txr,tyr
6860 real(kind=8) :: tk1,tk2
6861 real(kind=8) :: thec,the,t,tb
6862 real(kind=8) :: txk,tyk,tzk
6863 real(kind=8) :: t1,ti,tf,tt
6864 real(kind=8) :: txj,tyj,tzj
6865 real(kind=8) :: ccsq,cc,xysq
6866 real(kind=8) :: bsqk,bk,cosine
6867 real(kind=8) :: dsqj,gi,pix2
6868 real(kind=8) :: therk,dk,gk
6869 real(kind=8) :: risqk,rik
6870 real(kind=8) :: radius(2*nres) !(maxatm) (maxatm=maxres2)
6871 real(kind=8) :: ri(maxarc),risq(maxarc)
6872 real(kind=8) :: ux(maxarc),uy(maxarc),uz(maxarc)
6873 real(kind=8) :: xc(maxarc),yc(maxarc),zc(maxarc)
6874 real(kind=8) :: xc1(maxarc),yc1(maxarc),zc1(maxarc)
6875 real(kind=8) :: dsq(maxarc),bsq(maxarc)
6876 real(kind=8) :: dsq1(maxarc),bsq1(maxarc)
6877 real(kind=8) :: arci(maxarc),arcf(maxarc)
6878 real(kind=8) :: ex(maxarc),lt(maxarc),gr(maxarc)
6879 real(kind=8) :: b(maxarc),b1(maxarc),bg(maxarc)
6880 real(kind=8) :: kent(maxarc),kout(maxarc)
6881 real(kind=8) :: ther(maxarc)
6882 logical :: moved,top
6883 logical :: omit(maxarc)
6886 ! maxatm = 2*nres !maxres2 maxres2=2*maxres
6887 ! maxlight = 8*maxatm
6890 ! maxtors = 4*maxatm
6892 ! zero out the surface area for the sphere of interest
6895 ! write (2,*) "ir",ir," radius",radius(ir)
6896 if (radius(ir) .eq. 0.0d0) return
6898 ! set the overlap significance and connectivity shift
6902 delta2 = delta * delta
6907 ! store coordinates and radius of the sphere of interest
6915 ! initialize values of some counters and summations
6924 ! test each sphere to see if it overlaps the sphere of interest
6927 if (i.eq.ir .or. radius(i).eq.0.0d0) goto 30
6928 rplus = rr + radius(i)
6930 if (abs(tx) .ge. rplus) goto 30
6932 if (abs(ty) .ge. rplus) goto 30
6934 if (abs(tz) .ge. rplus) goto 30
6936 ! check for sphere overlap by testing distance against radii
6938 xysq = tx*tx + ty*ty
6939 if (xysq .lt. delta2) then
6946 if (rplus-cc .le. delta) goto 30
6947 rminus = rr - radius(i)
6949 ! check to see if sphere of interest is completely buried
6951 if (cc-abs(rminus) .le. delta) then
6952 if (rminus .le. 0.0d0) goto 170
6956 ! check for too many overlaps with sphere of interest
6958 if (io .ge. maxarc) then
6960 20 format (/,' SURFATOM -- Increase the Value of MAXARC')
6964 ! get overlap between current sphere and sphere of interest
6973 gr(io) = (ccsq+rplus*rminus) / (2.0d0*rr*b1(io))
6979 ! case where no other spheres overlap the sphere of interest
6982 area = 4.0d0 * pi * rrsq
6986 ! case where only one sphere overlaps the sphere of interest
6989 area = pix2 * (1.0d0 + gr(1))
6990 area = mod(area,4.0d0*pi) * rrsq
6994 ! case where many spheres intersect the sphere of interest;
6995 ! sort the intersecting spheres by their degree of overlap
6997 call sort2 (io,gr,key)
7000 intag(i) = intag1(k)
7009 ! get radius of each overlap circle on surface of the sphere
7014 risq(i) = rrsq - gi*gi
7015 ri(i) = sqrt(risq(i))
7016 ther(i) = 0.5d0*pi - asin(min(1.0d0,max(-1.0d0,gr(i))))
7019 ! find boundary of inaccessible area on sphere of interest
7022 if (.not. omit(k)) then
7029 ! check to see if J circle is intersecting K circle;
7030 ! get distance between circle centers and sum of radii
7033 if (omit(j)) goto 60
7034 cc = (txk*xc(j)+tyk*yc(j)+tzk*zc(j))/(bk*b(j))
7035 cc = acos(min(1.0d0,max(-1.0d0,cc)))
7036 td = therk + ther(j)
7038 ! check to see if circles enclose separate regions
7040 if (cc .ge. td) goto 60
7042 ! check for circle J completely inside circle K
7044 if (cc+ther(j) .lt. therk) goto 40
7046 ! check for circles that are essentially parallel
7048 if (cc .gt. delta) goto 50
7053 ! check to see if sphere of interest is completely buried
7056 if (pix2-cc .le. td) goto 170
7062 ! find T value of circle intersections
7065 if (omit(k)) goto 110
7080 ! rotation matrix elements
7092 if (.not. omit(j)) then
7097 ! rotate spheres so K vector colinear with z-axis
7099 uxj = txj*axx + tyj*axy - tzj*axz
7100 uyj = tyj*ayy - txj*ayx
7101 uzj = txj*azx + tyj*azy + tzj*azz
7102 cosine = min(1.0d0,max(-1.0d0,uzj/b(j)))
7103 if (acos(cosine) .lt. therk+ther(j)) then
7104 dsqj = uxj*uxj + uyj*uyj
7109 tr2 = risqk*dsqj - tb*tb
7115 ! get T values of intersection for K circle
7118 tb = min(1.0d0,max(-1.0d0,tb))
7120 if (tyb-txr .lt. 0.0d0) tk1 = pix2 - tk1
7122 tb = min(1.0d0,max(-1.0d0,tb))
7124 if (tyb+txr .lt. 0.0d0) tk2 = pix2 - tk2
7125 thec = (rrsq*uzj-gk*bg(j)) / (rik*ri(j)*b(j))
7126 if (abs(thec) .lt. 1.0d0) then
7128 else if (thec .ge. 1.0d0) then
7130 else if (thec .le. -1.0d0) then
7134 ! see if "tk1" is entry or exit point; check t=0 point;
7135 ! "ti" is exit point, "tf" is entry point
7137 cosine = min(1.0d0,max(-1.0d0, &
7138 (uzj*gk-uxj*rik)/(b(j)*rr)))
7139 if ((acos(cosine)-ther(j))*(tk2-tk1) .le. 0.0d0) then
7147 if (narc .ge. maxarc) then
7149 70 format (/,' SURFATOM -- Increase the Value',&
7153 if (tf .le. ti) then
7174 ! special case; K circle without intersections
7176 if (narc .le. 0) goto 90
7178 ! general case; sum up arclength and set connectivity code
7180 call sort2 (narc,arci,key)
7185 if (narc .gt. 1) then
7188 if (t .lt. arci(j)) then
7189 arcsum = arcsum + arci(j) - t
7190 exang = exang + ex(ni)
7192 if (jb .ge. maxarc) then
7194 80 format (/,' SURFATOM -- Increase the Value',&
7199 kent(jb) = maxarc*i + k
7201 kout(jb) = maxarc*k + i
7210 arcsum = arcsum + pix2 - t
7212 exang = exang + ex(ni)
7215 kent(jb) = maxarc*i + k
7217 kout(jb) = maxarc*k + i
7224 arclen = arclen + gr(k)*arcsum
7227 if (arclen .eq. 0.0d0) goto 170
7228 if (jb .eq. 0) goto 150
7230 ! find number of independent boundaries and check connectivity
7234 if (kout(k) .ne. 0) then
7241 if (m .eq. kent(ii)) then
7244 if (j .eq. jb) goto 150
7256 ! attempt to fix connectivity error by moving atom slightly
7260 140 format (/,' SURFATOM -- Connectivity Error at Atom',i6)
7269 ! compute the exposed surface area for the sphere of interest
7272 area = ib*pix2 + exang + arclen
7273 area = mod(area,4.0d0*pi) * rrsq
7275 ! attempt to fix negative area by moving atom slightly
7277 if (area .lt. 0.0d0) then
7280 160 format (/,' SURFATOM -- Negative Area at Atom',i6)
7291 end subroutine surfatom
7292 !----------------------------------------------------------------
7293 !----------------------------------------------------------------
7294 subroutine alloc_MD_arrays
7295 !EL Allocation of arrays used by MD module
7297 integer :: nres2,nres6
7300 !----------------------
7304 allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
7305 allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
7306 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
7307 allocate(fricvec(nres2,nres2))
7308 allocate(pfric_mat(nres2,nres2),vfric_mat(nres2,nres2))
7309 allocate(afric_mat(nres2,nres2),prand_mat(nres2,nres2))
7310 allocate(vrand_mat1(nres2,nres2),vrand_mat2(nres2,nres2)) !(MAXRES2,MAXRES2)
7311 allocate(pfric0_mat(nres2,nres2,0:maxflag_stoch))
7312 allocate(afric0_mat(nres2,nres2,0:maxflag_stoch))
7313 allocate(vfric0_mat(nres2,nres2,0:maxflag_stoch))
7314 allocate(prand0_mat(nres2,nres2,0:maxflag_stoch))
7315 allocate(vrand0_mat1(nres2,nres2,0:maxflag_stoch))
7316 allocate(vrand0_mat2(nres2,nres2,0:maxflag_stoch)) !(MAXRES2,MAXRES2,0:maxflag_stoch)
7317 allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
7319 allocate(mt1(nres2,nres2),mt2(nres2,nres2),mt3(nres2,nres2)) !(maxres2,maxres2)
7320 !----------------------
7322 ! commom.langevin.lang0
7324 allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
7325 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
7326 allocate(fricvec(nres2,nres2)) !(MAXRES2,MAXRES2)
7327 allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
7328 allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
7331 if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2))
7332 !----------------------
7333 ! commom.hairpin in CSA module
7334 !----------------------
7335 ! common.mce in MCM_MD module
7336 !----------------------
7338 ! common /mdgrad/ in module.energy
7339 ! common /back_constr/ in module.energy
7340 ! common /qmeas/ in module.energy
7343 allocate(potEcomp(0:n_ene+4)) !(0:n_ene+4)
7345 allocate(d_t(3,0:nres2),d_a(3,0:nres2),d_t_old(3,0:nres2)) !(3,0:MAXRES2)
7346 allocate(d_a_work(nres6)) !(6*MAXRES)
7348 allocate(DM(nres2),DU1(nres2),DU2(nres2))
7349 allocate(DMorig(nres2),DU1orig(nres2),DU2orig(nres2))
7350 allocate(Gvec(nres2,nres2))
7352 write (iout,*) "Before A Allocation",nfgtasks-1
7354 allocate(Gmat(nres2,nres2),A(nres2,nres2))
7355 if(.not.allocated(Ginv)) allocate(Ginv(nres2,nres2)) !in control: ergastulum
7356 allocate(Gsqrp(nres2,nres2),Gsqrm(nres2,nres2),Gvec(nres2,nres2)) !(maxres2,maxres2)
7358 allocate(Geigen(nres2)) !(maxres2)
7359 if(.not.allocated(vtot)) allocate(vtot(nres2)) !(maxres2)
7360 ! common /inertia/ in io_conf: parmread
7361 ! real(kind=8),dimension(:),allocatable :: ISC,msc !(ntyp+1)
7362 ! common /langevin/in io read_MDpar
7363 ! real(kind=8),dimension(:),allocatable :: gamsc !(ntyp1)
7364 ! real(kind=8),dimension(:),allocatable :: stdfsc !(ntyp)
7365 ! in io_conf: parmread
7366 ! real(kind=8),dimension(:),allocatable :: restok !(ntyp+1)
7367 ! common /mdpmpi/ in control: ergastulum
7368 if(.not.allocated(ng_start)) allocate(ng_start(0:nfgtasks-1))
7369 if(.not.allocated(ng_counts)) allocate(ng_counts(0:nfgtasks-1))
7370 if(.not.allocated(nginv_counts)) allocate(nginv_counts(0:nfgtasks-1)) !(0:MaxProcs-1)
7371 if(.not.allocated(nginv_start)) allocate(nginv_start(0:nfgtasks)) !(0:MaxProcs)
7372 !----------------------
7373 ! common.muca in read_muca
7374 ! common /double_muca/
7375 ! real(kind=8) :: elow,ehigh,factor,hbin,factor_min
7376 ! real(kind=8),dimension(:),allocatable :: emuca,nemuca,&
7377 ! nemuca2,hist !(4*maxres)
7378 ! common /integer_muca/
7379 ! integer :: nmuca,imtime,muca_smooth
7381 ! real(kind=8),dimension(:),allocatable :: elowi,ehighi !(maxprocs)
7382 !----------------------
7384 ! common /mdgrad/ in module.energy
7385 ! common /back_constr/ in module.energy
7386 ! common /qmeas/ in module.energy
7390 allocate(d_t_work(nres6),d_t_work_new(nres6),d_af_work(nres6))
7391 allocate(d_as_work(nres6),kinetic_force(nres6)) !(MAXRES6)
7392 allocate(d_t_new(3,0:nres2),d_a_old(3,0:nres2),d_a_short(3,0:nres2)) !,d_a !(3,0:MAXRES2)
7393 allocate(stdforcp(nres),stdforcsc(nres)) !(MAXRES)
7394 !----------------------
7396 allocate(D_ban(nres6)) !(MAXRES6) maxres6=6*maxres
7397 ! common /stochcalc/ stochforcvec
7398 allocate(stochforcvec(nres6)) !(MAXRES6) maxres6=6*maxres
7399 !----------------------
7401 end subroutine alloc_MD_arrays
7402 !-----------------------------------------------------------------------------
7403 !-----------------------------------------------------------------------------