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
145 real(kind=8),dimension(:,:),allocatable :: fcopy !(2*nres,2*nres)
147 !-----------------------------------------------------------------------------
150 !-----------------------------------------------------------------------------
152 !-----------------------------------------------------------------------------
154 !-----------------------------------------------------------------------------
155 subroutine brown_step(itime)
156 !------------------------------------------------
157 ! Perform a single Euler integration step of Brownian dynamics
158 !------------------------------------------------
159 ! implicit real*8 (a-h,o-z)
161 use control, only: tcpu
164 ! use io_conf, only:cartprint
165 ! include 'DIMENSIONS'
169 ! include 'COMMON.CONTROL'
170 ! include 'COMMON.VAR'
171 ! include 'COMMON.MD'
173 ! include 'COMMON.LANGEVIN'
175 ! include 'COMMON.LANGEVIN.lang0'
177 ! include 'COMMON.CHAIN'
178 ! include 'COMMON.DERIV'
179 ! include 'COMMON.GEO'
180 ! include 'COMMON.LOCAL'
181 ! include 'COMMON.INTERACT'
182 ! include 'COMMON.IOUNITS'
183 ! include 'COMMON.NAMES'
184 ! include 'COMMON.TIME1'
185 real(kind=8),dimension(6*nres) :: zapas !(MAXRES6) maxres6=6*maxres
186 integer :: rstcount !ilen,
188 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
189 real(kind=8),dimension(6*nres,2*nres) :: Bmat,GBmat,Tmat !(MAXRES6,MAXRES2) (maxres2=2*maxres,maxres6=6*maxres)
190 real(kind=8),dimension(2*nres,2*nres) :: Cmat_,Cinv !(maxres2,maxres2) maxres2=2*maxres
191 real(kind=8),dimension(6*nres,6*nres) :: Pmat !(maxres6,maxres6) maxres6=6*maxres
192 ! real(kind=8),dimension(:,:),allocatable :: Bmat,GBmat,Tmat !(MAXRES6,MAXRES2) (maxres2=2*maxres,maxres6=6*maxres)
193 ! real(kind=8),dimension(:,:),allocatable :: Cmat_,Cinv !(maxres2,maxres2) maxres2=2*maxres
194 ! real(kind=8),dimension(:,:),allocatable :: Pmat !(maxres6,maxres6) maxres6=6*maxres
195 real(kind=8),dimension(6*nres) :: Td !(maxres6) maxres6=6*maxres
196 real(kind=8),dimension(2*nres) :: ppvec !(maxres2) maxres2=2*maxres
197 !el common /stochcalc/ stochforcvec
198 !el real(kind=8),dimension(3) :: cm !el
199 !el common /gucio/ cm
201 logical :: lprn = .false.,lprn1 = .false.
202 integer :: maxiter = 5
203 real(kind=8) :: difftol = 1.0d-5
204 real(kind=8) :: xx,diffmax,blen2,diffbond,tt0
205 integer :: i,j,nbond,k,ind,ind1,iter
206 integer :: nres2,nres6
210 ! if (.not.allocated(Bmat)) allocate(Bmat(nres6,nres2))
211 ! if (.not.allocated(GBmat)) allocate (GBmat(nres6,nres2))
212 ! if (.not.allocated(Tmat)) allocate (Tmat(nres6,nres2))
213 ! if (.not.allocated(Cmat_)) allocate(Cmat_(nres2,nres2))
214 ! if (.not.allocated(Cinv)) allocate (Cinv(nres2,nres2))
215 ! if (.not.allocated(Pmat)) allocate(Pmat(6*nres,6*nres))
217 if (.not.allocated(stochforcvec)) allocate(stochforcvec(nres6)) !(MAXRES6) maxres6=6*maxres
221 if (itype(i,1).ne.10) nbond=nbond+1
225 write (iout,*) "Generalized inverse of fricmat"
226 call matout(dimen,dimen,nres6,nres6,fricmat)
238 Bmat(ind+j,ind1)=dC_norm(j,i)
243 if (itype(i,1).ne.10) then
246 Bmat(ind+j,ind1)=dC_norm(j,i+nres)
252 write (iout,*) "Matrix Bmat"
253 call MATOUT(nbond,dimen,nres6,nres6,Bmat)
259 GBmat(i,j)=GBmat(i,j)+fricmat(i,k)*Bmat(k,j)
264 write (iout,*) "Matrix GBmat"
265 call MATOUT(nbond,dimen,nres6,nres2,Gbmat)
271 Cmat_(i,j)=Cmat_(i,j)+Bmat(k,i)*GBmat(k,j)
276 write (iout,*) "Matrix Cmat"
277 call MATOUT(nbond,nbond,nres2,nres2,Cmat_)
279 call matinvert(nbond,nres2,Cmat_,Cinv,osob)
281 write (iout,*) "Matrix Cinv"
282 call MATOUT(nbond,nbond,nres2,nres2,Cinv)
288 Tmat(i,j)=Tmat(i,j)+GBmat(i,k)*Cinv(k,j)
293 write (iout,*) "Matrix Tmat"
294 call MATOUT(nbond,dimen,nres6,nres2,Tmat)
304 Pmat(i,j)=Pmat(i,j)-Tmat(i,k)*Bmat(j,k)
309 write (iout,*) "Matrix Pmat"
310 call MATOUT(dimen,dimen,nres6,nres6,Pmat)
317 Td(i)=Td(i)+vbl*Tmat(i,ind)
320 if (itype(k,1).ne.10) then
322 Td(i)=Td(i)+vbldsc0(1,itype(k,1))*Tmat(i,ind)
327 write (iout,*) "Vector Td"
329 write (iout,'(i5,f10.5)') i,Td(i)
332 call stochastic_force(stochforcvec)
334 write (iout,*) "stochforcvec"
336 write (iout,*) i,stochforcvec(i)
340 zapas(j)=-gcart(j,0)+stochforcvec(j)
342 dC_work(j)=dC_old(j,0)
348 zapas(ind)=-gcart(j,i)+stochforcvec(ind)
349 dC_work(ind)=dC_old(j,i)
353 if (itype(i,1).ne.10) then
356 zapas(ind)=-gxcart(j,i)+stochforcvec(ind)
357 dC_work(ind)=dC_old(j,i+nres)
363 write (iout,*) "Initial d_t_work"
365 write (iout,*) i,d_t_work(i)
372 d_t_work(i)=d_t_work(i)+fricmat(i,j)*zapas(j)
379 zapas(i)=zapas(i)+Pmat(i,j)*(dC_work(j)+d_t_work(j)*d_time)
383 write (iout,*) "Final d_t_work and zapas"
385 write (iout,*) i,d_t_work(i),zapas(i)
399 dc_work(ind+j)=dc(j,i)
405 d_t(j,i+nres)=d_t_work(ind+j)
406 dc(j,i+nres)=zapas(ind+j)
407 dc_work(ind+j)=dc(j,i+nres)
413 write (iout,*) "Before correction for rotational lengthening"
414 write (iout,*) "New coordinates",&
415 " and differences between actual and standard bond lengths"
420 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
424 if (itype(i,1).ne.10) then
426 xx=vbld(i+nres)-vbldsc0(1,itype(i,1))
427 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
428 i,(dC(j,i+nres),j=1,3),xx
432 ! Second correction (rotational lengthening)
438 blen2 = scalar(dc(1,i),dc(1,i))
439 ppvec(ind)=2*vbl**2-blen2
440 diffbond=dabs(vbl-dsqrt(blen2))
441 if (diffbond.gt.diffmax) diffmax=diffbond
442 if (ppvec(ind).gt.0.0d0) then
443 ppvec(ind)=dsqrt(ppvec(ind))
448 write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind)
452 if (itype(i,1).ne.10) then
454 blen2 = scalar(dc(1,i+nres),dc(1,i+nres))
455 ppvec(ind)=2*vbldsc0(1,itype(i,1))**2-blen2
456 diffbond=dabs(vbldsc0(1,itype(i,1))-dsqrt(blen2))
457 if (diffbond.gt.diffmax) diffmax=diffbond
458 if (ppvec(ind).gt.0.0d0) then
459 ppvec(ind)=dsqrt(ppvec(ind))
464 write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind)
468 if (lprn) write (iout,*) "iter",iter," diffmax",diffmax
469 if (diffmax.lt.difftol) goto 10
473 Td(i)=Td(i)+ppvec(j)*Tmat(i,j)
479 zapas(i)=zapas(i)+Pmat(i,j)*dc_work(j)
490 dc_work(ind+j)=zapas(ind+j)
495 if (itype(i,1).ne.10) then
497 dc(j,i+nres)=zapas(ind+j)
498 dc_work(ind+j)=zapas(ind+j)
503 ! Building the chain from the newly calculated coordinates
506 if (large.and. mod(itime,ntwe).eq.0) then
507 write (iout,*) "Cartesian and internal coordinates: step 1"
510 write (iout,'(a)') "Potential forces"
512 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(-gcart(j,i),j=1,3),&
515 write (iout,'(a)') "Stochastic forces"
517 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(stochforc(j,i),j=1,3),&
518 (stochforc(j,i+nres),j=1,3)
520 write (iout,'(a)') "Velocities"
522 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
523 (d_t(j,i+nres),j=1,3)
528 write (iout,*) "After correction for rotational lengthening"
529 write (iout,*) "New coordinates",&
530 " and differences between actual and standard bond lengths"
535 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
539 if (itype(i,1).ne.10) then
541 xx=vbld(i+nres)-vbldsc0(1,itype(i,1))
542 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
543 i,(dC(j,i+nres),j=1,3),xx
548 ! write (iout,*) "Too many attempts at correcting the bonds"
556 ! Calculate energy and forces
558 call etotal(potEcomp)
559 potE=potEcomp(0)-potEcomp(51)
563 ! Calculate the kinetic and total energy and the kinetic temperature
566 t_enegrad=t_enegrad+MPI_Wtime()-tt0
568 t_enegrad=t_enegrad+tcpu()-tt0
571 kinetic_T=2.0d0/(dimen*Rb)*EK
573 end subroutine brown_step
574 !-----------------------------------------------------------------------------
576 !-----------------------------------------------------------------------------
577 subroutine gauss(RO,AP,MT,M,N,*)
579 ! CALCULATES (RO**(-1))*AP BY GAUSS ELIMINATION
580 ! RO IS A SQUARE MATRIX
581 ! THE CALCULATED PRODUCT IS STORED IN AP
582 ! ABNORMAL EXIT IF RO IS SINGULAR
584 integer :: MT, M, N, M1,I,J,IM,&
586 real(kind=8) :: RO(MT,M),AP(MT,N),X,RM,PR,Y
592 if(dabs(X).le.1.0D-13) return 1
604 if(DABS(RO(J,I)).LE.RM) goto 2
618 if(dabs(X).le.1.0E-13) return 1
627 8 AP(J,K)=AP(J,K)-Y*AP(I,K)
629 9 RO(J,K)=RO(J,K)-Y*RO(I,K)
633 if(dabs(X).le.1.0E-13) return 1
643 15 X=X-AP(K,J)*RO(MI,K)
648 !-----------------------------------------------------------------------------
651 subroutine kinetic(KE_total)
652 !c----------------------------------------------------------------
653 !c This subroutine calculates the total kinetic energy of the chain
654 !c-----------------------------------------------------------------
655 !c 3/5/2020 AL Corrected for multichain systems, no fake peptide groups
656 !c inside, implemented with five-diagonal inertia matrix
659 real(kind=8):: KE_total,KEt_p,KEt_sc,KEr_p,KEr_sc,mag1,mag2
660 integer i,j,k,iti,mnum
661 real(kind=8),dimension(3) :: incr,v
667 !c write (iout,*) "ISC",(isc(itype(i)),i=1,nres)
668 !c The translational part for peptide virtual bonds
672 do i=nnt,nct-1 !czy na pewno nct-1??
674 !c write (iout,*) "Kinetic trp:",i,(incr(j),j=1,3
675 !c Skip dummy peptide groups
676 if (itype(i,mnum).ne.ntyp1_molec(mnum)&
677 .and. itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
679 v(j)=incr(j)+0.5d0*d_t(j,i)
681 if (mnum.eq.5) mp(mnum)=0.0d0
682 ! if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
683 !c write (iout,*) "Kinetic trp:",i,(v(j),j=1,3)
684 vtot(i)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
685 KEt_p=KEt_p+mp(mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
688 incr(j)=incr(j)+d_t(j,i)
691 !c write(iout,*) 'KEt_p', KEt_p
692 !c The translational part for the side chain virtual bond
693 !c Only now we can initialize incr with zeros. It must be equal
694 !c to the velocities of the first Calpha.
700 iti=iabs(itype(i,mnum))
701 if (mnum.eq.5) iti=itype(i,mnum)
702 if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
709 v(j)=incr(j)+d_t(j,nres+i)
712 ! if (mnum.ne.5) then
713 ! write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3)
714 ! write (iout,*) "i",i," msc",msc(iti,mnum)," v",(v(j),j=1,3)
715 KEt_sc=KEt_sc+msc(iti,mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
716 vtot(i+nres)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
719 incr(j)=incr(j)+d_t(j,i)
723 ! write(iout,*) 'KEt_sc', KEt_sc
724 ! The part due to stretching and rotation of the peptide groups
727 if (itype(i,mnum).ne.ntyp1_molec(mnum)&
728 .and.itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
729 if (mnum.eq.5) Ip(mnum)=Isc(itype(i,mnum),mnum)
730 ! write (iout,*) "i",i
731 ! write (iout,*) "i",i," mag1",mag1," mag2",mag2
735 !c write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3)
736 KEr_p=KEr_p+Ip(mnum)*(incr(1)*incr(1)+incr(2)*incr(2) &
741 !c write(iout,*) 'KEr_p', KEr_p
742 !c The rotational part of the side chain virtual bond
745 iti=iabs(itype(i,mnum))
746 if (itype(i,1).ne.10.and.itype(i,mnum).ne.ntyp1_molec(mnum)&
749 incr(j)=d_t(j,nres+i)
751 ! write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3)
752 KEr_sc=KEr_sc+Isc(iti,mnum)*(incr(1)*incr(1)+incr(2)*incr(2)+&
756 !c The total kinetic energy
758 ! write(iout,*) ' KEt_p',KEt_p,' KEt_sc',KEt_sc,' KEr_p',KEr_p, &
760 KE_total=0.5d0*(KEt_p+KEt_sc+0.25d0*KEr_p+KEr_sc)
761 !c write (iout,*) "KE_total",KE_total
763 end subroutine kinetic
766 !-----------------------------------------------------------------------------
767 subroutine kinetic(KE_total)
768 !----------------------------------------------------------------
769 ! This subroutine calculates the total kinetic energy of the chain
770 !-----------------------------------------------------------------
772 ! implicit real*8 (a-h,o-z)
773 ! include 'DIMENSIONS'
774 ! include 'COMMON.VAR'
775 ! include 'COMMON.CHAIN'
776 ! include 'COMMON.DERIV'
777 ! include 'COMMON.GEO'
778 ! include 'COMMON.LOCAL'
779 ! include 'COMMON.INTERACT'
780 ! include 'COMMON.MD'
781 ! include 'COMMON.IOUNITS'
782 real(kind=8) :: KE_total,mscab
784 integer :: i,j,k,iti,mnum,term
785 real(kind=8) :: KEt_p,KEt_sc,KEr_p,KEr_sc,incr(3),&
788 write (iout,*) "Velocities, kietic"
790 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
791 (d_t(j,i+nres),j=1,3)
796 ! write (iout,*) "ISC",(isc(itype(i,1)),i=1,nres)
797 ! The translational part for peptide virtual bonds
802 ! if (molnum(nct).gt.3) term=nct
805 if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
806 ! write (iout,*) "Kinetic trp:",i,(incr(j),j=1,3),mp(mnum)
809 v(j)=incr(j)+0.5d0*d_t(j,i)
813 v(j)=incr(j)+0.5d0*d_t(j,i)
816 vtot(i)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
817 KEt_p=KEt_p+mp(mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
819 incr(j)=incr(j)+d_t(j,i)
822 ! write(iout,*) 'KEt_p', KEt_p
823 ! The translational part for the side chain virtual bond
824 ! Only now we can initialize incr with zeros. It must be equal
825 ! to the velocities of the first Calpha.
831 iti=iabs(itype(i,mnum))
832 ! if (mnum.ge.4) then
837 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
838 if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
839 .or.(mnum.ge.4)) then
845 v(j)=incr(j)+d_t(j,nres+i)
848 ! write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3)
849 ! write (iout,*) "i",i," msc",msc(iti,mnum)," v",(v(j),j=1,3)
850 KEt_sc=KEt_sc+mscab*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
851 vtot(i+nres)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
853 incr(j)=incr(j)+d_t(j,i)
857 ! write(iout,*) 'KEt_sc', KEt_sc
858 ! The part due to stretching and rotation of the peptide groups
862 ! write (iout,*) "i",i
863 ! write (iout,*) "i",i," mag1",mag1," mag2",mag2
867 ! write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3)
868 KEr_p=KEr_p+Ip(mnum)*(incr(1)*incr(1)+incr(2)*incr(2) &
872 ! write(iout,*) 'KEr_p', KEr_p
873 ! The rotational part of the side chain virtual bond
877 iti=iabs(itype(i,mnum))
878 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
879 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
880 .and.(mnum.lt.4)) then
882 incr(j)=d_t(j,nres+i)
884 ! write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3)
885 KEr_sc=KEr_sc+Isc(iti,mnum)*(incr(1)*incr(1)+incr(2)*incr(2)+ &
889 ! The total kinetic energy
891 ! write(iout,*) 'KEr_sc', KEr_sc
892 KE_total=0.5d0*(KEt_p+KEt_sc+0.25d0*KEr_p+KEr_sc)
893 ! write (iout,*) "KE_total",KE_total
895 end subroutine kinetic
896 !-----------------------------------------------------------------------------
898 subroutine kinetic_CASC(KE_total)
899 !c----------------------------------------------------------------
900 !c Compute the kinetic energy of the system using the Calpha-SC
902 !c-----------------------------------------------------------------
904 double precision KE_total
906 integer i,j,k,iti,ichain,innt,inct,mnum
907 double precision KEt_p,KEt_sc,KEr_p,KEr_sc,incr(3),&
914 !c write (iout,*) "ISC",(isc(itype(i)),i=1,nres)
915 !c The translational part for peptide virtual bonds
918 innt=chain_border(1,ichain)
919 inct=chain_border(2,ichain)
920 !c write (iout,*) "Kinetic_CASC chain",ichain," innt",innt,
925 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
926 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
927 !c write (iout,*) i,(d_t(j,i),j=1,3),(d_t(j,i+1),j=1,3)
929 v(j)=0.5d0*(d_t(j,i)+d_t(j,i+1))
931 !c write (iout,*) "Kinetic trp i",i," v",(v(j),j=1,3)
932 KEt_p=KEt_p+mp(mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
934 !c write(iout,*) 'KEt_p', KEt_p
935 !c The translational part for the side chain virtual bond
936 !c Only now we can initialize incr with zeros. It must be equal
937 !c to the velocities of the first Calpha.
940 iti=iabs(itype(i,mnum))
941 if (itype(i,1).eq.10.or.mnum.ge.3.or. itype(i,mnum).eq.ntyp1_molec(mnum)) then
942 !c write (iout,*) i,iti,(d_t(j,i),j=1,3)
947 !c write (iout,*) i,iti,(d_t(j,nres+i),j=1,3)
952 !c write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3)
953 !c write (iout,*) "i",i," msc",msc(iti)," v",(v(j),j=1,3)
954 KEt_sc=KEt_sc+msc(iti,mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
957 !c write(iout,*) 'KEt_sc', KEt_sc
958 !c The part due to stretching and rotation of the peptide groups
962 incr(j)=d_t(j,i+1)-d_t(j,i)
964 if (mnum.eq.5) Ip(mnum)=Isc(itype(i,mnum),mnum)
965 !c write (iout,*) i,(incr(j),j=1,3)
966 !c write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3)
967 KEr_p=KEr_p+Ip(mnum)*(incr(1)*incr(1)+incr(2)*incr(2)&
971 !c write(iout,*) 'KEr_p', KEr_p
972 !c The rotational part of the side chain virtual bond
975 iti=iabs(itype(i,mnum))
976 ! if (iti.ne.10.and.mnum.lt.3) then
977 if (itype(i,1).ne.10.and.mnum.lt.3.and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
979 incr(j)=d_t(j,nres+i)-d_t(j,i)
981 !c write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3)
982 KEr_sc=KEr_sc+Isc(iti,mnum)*(incr(1)*incr(1)+incr(2)*incr(2)+&
988 !c The total kinetic energy
990 !c write(iout,*) ' KEt_p',KEt_p,' KEt_sc',KEt_sc,' KEr_p',KEr_p,
991 !c & ' KEr_sc', KEr_sc
992 KE_total=0.5d0*(KEt_p+KEt_sc+0.25d0*KEr_p+KEr_sc)
993 !c write (iout,*) "KE_total",KE_tota
995 write (iout,*) "Need to compile with -DFIVEDIAG to use this sub!"
999 end subroutine kinetic_CASC
1002 !-----------------------------------------------------------------------------
1004 !------------------------------------------------
1005 ! The driver for molecular dynamics subroutines
1006 !------------------------------------------------
1009 use control, only:tcpu,ovrtim
1010 ! use io_comm, only:ilen
1012 use compare, only:secondary2,hairpin
1013 use io, only:cartout,statout
1014 ! implicit real*8 (a-h,o-z)
1015 ! include 'DIMENSIONS'
1018 integer :: IERROR,ERRCODE
1020 ! include 'COMMON.SETUP'
1021 ! include 'COMMON.CONTROL'
1022 ! include 'COMMON.VAR'
1023 ! include 'COMMON.MD'
1025 ! include 'COMMON.LANGEVIN'
1027 ! include 'COMMON.LANGEVIN.lang0'
1029 ! include 'COMMON.CHAIN'
1030 ! include 'COMMON.DERIV'
1031 ! include 'COMMON.GEO'
1032 ! include 'COMMON.LOCAL'
1033 ! include 'COMMON.INTERACT'
1034 ! include 'COMMON.IOUNITS'
1035 ! include 'COMMON.NAMES'
1036 ! include 'COMMON.TIME1'
1037 ! include 'COMMON.HAIRPIN'
1038 real(kind=8),dimension(3) :: L,vcm,boxx
1040 real(kind=8),dimension(6*nres) :: v_work,v_transf !(maxres6) maxres6=6*maxres
1042 integer :: rstcount !ilen,
1044 character(len=50) :: tytul
1045 !el common /gucio/ cm
1046 integer :: i,j,nharp
1047 integer,dimension(4,nres) :: iharp !(4,nres/3)(4,maxres/3)
1050 real(kind=8) :: tt0,scalfac
1051 integer :: nres2,itime
1060 print *,"MY tmpdir",tmpdir,ilen(tmpdir)
1061 if (ilen(tmpdir).gt.0) &
1062 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_" &
1063 //liczba(:ilen(liczba))//'.rst')
1065 if (ilen(tmpdir).gt.0) &
1066 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"//'.rst')
1073 write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
1079 print *,"just befor setup matix",nres
1080 ! Determine the inverse of the inertia matrix.
1081 call setup_MD_matrices
1083 print *,"AFTER SETUP MATRICES"
1085 print *,"AFTER INIT MD"
1088 t_MDsetup = MPI_Wtime()-tt0
1090 t_MDsetup = tcpu()-tt0
1093 ! Entering the MD loop
1099 if (lang.eq.2 .or. lang.eq.3) then
1103 call sd_verlet_p_setup
1105 call sd_verlet_ciccotti_setup
1109 pfric0_mat(i,j,0)=pfric_mat(i,j)
1110 afric0_mat(i,j,0)=afric_mat(i,j)
1111 vfric0_mat(i,j,0)=vfric_mat(i,j)
1112 prand0_mat(i,j,0)=prand_mat(i,j)
1113 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
1114 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
1117 flag_stoch(0)=.true.
1118 do i=1,maxflag_stoch
1119 flag_stoch(i)=.false.
1123 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
1125 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1129 else if (lang.eq.1 .or. lang.eq.4) then
1130 print *,"before setup_fricmat"
1132 print *,"after setup_fricmat"
1135 t_langsetup=MPI_Wtime()-tt0
1138 t_langsetup=tcpu()-tt0
1141 do itime=1,n_timestep
1142 if (large) print *,itime,ntwe
1144 if (large.and. mod(itime,ntwe).eq.0) &
1145 write (iout,*) "itime",itime
1147 if (lang.gt.0 .and. surfarea .and. &
1148 mod(itime,reset_fricmat).eq.0) then
1149 if (lang.eq.2 .or. lang.eq.3) then
1153 call sd_verlet_p_setup
1155 call sd_verlet_ciccotti_setup
1159 pfric0_mat(i,j,0)=pfric_mat(i,j)
1160 afric0_mat(i,j,0)=afric_mat(i,j)
1161 vfric0_mat(i,j,0)=vfric_mat(i,j)
1162 prand0_mat(i,j,0)=prand_mat(i,j)
1163 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
1164 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
1167 flag_stoch(0)=.true.
1168 do i=1,maxflag_stoch
1169 flag_stoch(i)=.false.
1172 else if (lang.eq.1 .or. lang.eq.4) then
1173 print *,"before setup_fricmat"
1175 print *,"after setup_fricmat"
1177 write (iout,'(a,i10)') &
1178 "Friction matrix reset based on surface area, itime",itime
1180 if (reset_vel .and. tbf .and. lang.eq.0 &
1181 .and. mod(itime,count_reset_vel).eq.0) then
1184 if (molnum(i).eq.5) then
1185 call to_box(c(1,i),c(2,i),c(3,i))
1187 if (c(j,i).le.0) c(j,i)=c(j,i)+boxx(j)
1190 ! write(iout,*) "COORD",c(1,i),c(2,i),c(3,i)
1194 write(iout,'(a,f20.2)') &
1195 "Velocities reset to random values, time",totT
1198 d_t_old(j,i)=d_t(j,i)
1202 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
1206 d_t(j,0)=d_t(j,0)-vcm(j)
1209 kinetic_T=2.0d0/(dimen3*Rb)*EK
1210 scalfac=dsqrt(T_bath/kinetic_T)
1211 write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
1214 d_t_old(j,i)=scalfac*d_t(j,i)
1220 ! Time-reversible RESPA algorithm
1221 ! (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
1222 call RESPA_step(itime)
1224 ! Variable time step algorithm.
1225 if (large) print *,"before verlet_step"
1226 call velverlet_step(itime)
1227 if (large) print *,"after verlet_step"
1231 call brown_step(itime)
1233 print *,"Brown dynamics not here!"
1235 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1242 if (mod(itime,ntwe).eq.0) then
1246 ! call check_ecartint
1256 v_work(ind)=d_t(j,i)
1261 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum).and.mnum.lt.4) then
1264 v_work(ind)=d_t(j,i+nres)
1269 write (66,'(80f10.5)') &
1270 ((d_t(j,i),j=1,3),i=0,nres-1),((d_t(j,i+nres),j=1,3),i=1,nres)
1274 v_transf(i)=v_transf(i)+gvec(j,i)*v_work(j)
1276 v_transf(i)= v_transf(i)*dsqrt(geigen(i))
1278 write (67,'(80f10.5)') (v_transf(i),i=1,ind)
1281 if (mod(itime,ntwx).eq.0) then
1283 call enerprint(potEcomp)
1285 write (tytul,'("time",f8.2)') totT
1287 call hairpin(.true.,nharp,iharp)
1288 call secondary2(.true.)
1289 call pdbout(potE,tytul,ipdb)
1290 ! call enerprint(potEcomp)
1295 write(iout,*) "starting fodstep"
1296 call fodstep(nfodstep)
1297 write(iout,*) "after fodstep"
1300 call hairpin(.true.,nharp,iharp)
1301 call secondary2(.true.)
1302 call pdbout(potE,tytul,ipdb)
1309 if (rstcount.eq.1000.or.itime.eq.n_timestep) then
1310 open(irest2,file=rest2name,status='unknown')
1311 write(irest2,*) totT,EK,potE,totE,t_bath
1313 ! AL 4/17/17: Now writing d_t(0,:) too
1315 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
1317 ! AL 4/17/17: Now writing d_c(0,:) too
1319 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
1327 t_MD=MPI_Wtime()-tt0
1331 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))') &
1333 'MD calculations setup:',t_MDsetup,&
1334 'Energy & gradient evaluation:',t_enegrad,&
1335 'Stochastic MD setup:',t_langsetup,&
1336 'Stochastic MD step setup:',t_sdsetup,&
1338 write (iout,'(/28(1h=),a25,27(1h=))') &
1339 ' End of MD calculation '
1341 write (iout,*) "time for etotal",t_etotal," elong",t_elong,&
1343 write (iout,*) "time_fric",time_fric," time_stoch",time_stoch,&
1344 " time_fricmatmult",time_fricmatmult," time_fsample ",&
1349 !-----------------------------------------------------------------------------
1350 subroutine velverlet_step(itime)
1351 !-------------------------------------------------------------------------------
1352 ! Perform a single velocity Verlet step; the time step can be rescaled if
1353 ! increments in accelerations exceed the threshold
1354 !-------------------------------------------------------------------------------
1355 ! implicit real*8 (a-h,o-z)
1356 ! include 'DIMENSIONS'
1358 use control, only:tcpu
1360 use minimm, only:minim_dc
1363 integer :: ierror,ierrcode
1364 real(kind=8) :: errcode
1366 ! include 'COMMON.SETUP'
1367 ! include 'COMMON.VAR'
1368 ! include 'COMMON.MD'
1370 ! include 'COMMON.LANGEVIN'
1372 ! include 'COMMON.LANGEVIN.lang0'
1374 ! include 'COMMON.CHAIN'
1375 ! include 'COMMON.DERIV'
1376 ! include 'COMMON.GEO'
1377 ! include 'COMMON.LOCAL'
1378 ! include 'COMMON.INTERACT'
1379 ! include 'COMMON.IOUNITS'
1380 ! include 'COMMON.NAMES'
1381 ! include 'COMMON.TIME1'
1382 ! include 'COMMON.MUCA'
1383 real(kind=8),dimension(3) :: vcm,incr
1384 real(kind=8),dimension(3) :: L
1385 integer :: count,rstcount !ilen,
1387 character(len=50) :: tytul
1388 integer :: maxcount_scale = 30
1389 !el common /gucio/ cm
1390 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1391 !el common /stochcalc/ stochforcvec
1392 integer :: icount_scale,itime_scal,i,j,ifac_time,iretcode,itime
1394 real(kind=8) :: epdrift,tt0,fac_time
1396 if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres)) !(MAXRES6) maxres6=6*maxres
1402 if (large) print *,"after sddir_precalc"
1403 else if (lang.eq.2 .or. lang.eq.3) then
1405 call stochastic_force(stochforcvec)
1408 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
1410 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1417 icount_scale=icount_scale+1
1418 ! write(iout,*) "icount_scale",icount_scale,scalel
1419 if (icount_scale.gt.maxcount_scale) then
1421 "ERROR: too many attempts at scaling down the time step. ",&
1422 "amax=",amax,"epdrift=",epdrift,&
1423 "damax=",damax,"edriftmax=",edriftmax,&
1427 call MPI_Abort(MPI_COMM_WORLD,IERROR,IERRCODE)
1431 ! First step of the velocity Verlet algorithm
1436 else if (lang.eq.3) then
1438 call sd_verlet1_ciccotti
1440 else if (lang.eq.1) then
1445 ! Build the chain from the newly calculated coordinates
1446 call chainbuild_cart
1447 if (rattle) call rattle1
1449 if (large) then !.and. mod(itime,ntwe).eq.0) then
1450 write (iout,*) "Cartesian and internal coordinates: step 1"
1455 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),&
1456 (dc(j,i+nres),j=1,3)
1458 write (iout,*) "Accelerations"
1460 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1461 (d_a(j,i+nres),j=1,3)
1463 write (iout,*) "Velocities, step 1"
1465 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1466 (d_t(j,i+nres),j=1,3)
1475 ! Calculate energy and forces
1477 call etotal(potEcomp)
1478 ! AL 4/17/17: Reduce the steps if NaNs occurred.
1479 if (potEcomp(0).gt.0.99e18 .or. isnan(potEcomp(0)).gt.0) then
1480 call enerprint(potEcomp)
1482 if (icount_scale.gt.15) then
1483 write (iout,*) "Tu jest problem",potEcomp(0),d_time
1484 ! call gen_rand_conf(1,*335)
1485 ! call minim_dc(potEcomp(0),iretcode,100)
1488 ! call etotal(potEcomp)
1489 ! write(iout,*) "needed to repara,",potEcomp
1492 ! 335 write(iout,*) "Failed genrand"
1496 if (large.and. mod(itime,ntwe).eq.0) &
1497 call enerprint(potEcomp)
1500 t_etotal=t_etotal+MPI_Wtime()-tt0
1502 t_etotal=t_etotal+tcpu()-tt0
1505 potE=potEcomp(0)-potEcomp(51)
1507 ! Get the new accelerations
1510 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1512 t_enegrad=t_enegrad+tcpu()-tt0
1514 ! Determine maximum acceleration and scale down the timestep if needed
1516 amax=amax/(itime_scal+1)**2
1517 call predict_edrift(epdrift)
1518 ! write(iout,*) "amax=",amax,damax,epdrift,edriftmax,amax/(itime_scal+1)
1520 ! write (iout,*) "before enter if",scalel,icount_scale
1521 if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then
1522 ! write(iout,*) "I enter if"
1523 ! Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
1525 ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax)) &
1527 itime_scal=itime_scal+ifac_time
1528 ! fac_time=dmin1(damax/amax,0.5d0)
1529 fac_time=0.5d0**ifac_time
1530 d_time=d_time*fac_time
1531 if (lang.eq.2 .or. lang.eq.3) then
1533 ! write (iout,*) "Calling sd_verlet_setup: 1"
1534 ! Rescale the stochastic forces and recalculate or restore
1535 ! the matrices of tinker integrator
1536 if (itime_scal.gt.maxflag_stoch) then
1537 if (large) write (iout,'(a,i5,a)') &
1538 "Calculate matrices for stochastic step;",&
1539 " itime_scal ",itime_scal
1541 call sd_verlet_p_setup
1543 call sd_verlet_ciccotti_setup
1545 write (iout,'(2a,i3,a,i3,1h.)') &
1546 "Warning: cannot store matrices for stochastic",&
1547 " integration because the index",itime_scal,&
1548 " is greater than",maxflag_stoch
1549 write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",&
1550 " integration Langevin algorithm for better efficiency."
1551 else if (flag_stoch(itime_scal)) then
1552 if (large) write (iout,'(a,i5,a,l1)') &
1553 "Restore matrices for stochastic step; itime_scal ",&
1554 itime_scal," flag ",flag_stoch(itime_scal)
1557 pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
1558 afric_mat(i,j)=afric0_mat(i,j,itime_scal)
1559 vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
1560 prand_mat(i,j)=prand0_mat(i,j,itime_scal)
1561 vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
1562 vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
1566 if (large) write (iout,'(2a,i5,a,l1)') &
1567 "Calculate & store matrices for stochastic step;",&
1568 " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
1570 call sd_verlet_p_setup
1572 call sd_verlet_ciccotti_setup
1574 flag_stoch(ifac_time)=.true.
1577 pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
1578 afric0_mat(i,j,itime_scal)=afric_mat(i,j)
1579 vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
1580 prand0_mat(i,j,itime_scal)=prand_mat(i,j)
1581 vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
1582 vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
1586 fac_time=1.0d0/dsqrt(fac_time)
1588 stochforcvec(i)=fac_time*stochforcvec(i)
1591 else if (lang.eq.1) then
1592 ! Rescale the accelerations due to stochastic forces
1593 fac_time=1.0d0/dsqrt(fac_time)
1595 d_as_work(i)=d_as_work(i)*fac_time
1598 if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)') &
1599 "itime",itime," Timestep scaled down to ",&
1600 d_time," ifac_time",ifac_time," itime_scal",itime_scal
1602 ! Second step of the velocity Verlet algorithm
1607 else if (lang.eq.3) then
1609 call sd_verlet2_ciccotti
1611 else if (lang.eq.1) then
1616 if (rattle) call rattle2
1619 if (d_time.ne.d_time0) then
1622 if (lang.eq.2 .or. lang.eq.3) then
1623 if (large) write (iout,'(a)') &
1624 "Restore original matrices for stochastic step"
1625 ! write (iout,*) "Calling sd_verlet_setup: 2"
1626 ! Restore the matrices of tinker integrator if the time step has been restored
1629 pfric_mat(i,j)=pfric0_mat(i,j,0)
1630 afric_mat(i,j)=afric0_mat(i,j,0)
1631 vfric_mat(i,j)=vfric0_mat(i,j,0)
1632 prand_mat(i,j)=prand0_mat(i,j,0)
1633 vrand_mat1(i,j)=vrand0_mat1(i,j,0)
1634 vrand_mat2(i,j)=vrand0_mat2(i,j,0)
1642 ! Calculate the kinetic and the total energy and the kinetic temperature
1646 ! call kinetic1(EK1)
1647 ! write (iout,*) "step",itime," EK",EK," EK1",EK1
1649 ! Couple the system to Berendsen bath if needed
1650 if (tbf .and. lang.eq.0) then
1653 kinetic_T=2.0d0/(dimen3*Rb)*EK
1654 ! Backup the coordinates, velocities, and accelerations
1658 d_t_old(j,i)=d_t(j,i)
1659 d_a_old(j,i)=d_a(j,i)
1663 if (mod(itime,ntwe).eq.0 .and. large) then
1664 write (iout,*) "Velocities, step 2"
1666 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1667 (d_t(j,i+nres),j=1,3)
1672 end subroutine velverlet_step
1673 !-----------------------------------------------------------------------------
1674 subroutine RESPA_step(itime)
1675 !-------------------------------------------------------------------------------
1676 ! Perform a single RESPA step.
1677 !-------------------------------------------------------------------------------
1678 ! implicit real*8 (a-h,o-z)
1679 ! include 'DIMENSIONS'
1683 use control, only:tcpu
1685 ! use io_conf, only:cartprint
1688 integer :: IERROR,ERRCODE
1690 ! include 'COMMON.SETUP'
1691 ! include 'COMMON.CONTROL'
1692 ! include 'COMMON.VAR'
1693 ! include 'COMMON.MD'
1695 ! include 'COMMON.LANGEVIN'
1697 ! include 'COMMON.LANGEVIN.lang0'
1699 ! include 'COMMON.CHAIN'
1700 ! include 'COMMON.DERIV'
1701 ! include 'COMMON.GEO'
1702 ! include 'COMMON.LOCAL'
1703 ! include 'COMMON.INTERACT'
1704 ! include 'COMMON.IOUNITS'
1705 ! include 'COMMON.NAMES'
1706 ! include 'COMMON.TIME1'
1707 real(kind=8),dimension(0:n_ene) :: energia_short,energia_long
1708 real(kind=8),dimension(3) :: L,vcm,incr
1709 real(kind=8),dimension(3,0:2*nres) :: dc_old0,d_t_old0,d_a_old0 !(3,0:maxres2) maxres2=2*maxres
1710 logical :: PRINT_AMTS_MSG = .false.
1711 integer :: count,rstcount !ilen,
1713 character(len=50) :: tytul
1714 integer :: maxcount_scale = 10
1715 !el common /gucio/ cm
1716 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1717 !el common /stochcalc/ stochforcvec
1718 integer :: itt,i,j,itsplit,itime
1720 !el common /cipiszcze/ itt
1722 real(kind=8) :: epdrift,tt0,epdriftmax
1725 if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres)) !(MAXRES6) maxres6=6*maxres
1729 if (large.and. mod(itime,ntwe).eq.0) then
1730 write (iout,*) "***************** RESPA itime",itime
1731 write (iout,*) "Cartesian and internal coordinates: step 0"
1733 call pdbout(0.0d0,"cipiszcze",iout)
1735 write (iout,*) "Accelerations from long-range forces"
1737 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1738 (d_a(j,i+nres),j=1,3)
1740 write (iout,*) "Velocities, step 0"
1742 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1743 (d_t(j,i+nres),j=1,3)
1748 ! Perform the initial RESPA step (increment velocities)
1749 ! write (iout,*) "*********************** RESPA ini"
1752 if (mod(itime,ntwe).eq.0 .and. large) then
1753 write (iout,*) "Velocities, end"
1755 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1756 (d_t(j,i+nres),j=1,3)
1760 ! Compute the short-range forces
1766 ! 7/2/2009 commented out
1768 ! call etotal_short(energia_short)
1771 ! 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
1774 d_a(j,i)=d_a_short(j,i)
1778 if (large.and. mod(itime,ntwe).eq.0) then
1779 write (iout,*) "energia_short",energia_short(0)
1780 write (iout,*) "Accelerations from short-range forces"
1782 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1783 (d_a(j,i+nres),j=1,3)
1788 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1790 t_enegrad=t_enegrad+tcpu()-tt0
1795 d_t_old(j,i)=d_t(j,i)
1796 d_a_old(j,i)=d_a(j,i)
1799 ! 6/30/08 A-MTS: attempt at increasing the split number
1802 dc_old0(j,i)=dc_old(j,i)
1803 d_t_old0(j,i)=d_t_old(j,i)
1804 d_a_old0(j,i)=d_a_old(j,i)
1807 if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
1808 if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
1815 ! write (iout,*) "itime",itime," ntime_split",ntime_split
1816 ! Split the time step
1817 d_time=d_time0/ntime_split
1818 ! Perform the short-range RESPA steps (velocity Verlet increments of
1819 ! positions and velocities using short-range forces)
1820 ! write (iout,*) "*********************** RESPA split"
1821 do itsplit=1,ntime_split
1824 else if (lang.eq.2 .or. lang.eq.3) then
1826 call stochastic_force(stochforcvec)
1829 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
1831 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1836 ! First step of the velocity Verlet algorithm
1841 else if (lang.eq.3) then
1843 call sd_verlet1_ciccotti
1845 else if (lang.eq.1) then
1850 ! Build the chain from the newly calculated coordinates
1851 call chainbuild_cart
1852 if (rattle) call rattle1
1854 if (large.and. mod(itime,ntwe).eq.0) then
1855 write (iout,*) "***** ITSPLIT",itsplit
1856 write (iout,*) "Cartesian and internal coordinates: step 1"
1857 call pdbout(0.0d0,"cipiszcze",iout)
1860 write (iout,*) "Velocities, step 1"
1862 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1863 (d_t(j,i+nres),j=1,3)
1872 ! Calculate energy and forces
1874 call etotal_short(energia_short)
1875 ! AL 4/17/17: Exit itime_split loop when energy goes infinite
1876 if (energia_short(0).gt.0.99e20 .or. isnan(energia_short(0)) ) then
1877 if (PRINT_AMTS_MSG) &
1878 write (iout,*) "Infinities/NaNs in energia_short",energia_short(0),"; increasing ntime_split to",ntime_split
1879 ntime_split=ntime_split*2
1880 if (ntime_split.gt.maxtime_split) then
1883 "Cannot rescue the run; aborting job. Retry with a smaller time step"
1885 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1888 "Cannot rescue the run; terminating. Retry with a smaller time step"
1894 if (large.and. mod(itime,ntwe).eq.0) &
1895 call enerprint(energia_short)
1898 t_eshort=t_eshort+MPI_Wtime()-tt0
1900 t_eshort=t_eshort+tcpu()-tt0
1904 ! Get the new accelerations
1906 ! 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
1909 d_a_short(j,i)=d_a(j,i)
1913 if (large.and. mod(itime,ntwe).eq.0) then
1914 write (iout,*)"energia_short",energia_short(0)
1915 write (iout,*) "Accelerations from short-range forces"
1917 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1918 (d_a(j,i+nres),j=1,3)
1923 ! Determine maximum acceleration and scale down the timestep if needed
1925 amax=amax/ntime_split**2
1926 call predict_edrift(epdrift)
1927 if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0) &
1928 write (iout,*) "amax",amax," damax",damax,&
1929 " epdrift",epdrift," epdriftmax",epdriftmax
1930 ! Exit loop and try with increased split number if the change of
1931 ! acceleration is too big
1932 if (amax.gt.damax .or. epdrift.gt.edriftmax) then
1933 if (ntime_split.lt.maxtime_split) then
1935 ntime_split=ntime_split*2
1936 ! AL 4/17/17: We should exit the itime_split loop when acceleration change is too big
1940 dc_old(j,i)=dc_old0(j,i)
1941 d_t_old(j,i)=d_t_old0(j,i)
1942 d_a_old(j,i)=d_a_old0(j,i)
1945 if (PRINT_AMTS_MSG) then
1946 write (iout,*) "acceleration/energy drift too large",amax,&
1947 epdrift," split increased to ",ntime_split," itime",itime,&
1953 "Uh-hu. Bumpy landscape. Maximum splitting number",&
1955 " already reached!!! Trying to carry on!"
1959 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1961 t_enegrad=t_enegrad+tcpu()-tt0
1963 ! Second step of the velocity Verlet algorithm
1968 else if (lang.eq.3) then
1970 call sd_verlet2_ciccotti
1972 else if (lang.eq.1) then
1977 if (rattle) call rattle2
1978 ! Backup the coordinates, velocities, and accelerations
1982 d_t_old(j,i)=d_t(j,i)
1983 d_a_old(j,i)=d_a(j,i)
1990 ! Restore the time step
1992 ! Compute long-range forces
1999 call etotal_long(energia_long)
2000 if (energia_long(0).gt.0.99e20 .or. isnan(energia_long(0))) then
2003 "Infinitied/NaNs in energia_long, Aborting MPI job."
2005 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
2007 write (iout,*) "Infinitied/NaNs in energia_long, terminating."
2011 if (large.and. mod(itime,ntwe).eq.0) &
2012 call enerprint(energia_long)
2015 t_elong=t_elong+MPI_Wtime()-tt0
2017 t_elong=t_elong+tcpu()-tt0
2020 potE=potEcomp(0)-potEcomp(51)
2024 t_enegrad=t_enegrad+MPI_Wtime()-tt0
2026 t_enegrad=t_enegrad+tcpu()-tt0
2028 ! Compute accelerations from long-range forces
2030 if (large.and. mod(itime,ntwe).eq.0) then
2031 write (iout,*) "energia_long",energia_long(0)
2032 write (iout,*) "Cartesian and internal coordinates: step 2"
2034 call pdbout(0.0d0,"cipiszcze",iout)
2036 write (iout,*) "Accelerations from long-range forces"
2038 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
2039 (d_a(j,i+nres),j=1,3)
2041 write (iout,*) "Velocities, step 2"
2043 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2044 (d_t(j,i+nres),j=1,3)
2048 ! Compute the final RESPA step (increment velocities)
2049 ! write (iout,*) "*********************** RESPA fin"
2051 ! Compute the complete potential energy
2053 potEcomp(i)=energia_short(i)+energia_long(i)
2055 potE=potEcomp(0)-potEcomp(51)
2056 ! potE=energia_short(0)+energia_long(0)
2059 ! Calculate the kinetic and the total energy and the kinetic temperature
2062 ! Couple the system to Berendsen bath if needed
2063 if (tbf .and. lang.eq.0) then
2066 kinetic_T=2.0d0/(dimen3*Rb)*EK
2067 ! Backup the coordinates, velocities, and accelerations
2069 if (mod(itime,ntwe).eq.0 .and. large) then
2070 write (iout,*) "Velocities, end"
2072 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2073 (d_t(j,i+nres),j=1,3)
2078 end subroutine RESPA_step
2079 !-----------------------------------------------------------------------------
2080 subroutine RESPA_vel
2081 ! First and last RESPA step (incrementing velocities using long-range
2084 ! implicit real*8 (a-h,o-z)
2085 ! include 'DIMENSIONS'
2086 ! include 'COMMON.CONTROL'
2087 ! include 'COMMON.VAR'
2088 ! include 'COMMON.MD'
2089 ! include 'COMMON.CHAIN'
2090 ! include 'COMMON.DERIV'
2091 ! include 'COMMON.GEO'
2092 ! include 'COMMON.LOCAL'
2093 ! include 'COMMON.INTERACT'
2094 ! include 'COMMON.IOUNITS'
2095 ! include 'COMMON.NAMES'
2096 integer :: i,j,inres,mnum
2099 d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
2103 d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
2108 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2109 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2110 .and.(mnum.lt.4)) then
2113 d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
2118 end subroutine RESPA_vel
2119 !-----------------------------------------------------------------------------
2121 ! Applying velocity Verlet algorithm - step 1 to coordinates
2123 ! implicit real*8 (a-h,o-z)
2124 ! include 'DIMENSIONS'
2125 ! include 'COMMON.CONTROL'
2126 ! include 'COMMON.VAR'
2127 ! include 'COMMON.MD'
2128 ! include 'COMMON.CHAIN'
2129 ! include 'COMMON.DERIV'
2130 ! include 'COMMON.GEO'
2131 ! include 'COMMON.LOCAL'
2132 ! include 'COMMON.INTERACT'
2133 ! include 'COMMON.IOUNITS'
2134 ! include 'COMMON.NAMES'
2135 real(kind=8) :: adt,adt2
2136 integer :: i,j,inres,mnum
2139 write (iout,*) "VELVERLET1 START: DC"
2141 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),&
2142 (dc(j,i+nres),j=1,3)
2146 adt=d_a_old(j,0)*d_time
2148 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
2149 d_t_new(j,0)=d_t_old(j,0)+adt2
2150 d_t(j,0)=d_t_old(j,0)+adt
2154 adt=d_a_old(j,i)*d_time
2156 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
2157 d_t_new(j,i)=d_t_old(j,i)+adt2
2158 d_t(j,i)=d_t_old(j,i)+adt
2163 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2164 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2165 .and.(mnum.lt.4)) then
2168 adt=d_a_old(j,inres)*d_time
2170 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
2171 d_t_new(j,inres)=d_t_old(j,inres)+adt2
2172 d_t(j,inres)=d_t_old(j,inres)+adt
2177 write (iout,*) "VELVERLET1 END: DC"
2179 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),&
2180 (dc(j,i+nres),j=1,3)
2184 end subroutine verlet1
2185 !-----------------------------------------------------------------------------
2187 ! Step 2 of the velocity Verlet algorithm: update velocities
2189 ! implicit real*8 (a-h,o-z)
2190 ! include 'DIMENSIONS'
2191 ! include 'COMMON.CONTROL'
2192 ! include 'COMMON.VAR'
2193 ! include 'COMMON.MD'
2194 ! include 'COMMON.CHAIN'
2195 ! include 'COMMON.DERIV'
2196 ! include 'COMMON.GEO'
2197 ! include 'COMMON.LOCAL'
2198 ! include 'COMMON.INTERACT'
2199 ! include 'COMMON.IOUNITS'
2200 ! include 'COMMON.NAMES'
2201 integer :: i,j,inres,mnum
2204 d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
2208 d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
2213 ! iti=iabs(itype(i,mnum))
2214 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2215 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2216 .and.(mnum.lt.4)) then
2219 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
2224 end subroutine verlet2
2225 !-----------------------------------------------------------------------------
2226 subroutine sddir_precalc
2227 ! Applying velocity Verlet algorithm - step 1 to coordinates
2228 ! implicit real*8 (a-h,o-z)
2229 ! include 'DIMENSIONS'
2235 ! include 'COMMON.CONTROL'
2236 ! include 'COMMON.VAR'
2237 ! include 'COMMON.MD'
2239 ! include 'COMMON.LANGEVIN'
2241 ! include 'COMMON.LANGEVIN.lang0'
2243 ! include 'COMMON.CHAIN'
2244 ! include 'COMMON.DERIV'
2245 ! include 'COMMON.GEO'
2246 ! include 'COMMON.LOCAL'
2247 ! include 'COMMON.INTERACT'
2248 ! include 'COMMON.IOUNITS'
2249 ! include 'COMMON.NAMES'
2250 ! include 'COMMON.TIME1'
2251 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
2252 !el common /stochcalc/ stochforcvec
2253 real(kind=8) :: time00
2256 ! Compute friction and stochastic forces
2260 if (large) print *,"before friction_force"
2262 if (large) print *,"after friction_force"
2263 time_fric=time_fric+MPI_Wtime()-time00
2265 call stochastic_force(stochforcvec)
2266 time_stoch=time_stoch+MPI_Wtime()-time00
2269 ! Compute the acceleration due to friction forces (d_af_work) and stochastic
2270 ! forces (d_as_work)
2272 ! call ginv_mult(fric_work, d_af_work)
2273 ! call ginv_mult(stochforcvec, d_as_work)
2275 call fivediaginv_mult(dimen,fric_work, d_af_work)
2276 call fivediaginv_mult(dimen,stochforcvec, d_as_work)
2278 write(iout,*),"dimen",dimen
2280 write(iout,*) "fricwork",fric_work(i), d_af_work(i)
2281 write(iout,*) "stochforcevec", stochforcvec(i), d_as_work(i)
2285 call ginv_mult(fric_work, d_af_work)
2286 call ginv_mult(stochforcvec, d_as_work)
2290 end subroutine sddir_precalc
2291 !-----------------------------------------------------------------------------
2292 subroutine sddir_verlet1
2293 ! Applying velocity Verlet algorithm - step 1 to velocities
2296 ! implicit real*8 (a-h,o-z)
2297 ! include 'DIMENSIONS'
2298 ! include 'COMMON.CONTROL'
2299 ! include 'COMMON.VAR'
2300 ! include 'COMMON.MD'
2302 ! include 'COMMON.LANGEVIN'
2304 ! include 'COMMON.LANGEVIN.lang0'
2306 ! include 'COMMON.CHAIN'
2307 ! include 'COMMON.DERIV'
2308 ! include 'COMMON.GEO'
2309 ! include 'COMMON.LOCAL'
2310 ! include 'COMMON.INTERACT'
2311 ! include 'COMMON.IOUNITS'
2312 ! include 'COMMON.NAMES'
2313 ! Revised 3/31/05 AL: correlation between random contributions to
2314 ! position and velocity increments included.
2315 real(kind=8) :: sqrt13 = 0.57735026918962576451d0 ! 1/sqrt(3)
2316 real(kind=8) :: adt,adt2
2317 integer :: i,j,ind,inres,mnum
2319 ! Add the contribution from BOTH friction and stochastic force to the
2320 ! coordinates, but ONLY the contribution from the friction forces to velocities
2323 adt=(d_a_old(j,0)+d_af_work(j))*d_time
2324 adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
2325 ! write(iout,*) i,"adt",adt,"ads",adt2,d_a_old(j,0),d_af_work(j),d_time
2326 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
2327 d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
2328 d_t(j,0)=d_t_old(j,0)+adt
2333 adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
2334 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
2335 ! write(iout,*) i,"adt",adt,"ads",adt2,d_a_old(j,i),d_af_work(ind+j)
2336 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
2337 d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
2338 d_t(j,i)=d_t_old(j,i)+adt
2344 ! iti=iabs(itype(i,mnum))
2345 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2346 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2347 .and.(mnum.lt.4)) then
2350 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
2351 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
2352 ! write(iout,*) i,"adt",adt,"ads",adt2,d_a_old(j,inres),d_af_work(ind+j)
2353 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
2354 d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
2355 d_t(j,inres)=d_t_old(j,inres)+adt
2362 end subroutine sddir_verlet1
2363 !-----------------------------------------------------------------------------
2364 subroutine sddir_verlet2
2365 ! Calculating the adjusted velocities for accelerations
2368 ! implicit real*8 (a-h,o-z)
2369 ! include 'DIMENSIONS'
2370 ! include 'COMMON.CONTROL'
2371 ! include 'COMMON.VAR'
2372 ! include 'COMMON.MD'
2374 ! include 'COMMON.LANGEVIN'
2376 ! include 'COMMON.LANGEVIN.lang0'
2378 ! include 'COMMON.CHAIN'
2379 ! include 'COMMON.DERIV'
2380 ! include 'COMMON.GEO'
2381 ! include 'COMMON.LOCAL'
2382 ! include 'COMMON.INTERACT'
2383 ! include 'COMMON.IOUNITS'
2384 ! include 'COMMON.NAMES'
2385 real(kind=8),dimension(6*nres) :: stochforcvec,d_as_work1 !(MAXRES6) maxres6=6*maxres
2386 real(kind=8) :: cos60 = 0.5d0, sin60 = 0.86602540378443864676d0
2387 integer :: i,j,ind,inres,mnum
2388 ! Revised 3/31/05 AL: correlation between random contributions to
2389 ! position and velocity increments included.
2390 ! The correlation coefficients are calculated at low-friction limit.
2391 ! Also, friction forces are now not calculated with new velocities.
2393 ! call friction_force
2394 call stochastic_force(stochforcvec)
2396 ! Compute the acceleration due to friction forces (d_af_work) and stochastic
2397 ! forces (d_as_work)
2400 call fivediaginv_mult(6*nres,stochforcvec, d_as_work1)
2402 call ginv_mult(stochforcvec, d_as_work1)
2409 d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j)) &
2410 +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
2415 d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j)) &
2416 +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
2422 ! iti=iabs(itype(i,mnum))
2423 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2424 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2425 .and.(mnum.lt.4)) then
2428 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres) &
2429 +d_af_work(ind+j))+sin60*d_as_work(ind+j) &
2430 +cos60*d_as_work1(ind+j))*d_time
2436 end subroutine sddir_verlet2
2437 !-----------------------------------------------------------------------------
2438 subroutine max_accel
2440 ! Find the maximum difference in the accelerations of the the sites
2441 ! at the beginning and the end of the time step.
2444 ! implicit real*8 (a-h,o-z)
2445 ! include 'DIMENSIONS'
2446 ! include 'COMMON.CONTROL'
2447 ! include 'COMMON.VAR'
2448 ! include 'COMMON.MD'
2449 ! include 'COMMON.CHAIN'
2450 ! include 'COMMON.DERIV'
2451 ! include 'COMMON.GEO'
2452 ! include 'COMMON.LOCAL'
2453 ! include 'COMMON.INTERACT'
2454 ! include 'COMMON.IOUNITS'
2455 real(kind=8),dimension(3) :: aux,accel,accel_old
2456 real(kind=8) :: dacc
2460 ! aux(j)=d_a(j,0)-d_a_old(j,0)
2461 accel_old(j)=d_a_old(j,0)
2468 ! 7/3/08 changed to asymmetric difference
2470 ! accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
2471 accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
2472 accel(j)=accel(j)+0.5d0*d_a(j,i)
2473 ! if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
2474 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
2475 dacc=dabs(accel(j)-accel_old(j))
2476 ! write (iout,*) i,dacc
2477 if (dacc.gt.amax) amax=dacc
2485 accel_old(j)=d_a_old(j,0)
2490 accel_old(j)=accel_old(j)+d_a_old(j,1)
2491 accel(j)=accel(j)+d_a(j,1)
2496 ! iti=iabs(itype(i,mnum))
2497 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2498 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2499 .and.(mnum.lt.4)) then
2501 ! accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
2502 accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
2503 accel(j)=accel(j)+d_a(j,i+nres)
2507 ! if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
2508 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
2509 dacc=dabs(accel(j)-accel_old(j))
2510 ! write (iout,*) "side-chain",i,dacc
2511 if (dacc.gt.amax) amax=dacc
2515 accel_old(j)=accel_old(j)+d_a_old(j,i)
2516 accel(j)=accel(j)+d_a(j,i)
2517 ! aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
2521 end subroutine max_accel
2522 !-----------------------------------------------------------------------------
2523 subroutine predict_edrift(epdrift)
2525 ! Predict the drift of the potential energy
2528 use control_data, only: lmuca
2529 ! implicit real*8 (a-h,o-z)
2530 ! include 'DIMENSIONS'
2531 ! include 'COMMON.CONTROL'
2532 ! include 'COMMON.VAR'
2533 ! include 'COMMON.MD'
2534 ! include 'COMMON.CHAIN'
2535 ! include 'COMMON.DERIV'
2536 ! include 'COMMON.GEO'
2537 ! include 'COMMON.LOCAL'
2538 ! include 'COMMON.INTERACT'
2539 ! include 'COMMON.IOUNITS'
2540 ! include 'COMMON.MUCA'
2541 real(kind=8) :: epdrift,epdriftij
2543 ! Drift of the potential energy
2549 epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
2550 if (lmuca) epdriftij=epdriftij*factor
2551 ! write (iout,*) "back",i,j,epdriftij
2552 if (epdriftij.gt.epdrift) epdrift=epdriftij
2556 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1.and.&
2557 molnum(i).lt.4) then
2560 dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
2561 if (lmuca) epdriftij=epdriftij*factor
2562 ! write (iout,*) "side",i,j,epdriftij
2563 if (epdriftij.gt.epdrift) epdrift=epdriftij
2567 epdrift=0.5d0*epdrift*d_time*d_time
2568 ! write (iout,*) "epdrift",epdrift
2570 end subroutine predict_edrift
2571 !-----------------------------------------------------------------------------
2572 subroutine verlet_bath
2574 ! Coupling to the thermostat by using the Berendsen algorithm
2577 ! implicit real*8 (a-h,o-z)
2578 ! include 'DIMENSIONS'
2579 ! include 'COMMON.CONTROL'
2580 ! include 'COMMON.VAR'
2581 ! include 'COMMON.MD'
2582 ! include 'COMMON.CHAIN'
2583 ! include 'COMMON.DERIV'
2584 ! include 'COMMON.GEO'
2585 ! include 'COMMON.LOCAL'
2586 ! include 'COMMON.INTERACT'
2587 ! include 'COMMON.IOUNITS'
2588 ! include 'COMMON.NAMES'
2589 real(kind=8) :: T_half,fact
2590 integer :: i,j,inres,mnum
2592 T_half=2.0d0/(dimen3*Rb)*EK
2593 fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
2594 ! write(iout,*) "T_half", T_half
2595 ! write(iout,*) "EK", EK
2596 ! write(iout,*) "fact", fact
2598 d_t(j,0)=fact*d_t(j,0)
2602 d_t(j,i)=fact*d_t(j,i)
2607 ! iti=iabs(itype(i,mnum))
2608 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2609 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2610 .and.(mnum.lt.4)) then
2613 d_t(j,inres)=fact*d_t(j,inres)
2618 end subroutine verlet_bath
2619 !-----------------------------------------------------------------------------
2621 ! Set up the initial conditions of a MD simulation
2624 use control, only:tcpu
2625 !el use io_basic, only:ilen
2628 use minimm, only:minim_dc,minimize,sc_move
2629 use io_config, only:readrst
2630 use io, only:statout
2631 use random, only: iran_num
2632 ! implicit real*8 (a-h,o-z)
2633 ! include 'DIMENSIONS'
2636 character(len=16) :: form
2637 integer :: IERROR,ERRCODE
2639 integer :: iranmin,itrial,itmp,n_model_try,k, &
2641 integer, dimension(:),allocatable :: list_model_try
2642 integer, dimension(0:nodes-1) :: i_start_models
2643 ! include 'COMMON.SETUP'
2644 ! include 'COMMON.CONTROL'
2645 ! include 'COMMON.VAR'
2646 ! include 'COMMON.MD'
2648 ! include 'COMMON.LANGEVIN'
2650 ! include 'COMMON.LANGEVIN.lang0'
2652 ! include 'COMMON.CHAIN'
2653 ! include 'COMMON.DERIV'
2654 ! include 'COMMON.GEO'
2655 ! include 'COMMON.LOCAL'
2656 ! include 'COMMON.INTERACT'
2657 ! include 'COMMON.IOUNITS'
2658 ! include 'COMMON.NAMES'
2659 ! include 'COMMON.REMD'
2660 real(kind=8),dimension(0:n_ene) :: energia_long,energia_short,energia
2661 real(kind=8),dimension(3) :: vcm,incr,L
2662 real(kind=8) :: xv,sigv,lowb,highb
2663 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
2664 character(len=256) :: qstr
2667 character(len=50) :: tytul
2668 logical :: file_exist
2669 !el common /gucio/ cm
2670 integer :: i,j,ipos,iq,iw,nft_sc,iretcode,ierr,mnum,itime
2674 real(kind=8) :: etot,tt0
2678 ! write(iout,*) "d_time", d_time
2679 ! Compute the standard deviations of stochastic forces for Langevin dynamics
2680 ! if the friction coefficients do not depend on surface area
2681 if (lang.gt.0 .and. .not.surfarea) then
2684 stdforcp(i)=stdfp(mnum)*dsqrt(gamp(mnum))
2688 stdforcsc(i)=stdfsc(iabs(itype(i,mnum)),mnum) &
2689 *dsqrt(gamsc(iabs(itype(i,mnum)),mnum))
2693 ! Open the pdb file for snapshotshots
2696 if (ilen(tmpdir).gt.0) &
2697 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2698 liczba(:ilen(liczba))//".pdb")
2700 file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2704 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) &
2705 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2706 liczba(:ilen(liczba))//".x")
2707 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2710 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) &
2711 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2712 liczba(:ilen(liczba))//".cx")
2713 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2719 if (ilen(tmpdir).gt.0) &
2720 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
2721 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
2723 if (ilen(tmpdir).gt.0) &
2724 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
2725 cartname=prefix(:ilen(prefix))//"_MD.cx"
2729 write (qstr,'(256(1h ))')
2732 iq = qinfrag(i,iset)*10
2733 iw = wfrag(i,iset)/100
2735 if(me.eq.king.or..not.out1file) &
2736 write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
2737 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
2742 iq = qinpair(i,iset)*10
2743 iw = wpair(i,iset)/100
2745 if(me.eq.king.or..not.out1file) &
2746 write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
2747 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
2751 ! pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
2753 ! cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
2755 ! cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
2757 ! statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
2761 if (restart1file) then
2763 inquire(file=mremd_rst_name,exist=file_exist)
2764 write (*,*) me," Before broadcast: file_exist",file_exist
2766 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,&
2769 write (*,*) me," After broadcast: file_exist",file_exist
2770 ! inquire(file=mremd_rst_name,exist=file_exist)
2771 if(me.eq.king.or..not.out1file) &
2772 write(iout,*) "Initial state read by master and distributed"
2774 if (ilen(tmpdir).gt.0) &
2775 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_' &
2776 //liczba(:ilen(liczba))//'.rst')
2777 inquire(file=rest2name,exist=file_exist)
2780 if(.not.restart1file) then
2781 if(me.eq.king.or..not.out1file) &
2782 write(iout,*) "Initial state will be read from file ",&
2783 rest2name(:ilen(rest2name))
2786 call rescale_weights(t_bath)
2788 if(me.eq.king.or..not.out1file)then
2789 if (restart1file) then
2790 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),&
2793 write(iout,*) "File ",rest2name(:ilen(rest2name)),&
2796 write(iout,*) "Initial velocities randomly generated"
2803 ! Generate initial velocities
2804 if(me.eq.king.or..not.out1file) &
2805 write(iout,*) "Initial velocities randomly generated"
2810 ! rest2name = prefix(:ilen(prefix))//'.rst'
2811 if(me.eq.king.or..not.out1file)then
2812 write (iout,*) "Initial velocities"
2814 write (iout,'(i6,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2815 (d_t(j,i+nres),j=1,3)
2817 ! Zeroing the total angular momentum of the system
2818 write(iout,*) "Calling the zero-angular momentum subroutine"
2821 ! Getting the potential energy and forces and velocities and accelerations
2823 ! write (iout,*) "velocity of the center of the mass:"
2824 ! write (iout,*) (vcm(j),j=1,3)
2826 d_t(j,0)=d_t(j,0)-vcm(j)
2828 ! Removing the velocity of the center of mass
2830 if(me.eq.king.or..not.out1file)then
2831 write (iout,*) "vcm right after adjustment:"
2832 write (iout,*) (vcm(j),j=1,3)
2839 if ((.not.rest).or.(forceminim)) then
2840 if (forceminim) call chainbuild_cart
2842 if(iranconf.ne.0 .or.indpdb.gt.0.and..not.unres_pdb .or.preminim) then
2844 print *, 'Calling OVERLAP_SC'
2845 call overlap_sc(fail)
2846 print *,'after OVERLAP'
2849 print *,'call SC_MOVE'
2850 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
2851 print *,'SC_move',nft_sc,etot
2852 if(me.eq.king.or..not.out1file) &
2853 write(iout,*) 'SC_move',nft_sc,etot
2857 print *, 'Calling MINIM_DC'
2858 call minim_dc(etot,iretcode,nfun)
2860 call geom_to_var(nvar,varia)
2861 print *,'Calling MINIMIZE.'
2862 call minimize(etot,varia,iretcode,nfun)
2863 call var_to_geom(nvar,varia)
2865 write(iout,*) "just before minimin"
2867 if(me.eq.king.or..not.out1file) &
2868 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2870 write(iout,*) "just after minimin"
2872 if(iranconf.ne.0) then
2873 !c 8/22/17 AL Loop to produce a low-energy random conformation
2876 if(me.eq.king.or..not.out1file) &
2877 write (iout,*) 'Calling OVERLAP_SC'
2878 call overlap_sc(fail)
2879 endif !endif overlap
2882 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
2883 print *,'SC_move',nft_sc,etot
2884 if(me.eq.king.or..not.out1file) &
2885 write(iout,*) 'SC_move',nft_sc,etot
2889 print *, 'Calling MINIM_DC'
2890 call minim_dc(etot,iretcode,nfun)
2891 call int_from_cart1(.false.)
2893 call geom_to_var(nvar,varia)
2894 print *,'Calling MINIMIZE.'
2895 call minimize(etot,varia,iretcode,nfun)
2896 call var_to_geom(nvar,varia)
2898 if(me.eq.king.or..not.out1file) &
2899 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2900 write(iout,*) "just after minimin"
2902 if (isnan(etot) .or. etot.gt.4.0d6) then
2903 write (iout,*) "Energy too large",etot, &
2904 " trying another random conformation"
2907 call gen_rand_conf(itmp,*30)
2909 30 write (iout,*) 'Failed to generate random conformation', &
2911 write (*,*) 'Processor:',me, &
2912 ' Failed to generate random conformation',&
2921 write (iout,'(a,i3,a)') 'Processor:',me, &
2922 ' error in generating random conformation.'
2923 write (*,'(a,i3,a)') 'Processor:',me, &
2924 ' error in generating random conformation.'
2927 ! call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
2928 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
2938 write (iout,'(a,i3,a)') 'Processor:',me, &
2939 ' failed to generate a low-energy random conformation.'
2940 write (*,'(a,i3,a,f10.3)') 'Processor:',me, &
2941 ' failed to generate a low-energy random conformation.',etot
2945 ! call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
2946 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
2951 else if (preminim) then
2952 if (start_from_model) then
2956 do while (fail .and. n_model_try.lt.nmodel_start)
2957 write (iout,*) "n_model_try",n_model_try
2959 i_model=iran_num(1,nmodel_start)
2961 if (i_model.eq.list_model_try(k)) exit
2963 if (k.gt.n_model_try) exit
2965 n_model_try=n_model_try+1
2966 list_model_try(n_model_try)=i_model
2967 if (me.eq.king .or. .not. out1file) &
2968 write (iout,*) 'Trying to start from model ',&
2969 pdbfiles_chomo(i_model)(:ilen(pdbfiles_chomo(i_model)))
2972 c(j,i)=chomo(j,i,i_model)
2975 call int_from_cart(.true.,.false.)
2976 call sc_loc_geom(.false.)
2980 dc(j,i)=c(j,i+1)-c(j,i)
2981 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
2986 dc(j,i+nres)=c(j,i+nres)-c(j,i)
2987 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
2990 if (me.eq.king.or..not.out1file) then
2991 write (iout,*) "Energies before removing overlaps"
2992 call etotal(energia(0))
2993 call enerprint(energia(0))
2995 ! Remove SC overlaps if requested
2997 write (iout,*) 'Calling OVERLAP_SC'
2998 call overlap_sc(fail)
3001 "Failed to remove overlap from model",i_model
3005 if (me.eq.king.or..not.out1file) then
3006 write (iout,*) "Energies after removing overlaps"
3007 call etotal(energia(0))
3008 call enerprint(energia(0))
3011 ! Search for better SC rotamers if requested
3013 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
3014 print *,'SC_move',nft_sc,etot
3015 if (me.eq.king.or..not.out1file)&
3016 write(iout,*) 'SC_move',nft_sc,etot
3018 call etotal(energia(0))
3021 call MPI_Gather(i_model,1,MPI_INTEGER,i_start_models(0),&
3022 1,MPI_INTEGER,king,CG_COMM,IERROR)
3023 if (n_model_try.gt.nmodel_start .and.&
3024 (me.eq.king .or. out1file)) then
3026 "All models have irreparable overlaps. Trying randoms starts."
3028 i_model=nmodel_start+1
3032 ! Remove SC overlaps if requested
3034 write (iout,*) 'Calling OVERLAP_SC'
3035 call overlap_sc(fail)
3038 "Failed to remove overlap"
3041 if (me.eq.king.or..not.out1file) then
3042 write (iout,*) "Energies after removing overlaps"
3043 call etotal(energia(0))
3044 call enerprint(energia(0))
3047 ! 8/22/17 AL Minimize initial structure
3049 if (me.eq.king.or..not.out1file) write(iout,*)&
3050 'Minimizing initial PDB structure: Calling MINIM_DC'
3051 call minim_dc(etot,iretcode,nfun)
3053 if (me.eq.king.or..not.out1file)&
3054 write(iout,*) 'LBFGS return code is ',statusbf,' eval ',nfun
3057 call geom_to_var(nvar,varia)
3058 if(me.eq.king.or..not.out1file) write (iout,*)&
3059 'Minimizing initial PDB structure: Calling MINIMIZE.'
3060 call minimize(etot,varia,iretcode,nfun)
3061 call var_to_geom(nvar,varia)
3063 if (me.eq.king.or..not.out1file)&
3064 write(iout,*) 'LBFGS return code is ',statusbf,' eval ',nfun
3065 if(me.eq.king.or..not.out1file)&
3066 write(iout,*) 'LBFGS return code is ',statusbf,' eval ',nfun
3068 if (me.eq.king.or..not.out1file)&
3069 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
3070 if(me.eq.king.or..not.out1file)&
3071 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
3075 if (nmodel_start.gt.0 .and. me.eq.king) then
3076 write (iout,'(a)') "Task Starting model"
3078 if (i_start_models(i).gt.nmodel_start) then
3079 write (iout,'(i4,2x,a)') i,"RANDOM STRUCTURE"
3081 write(iout,'(i4,2x,a)')i,pdbfiles_chomo(i_start_models(i)) &
3082 (:ilen(pdbfiles_chomo(i_start_models(i))))
3087 call chainbuild_cart
3089 write(iout,*) "just after kinetic"
3094 kinetic_T=2.0d0/(dimen3*Rb)*EK
3095 if(me.eq.king.or..not.out1file)then
3096 write(iout,*) "just after verlet_bath"
3106 write(iout,*) "before ETOTAL"
3107 call etotal(potEcomp)
3108 if (large) call enerprint(potEcomp)
3111 t_etotal=t_etotal+MPI_Wtime()-tt0
3113 t_etotal=t_etotal+tcpu()-tt0
3118 write(iout,*) "before lagrangian"
3120 write(iout,*) "before max_accel"
3122 if (amax*d_time .gt. dvmax) then
3123 d_time=d_time*dvmax/amax
3124 if(me.eq.king.or..not.out1file) write (iout,*) &
3125 "Time step reduced to",d_time,&
3126 " because of too large initial acceleration."
3128 if(me.eq.king.or..not.out1file)then
3129 write(iout,*) "Potential energy and its components"
3130 call enerprint(potEcomp)
3131 ! write(iout,*) (potEcomp(i),i=0,n_ene)
3133 potE=potEcomp(0)-potEcomp(51)
3137 if (ntwe.ne.0) call statout(itime)
3138 if(me.eq.king.or..not.out1file) &
3139 write (iout,'(/a/3(a25,1pe14.5/))') "Initial:", &
3140 " Kinetic energy",EK," Potential energy",potE, &
3141 " Total energy",totE," Maximum acceleration ", &
3144 write (iout,*) "Initial coordinates"
3146 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),&
3149 write (iout,*) "Initial dC"
3151 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),&
3152 (dc(j,i+nres),j=1,3)
3154 write (iout,*) "Initial velocities"
3155 write (iout,"(13x,' backbone ',23x,' side chain')")
3157 write (iout,'(i6,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
3158 (d_t(j,i+nres),j=1,3)
3160 write (iout,*) "Initial accelerations"
3162 ! write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
3163 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),&
3164 (d_a(j,i+nres),j=1,3)
3170 d_t_old(j,i)=d_t(j,i)
3171 d_a_old(j,i)=d_a(j,i)
3173 ! write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
3182 call etotal_short(energia_short)
3183 if (large) call enerprint(potEcomp)
3186 t_eshort=t_eshort+MPI_Wtime()-tt0
3188 t_eshort=t_eshort+tcpu()-tt0
3193 if(.not.out1file .and. large) then
3194 write (iout,*) "energia_long",energia_long(0),&
3195 " energia_short",energia_short(0),&
3196 " total",energia_long(0)+energia_short(0)
3197 write (iout,*) "Initial fast-force accelerations"
3199 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
3200 (d_a(j,i+nres),j=1,3)
3203 ! 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
3206 d_a_short(j,i)=d_a(j,i)
3215 call etotal_long(energia_long)
3216 if (large) call enerprint(potEcomp)
3219 t_elong=t_elong+MPI_Wtime()-tt0
3221 t_elong=t_elong+tcpu()-tt0
3226 if(.not.out1file .and. large) then
3227 write (iout,*) "energia_long",energia_long(0)
3228 write (iout,*) "Initial slow-force accelerations"
3230 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
3231 (d_a(j,i+nres),j=1,3)
3235 t_enegrad=t_enegrad+MPI_Wtime()-tt0
3237 t_enegrad=t_enegrad+tcpu()-tt0
3241 end subroutine init_MD
3242 !-----------------------------------------------------------------------------
3243 subroutine random_vel
3245 ! implicit real*8 (a-h,o-z)
3247 use random, only:anorm_distr
3249 ! include 'DIMENSIONS'
3250 ! include 'COMMON.CONTROL'
3251 ! include 'COMMON.VAR'
3252 ! include 'COMMON.MD'
3254 ! include 'COMMON.LANGEVIN'
3256 ! include 'COMMON.LANGEVIN.lang0'
3258 ! include 'COMMON.CHAIN'
3259 ! include 'COMMON.DERIV'
3260 ! include 'COMMON.GEO'
3261 ! include 'COMMON.LOCAL'
3262 ! include 'COMMON.INTERACT'
3263 ! include 'COMMON.IOUNITS'
3264 ! include 'COMMON.NAMES'
3265 ! include 'COMMON.TIME1'
3266 real(kind=8) :: xv,sigv,lowb,highb ,Ek1
3268 integer ichain,n,innt,inct,ibeg,ierr
3269 real(kind=8) ,allocatable, dimension(:):: work
3270 integer,allocatable,dimension(:) :: iwork
3271 ! double precision Ghalf(mmaxres2_chain),Geigen(maxres2_chain),&
3272 ! Gvec(maxres2_chain,maxres2_chain)
3273 ! common /przechowalnia/Ghalf,Geigen,Gvec
3275 ! double precision inertia(maxres2_chain,maxres2_chain)
3280 real(kind=8) ,allocatable, dimension(:) :: xsolv,DML,rs
3281 real(kind=8) :: sumx,Ek2,Ek3,aux,masinv
3283 real(kind=8) ,allocatable, dimension(:) :: rsold
3284 real (kind=8),allocatable,dimension(:,:) :: matold,inertia
3288 integer :: i,j,ii,k,mark,imark,mnum,nres2
3289 integer(kind=8) :: ind
3290 ! Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
3292 ! First generate velocities in the eigenspace of the G matrix
3293 ! write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
3296 if(.not.allocated(work)) then
3297 allocate(work(48*nres))
3298 allocate(iwork(6*nres))
3300 print *,"IN RANDOM VEL"
3302 ! print *,size(ghalf)
3305 write (iout,*) "Random_vel, fivediag"
3307 allocate(inertia(2*nres,2*nres))
3314 write(iout,*), nchain
3318 ! if(.not.allocated(ghalf)) print *,"COCO"
3319 ! if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2))
3321 n=dimen_chain(ichain)
3322 innt=iposd_chain(ichain)
3323 if (molnum(innt).eq.5) go to 137
3324 if(.not.allocated(ghalf)) print *,"COCO"
3325 if(.not.allocated(Ghalf)) allocate(Ghalf(1300*(1300+1)/2))
3329 write (iout,*) "Chain",ichain," n",n," start",innt
3331 if (i.lt.inct-1) then
3332 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i),DU1orig(i),&
3334 else if (i.eq.inct-1) then
3335 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i),DU1orig(i)
3337 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i)
3342 ghalf(ind+1)=dmorig(innt)
3343 ghalf(ind+2)=du1orig(innt)
3344 ghalf(ind+3)=dmorig(innt+1)
3348 write (iout,*) "i",i," ind",ind," indu2",innt+i-2,&
3349 " indu1",innt+i-1," indm",innt+i
3350 ghalf(ind+1)=du2orig(innt-1+i-2)
3351 ghalf(ind+2)=du1orig(innt-1+i-1)
3352 ghalf(ind+3)=dmorig(innt-1+i)
3353 !c write (iout,'(3(a,i2,1x))') "DU2",innt-1+i-2,
3354 !c & "DU1",innt-1+i-1,"DM ",innt-1+i
3362 inertia(i,j)=ghalf(ind)
3363 inertia(j,i)=ghalf(ind)
3368 write (iout,*) "Chain ",ichain," ind",ind," dim",n*(n+1)/2
3369 write (iout,*) "Five-diagonal inertia matrix, lower triangle"
3370 ! call matoutr(n,ghalf)
3372 call gldiag(nres*2,n,n,Ghalf,work,Geigen,Gvec,ierr,iwork)
3374 write (iout,'(//a,i3)')&
3375 "Eigenvectors and eigenvalues of the G matrix chain",ichain
3376 call eigout(n,n,nres*2,nres*2,Gvec,Geigen)
3379 !c check diagonalization
3385 aux=aux+gvec(k,i)*gvec(l,j)*inertia(k,l)
3389 write (iout,*) i,j,aux,geigen(i)
3391 write (iout,*) i,j,aux
3397 write(iout,*) "HERE,",n
3403 if (molnum(innt).eq.5) geigen(i)=1.0/msc(itype(innt+i-1,5),5)
3405 sigv=dsqrt((Rb*t_bath)/geigen(i))
3408 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
3409 EK=EK+0.5d0*geigen(i)*d_t_work_new(ii)**2
3410 ! write (iout,*) "i",i," ii",ii," geigen",geigen(i), &
3411 ! " d_t_work_new",d_t_work_new(ii)
3414 if (molnum(innt).eq.5) then
3420 masinv=1.0/msc(itype(innt+i-1,5),5)
3421 d_t_work(ind)=d_t_work(ind)&
3422 +masinv*d_t_work_new((i-1)*3+k)
3432 d_t_work(ind)=d_t_work(ind)&
3433 +Gvec(i,j)*d_t_work_new((j-1)*3+k)
3443 aux=aux+inertia(i,j)*d_t_work(3*(i-1)+k)*d_t_work(3*(j-1)+k)
3449 !c Transfer to the d_t vector
3450 innt=chain_border(1,ichain)
3451 inct=chain_border(2,ichain)
3453 !c write (iout,*) "ichain",ichain," innt",innt," inct",inct
3457 d_t(j,i)=d_t_work(ind)
3460 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum).and.mnum.le.2) then
3463 d_t(j,i+nres)=d_t_work(ind)
3470 write (iout,*) "Random velocities in the Calpha,SC space"
3473 write (iout,'(a3,1h(,i5,1h),3f10.5,3x,3f10.5)')&
3474 restyp(itype(i,mnum),mnum),i,(d_t(j,i),j=1,3),(d_t(j,i+nres),j=1,3)
3477 call kinetic_CASC(Ek1)
3479 ! Transform the velocities to virtual-bond space
3488 if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum).or.mnum.ge.3) then
3490 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
3494 d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
3495 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
3506 !c d_a(:,0)=d_a(:,1)
3508 !c write (iout,*) "Shifting accelerations"
3510 !c write (iout,*) "ichain",chain_border1(1,ichain)-1,
3511 !c & chain_border1(1,ichain)
3512 d_t(:,chain_border1(1,ichain)-1)=d_t(:,chain_border1(1,ichain))
3513 d_t(:,chain_border1(1,ichain))=0.0d0
3515 !c write (iout,*) "Adding accelerations"
3517 !c write (iout,*) "chain",ichain,chain_border1(1,ichain)-1,
3518 !c & chain_border(2,ichain-1)
3519 d_t(:,chain_border1(1,ichain)-1)=&
3520 d_t(:,chain_border1(1,ichain)-1)+d_t(:,chain_border(2,ichain-1))
3521 d_t(:,chain_border(2,ichain-1))=0.0d0
3524 write (iout,*) "chain",ichain,chain_border1(1,ichain)-1,&
3525 chain_border(2,ichain-1)
3526 d_t(:,chain_border1(1,ichain)-1)=&
3527 d_t(:,chain_border1(1,ichain)-1)+d_t(:,chain_border(2,ichain-1))
3528 d_t(:,chain_border(2,ichain-1))=0.0d0
3533 !c d_t(j,0)=d_t(j,nnt)
3536 innt=chain_border(1,ichain)
3537 inct=chain_border(2,ichain)
3538 !c write (iout,*) "ichain",ichain," innt",innt," inct",inct
3539 !c write (iout,*) "ibeg",ibeg
3541 d_t(j,ibeg)=d_t(j,innt)
3546 if (iabs(itype(i,1).eq.10).or.mnum.ge.3) then
3547 !c write (iout,*) "i",i,(d_t(j,i),j=1,3),(d_t(j,i+1),j=1,3)
3549 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
3553 d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
3554 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
3563 "Random velocities in the virtual-bond-vector space"
3564 write (iout,'(3hORG,1h(,i5,1h),3f10.5)') 0,(d_t(j,0),j=1,3)
3566 write (iout,'(a3,1h(,i5,1h),3f10.5,3x,3f10.5)')&
3567 restyp(itype(i,mnum),mnum),i,(d_t(j,i),j=1,3),(d_t(j,i+nres),j=1,3)
3570 write (iout,*) "Kinetic energy from inertia matrix eigenvalues",&
3573 "Kinetic temperatures from inertia matrix eigenvalues",&
3576 write (iout,*) "Kinetic energy from inertia matrix",Ek3
3577 write (iout,*) "Kinetic temperatures from inertia",&
3580 write (iout,*) "Kinetic energy from velocities in CA-SC space",&
3583 "Kinetic temperatures from velovities in CA-SC space",&
3587 "Kinetic energy from virtual-bond-vector velocities",Ek1
3589 "Kinetic temperature from virtual-bond-vector velocities ",&
3598 sigv=dsqrt((Rb*t_bath)/geigen(i))
3601 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
3603 write (iout,*) "i",i," ii",ii," geigen",geigen(i),&
3604 " d_t_work_new",d_t_work_new(ii)
3615 Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
3618 write (iout,*) "Ek from eigenvectors",Ek1
3619 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
3628 d_t_work(ind)=d_t_work(ind) &
3629 +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
3631 ! write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
3635 ! Transfer to the d_t vector
3637 d_t(j,0)=d_t_work(j)
3643 d_t(j,i)=d_t_work(ind)
3648 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3649 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3650 .and.(mnum.lt.4)) then
3653 d_t(j,i+nres)=d_t_work(ind)
3659 ! write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",&
3660 ! 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
3662 ! write(iout,*) "end init MD"
3665 end subroutine random_vel
3666 !-----------------------------------------------------------------------------
3668 subroutine sd_verlet_p_setup
3669 ! Sets up the parameters of stochastic Verlet algorithm
3670 ! implicit real*8 (a-h,o-z)
3671 ! include 'DIMENSIONS'
3672 use control, only: tcpu
3677 ! include 'COMMON.CONTROL'
3678 ! include 'COMMON.VAR'
3679 ! include 'COMMON.MD'
3681 ! include 'COMMON.LANGEVIN'
3683 ! include 'COMMON.LANGEVIN.lang0'
3685 ! include 'COMMON.CHAIN'
3686 ! include 'COMMON.DERIV'
3687 ! include 'COMMON.GEO'
3688 ! include 'COMMON.LOCAL'
3689 ! include 'COMMON.INTERACT'
3690 ! include 'COMMON.IOUNITS'
3691 ! include 'COMMON.NAMES'
3692 ! include 'COMMON.TIME1'
3693 real(kind=8),dimension(6*nres) :: emgdt !(MAXRES6) maxres6=6*maxres
3694 real(kind=8) :: pterm,vterm,rho,rhoc,vsig
3695 real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
3696 prand_vec,vrand_vec1,vrand_vec2 !(MAXRES6) maxres6=6*maxres
3697 logical :: lprn = .false.
3698 real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
3699 real(kind=8) :: ktm,gdt,egdt,gdt2,gdt3,gdt4,gdt5,gdt6,gdt7,gdt8,&
3701 integer :: i,maxres2
3708 ! AL 8/17/04 Code adapted from tinker
3710 ! Get the frictional and random terms for stochastic dynamics in the
3711 ! eigenspace of mass-scaled UNRES friction matrix
3715 gdt = fricgam(i) * d_time
3717 ! Stochastic dynamics reduces to simple MD for zero friction
3719 if (gdt .le. zero) then
3720 pfric_vec(i) = 1.0d0
3721 vfric_vec(i) = d_time
3722 afric_vec(i) = 0.5d0 * d_time * d_time
3723 prand_vec(i) = 0.0d0
3724 vrand_vec1(i) = 0.0d0
3725 vrand_vec2(i) = 0.0d0
3727 ! Analytical expressions when friction coefficient is large
3730 if (gdt .ge. gdt_radius) then
3733 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
3734 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
3735 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
3736 vterm = 1.0d0 - egdt**2
3737 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
3739 ! Use series expansions when friction coefficient is small
3750 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0 &
3751 - gdt5/120.0d0 + gdt6/720.0d0 &
3752 - gdt7/5040.0d0 + gdt8/40320.0d0 &
3753 - gdt9/362880.0d0) / fricgam(i)**2
3754 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
3755 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
3756 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0 &
3757 + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0 &
3758 + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0 &
3759 + 127.0d0*gdt9/90720.0d0
3760 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0 &
3761 - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0 &
3762 - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0 &
3763 - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
3764 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0 &
3765 - 17.0d0*gdt2/1280.0d0 &
3766 + 17.0d0*gdt3/6144.0d0 &
3767 + 40967.0d0*gdt4/34406400.0d0 &
3768 - 57203.0d0*gdt5/275251200.0d0 &
3769 - 1429487.0d0*gdt6/13212057600.0d0)
3772 ! Compute the scaling factors of random terms for the nonzero friction case
3774 ktm = 0.5d0*d_time/fricgam(i)
3775 psig = dsqrt(ktm*pterm) / fricgam(i)
3776 vsig = dsqrt(ktm*vterm)
3777 rhoc = dsqrt(1.0d0 - rho*rho)
3779 vrand_vec1(i) = vsig * rho
3780 vrand_vec2(i) = vsig * rhoc
3785 "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",&
3788 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),&
3789 afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
3793 ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
3796 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
3797 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
3798 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
3799 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
3800 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
3801 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
3804 t_sdsetup=t_sdsetup+MPI_Wtime()
3806 t_sdsetup=t_sdsetup+tcpu()-tt0
3809 end subroutine sd_verlet_p_setup
3810 !-----------------------------------------------------------------------------
3811 subroutine eigtransf1(n,ndim,ab,d,c)
3815 real(kind=8) :: ab(ndim,ndim,n),c(ndim,n),d(ndim)
3821 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
3826 end subroutine eigtransf1
3827 !-----------------------------------------------------------------------------
3828 subroutine eigtransf(n,ndim,a,b,d,c)
3832 real(kind=8) :: a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
3838 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
3843 end subroutine eigtransf
3844 !-----------------------------------------------------------------------------
3845 subroutine sd_verlet1
3847 ! Applying stochastic velocity Verlet algorithm - step 1 to velocities
3849 ! implicit real*8 (a-h,o-z)
3850 ! include 'DIMENSIONS'
3851 ! include 'COMMON.CONTROL'
3852 ! include 'COMMON.VAR'
3853 ! include 'COMMON.MD'
3855 ! include 'COMMON.LANGEVIN'
3857 ! include 'COMMON.LANGEVIN.lang0'
3859 ! include 'COMMON.CHAIN'
3860 ! include 'COMMON.DERIV'
3861 ! include 'COMMON.GEO'
3862 ! include 'COMMON.LOCAL'
3863 ! include 'COMMON.INTERACT'
3864 ! include 'COMMON.IOUNITS'
3865 ! include 'COMMON.NAMES'
3866 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
3867 !el common /stochcalc/ stochforcvec
3868 logical :: lprn = .false.
3869 real(kind=8) :: ddt1,ddt2
3870 integer :: i,j,ind,inres
3872 ! write (iout,*) "dc_old"
3874 ! write (iout,'(i5,3f10.5,5x,3f10.5)')
3875 ! & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3878 dc_work(j)=dc_old(j,0)
3879 d_t_work(j)=d_t_old(j,0)
3880 d_a_work(j)=d_a_old(j,0)
3885 dc_work(ind+j)=dc_old(j,i)
3886 d_t_work(ind+j)=d_t_old(j,i)
3887 d_a_work(ind+j)=d_a_old(j,i)
3893 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3894 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3895 .and.(mnum.lt.4)) then
3897 dc_work(ind+j)=dc_old(j,i+nres)
3898 d_t_work(ind+j)=d_t_old(j,i+nres)
3899 d_a_work(ind+j)=d_a_old(j,i+nres)
3907 "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",&
3911 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
3912 vfric_mat(i,j),afric_mat(i,j),&
3913 prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3921 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) &
3922 +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3923 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3924 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3926 d_t_work_new(i)=ddt1+0.5d0*ddt2
3927 d_t_work(i)=ddt1+ddt2
3932 d_t(j,0)=d_t_work(j)
3937 dc(j,i)=dc_work(ind+j)
3938 d_t(j,i)=d_t_work(ind+j)
3944 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3945 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3946 .and.(mnum.lt.4)) then
3949 dc(j,inres)=dc_work(ind+j)
3950 d_t(j,inres)=d_t_work(ind+j)
3956 end subroutine sd_verlet1
3957 !-----------------------------------------------------------------------------
3958 subroutine sd_verlet2
3960 ! Calculating the adjusted velocities for accelerations
3962 ! implicit real*8 (a-h,o-z)
3963 ! include 'DIMENSIONS'
3964 ! include 'COMMON.CONTROL'
3965 ! include 'COMMON.VAR'
3966 ! include 'COMMON.MD'
3968 ! include 'COMMON.LANGEVIN'
3970 ! include 'COMMON.LANGEVIN.lang0'
3972 ! include 'COMMON.CHAIN'
3973 ! include 'COMMON.DERIV'
3974 ! include 'COMMON.GEO'
3975 ! include 'COMMON.LOCAL'
3976 ! include 'COMMON.INTERACT'
3977 ! include 'COMMON.IOUNITS'
3978 ! include 'COMMON.NAMES'
3979 !el real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV !(MAXRES6) maxres6=6*maxres
3980 real(kind=8),dimension(6*nres) :: stochforcvecV !(MAXRES6) maxres6=6*maxres
3981 !el common /stochcalc/ stochforcvec
3983 real(kind=8) :: ddt1,ddt2
3984 integer :: i,j,ind,inres
3985 ! Compute the stochastic forces which contribute to velocity change
3987 call stochastic_force(stochforcvecV)
3994 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3995 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+ &
3996 vrand_mat2(i,j)*stochforcvecV(j)
3998 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
4002 d_t(j,0)=d_t_work(j)
4007 d_t(j,i)=d_t_work(ind+j)
4013 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
4014 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4015 .and.(mnum.lt.4)) then
4018 d_t(j,inres)=d_t_work(ind+j)
4024 end subroutine sd_verlet2
4025 !-----------------------------------------------------------------------------
4026 subroutine sd_verlet_ciccotti_setup
4028 ! Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
4030 ! implicit real*8 (a-h,o-z)
4031 ! include 'DIMENSIONS'
4032 use control, only: tcpu
4037 ! include 'COMMON.CONTROL'
4038 ! include 'COMMON.VAR'
4039 ! include 'COMMON.MD'
4041 ! include 'COMMON.LANGEVIN'
4043 ! include 'COMMON.LANGEVIN.lang0'
4045 ! include 'COMMON.CHAIN'
4046 ! include 'COMMON.DERIV'
4047 ! include 'COMMON.GEO'
4048 ! include 'COMMON.LOCAL'
4049 ! include 'COMMON.INTERACT'
4050 ! include 'COMMON.IOUNITS'
4051 ! include 'COMMON.NAMES'
4052 ! include 'COMMON.TIME1'
4053 real(kind=8),dimension(6*nres) :: emgdt !(MAXRES6) maxres6=6*maxres
4054 real(kind=8) :: pterm,vterm,rho,rhoc,vsig
4055 real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
4056 prand_vec,vrand_vec1,vrand_vec2 !(MAXRES6) maxres6=6*maxres
4057 logical :: lprn = .false.
4058 real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
4059 real(kind=8) :: ktm,gdt,egdt,tt0
4060 integer :: i,maxres2
4067 ! AL 8/17/04 Code adapted from tinker
4069 ! Get the frictional and random terms for stochastic dynamics in the
4070 ! eigenspace of mass-scaled UNRES friction matrix
4074 write (iout,*) "i",i," fricgam",fricgam(i)
4075 gdt = fricgam(i) * d_time
4077 ! Stochastic dynamics reduces to simple MD for zero friction
4079 if (gdt .le. zero) then
4080 pfric_vec(i) = 1.0d0
4081 vfric_vec(i) = d_time
4082 afric_vec(i) = 0.5d0*d_time*d_time
4083 prand_vec(i) = afric_vec(i)
4084 vrand_vec2(i) = vfric_vec(i)
4086 ! Analytical expressions when friction coefficient is large
4091 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
4092 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
4093 prand_vec(i) = afric_vec(i)
4094 vrand_vec2(i) = vfric_vec(i)
4096 ! Compute the scaling factors of random terms for the nonzero friction case
4098 ! ktm = 0.5d0*d_time/fricgam(i)
4099 ! psig = dsqrt(ktm*pterm) / fricgam(i)
4100 ! vsig = dsqrt(ktm*vterm)
4101 ! prand_vec(i) = psig*afric_vec(i)
4102 ! vrand_vec2(i) = vsig*vfric_vec(i)
4107 "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",&
4110 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),&
4111 afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
4115 ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
4117 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
4118 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
4119 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
4120 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
4121 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
4123 t_sdsetup=t_sdsetup+MPI_Wtime()
4125 t_sdsetup=t_sdsetup+tcpu()-tt0
4128 end subroutine sd_verlet_ciccotti_setup
4129 !-----------------------------------------------------------------------------
4130 subroutine sd_verlet1_ciccotti
4132 ! Applying stochastic velocity Verlet algorithm - step 1 to velocities
4133 ! implicit real*8 (a-h,o-z)
4135 ! include 'DIMENSIONS'
4139 ! include 'COMMON.CONTROL'
4140 ! include 'COMMON.VAR'
4141 ! include 'COMMON.MD'
4143 ! include 'COMMON.LANGEVIN'
4145 ! include 'COMMON.LANGEVIN.lang0'
4147 ! include 'COMMON.CHAIN'
4148 ! include 'COMMON.DERIV'
4149 ! include 'COMMON.GEO'
4150 ! include 'COMMON.LOCAL'
4151 ! include 'COMMON.INTERACT'
4152 ! include 'COMMON.IOUNITS'
4153 ! include 'COMMON.NAMES'
4154 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
4155 !el common /stochcalc/ stochforcvec
4156 logical :: lprn = .false.
4157 real(kind=8) :: ddt1,ddt2
4158 integer :: i,j,ind,inres
4159 ! write (iout,*) "dc_old"
4161 ! write (iout,'(i5,3f10.5,5x,3f10.5)')
4162 ! & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
4165 dc_work(j)=dc_old(j,0)
4166 d_t_work(j)=d_t_old(j,0)
4167 d_a_work(j)=d_a_old(j,0)
4172 dc_work(ind+j)=dc_old(j,i)
4173 d_t_work(ind+j)=d_t_old(j,i)
4174 d_a_work(ind+j)=d_a_old(j,i)
4179 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
4181 dc_work(ind+j)=dc_old(j,i+nres)
4182 d_t_work(ind+j)=d_t_old(j,i+nres)
4183 d_a_work(ind+j)=d_a_old(j,i+nres)
4192 "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",&
4196 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
4197 vfric_mat(i,j),afric_mat(i,j),&
4198 prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
4206 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) &
4207 +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
4208 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
4209 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
4211 d_t_work_new(i)=ddt1+0.5d0*ddt2
4212 d_t_work(i)=ddt1+ddt2
4217 d_t(j,0)=d_t_work(j)
4222 dc(j,i)=dc_work(ind+j)
4223 d_t(j,i)=d_t_work(ind+j)
4229 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
4230 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4231 .and.(mnum.lt.4)) then
4234 dc(j,inres)=dc_work(ind+j)
4235 d_t(j,inres)=d_t_work(ind+j)
4241 end subroutine sd_verlet1_ciccotti
4242 !-----------------------------------------------------------------------------
4243 subroutine sd_verlet2_ciccotti
4245 ! Calculating the adjusted velocities for accelerations
4247 ! implicit real*8 (a-h,o-z)
4248 ! include 'DIMENSIONS'
4249 ! include 'COMMON.CONTROL'
4250 ! include 'COMMON.VAR'
4251 ! include 'COMMON.MD'
4253 ! include 'COMMON.LANGEVIN'
4255 ! include 'COMMON.LANGEVIN.lang0'
4257 ! include 'COMMON.CHAIN'
4258 ! include 'COMMON.DERIV'
4259 ! include 'COMMON.GEO'
4260 ! include 'COMMON.LOCAL'
4261 ! include 'COMMON.INTERACT'
4262 ! include 'COMMON.IOUNITS'
4263 ! include 'COMMON.NAMES'
4264 !el real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV !(MAXRES6) maxres6=6*maxres
4265 real(kind=8),dimension(6*nres) :: stochforcvecV !(MAXRES6) maxres6=6*maxres
4266 !el common /stochcalc/ stochforcvec
4267 real(kind=8) :: ddt1,ddt2
4268 integer :: i,j,ind,inres
4270 ! Compute the stochastic forces which contribute to velocity change
4272 call stochastic_force(stochforcvecV)
4279 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
4280 ! ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
4281 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
4283 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
4287 d_t(j,0)=d_t_work(j)
4292 d_t(j,i)=d_t_work(ind+j)
4298 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4300 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
4303 d_t(j,inres)=d_t_work(ind+j)
4309 end subroutine sd_verlet2_ciccotti
4311 !-----------------------------------------------------------------------------
4313 !-----------------------------------------------------------------------------
4315 subroutine inertia_tensor
4318 real(kind=8) Im(3,3),Imcp(3,3),pr(3),M_SC,&
4319 eigvec(3,3),Id(3,3),eigval(3),L(3),vp(3),vrot(3),&
4320 vpp(3,0:MAXRES),vs_p(3),pr1(3,3),&
4321 pr2(3,3),pp(3),incr(3),v(3),mag,mag2,M_PEP
4322 integer iti,inres,i,j,k,mnum,mnum1
4335 !c caulating the center of the mass of the protein
4339 if (itype(i,mnum).eq.ntyp1_molec(mnum)&
4340 .or. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) cycle
4341 ! if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
4342 if (mnum.ge.5) mp(mnum)=0.0d0
4343 M_PEP=M_PEP+mp(mnum)
4346 cm(j)=cm(j)+(c(j,i)+0.5d0*dc(j,i))*mp(mnum)
4356 iti=iabs(itype(i,mnum))
4357 if (iti.eq.ntyp1_molec(mnum)) cycle
4358 M_SC=M_SC+msc(iabs(iti),mnum)
4361 cm(j)=cm(j)+msc(iabs(iti),mnum)*c(j,inres)
4365 cm(j)=cm(j)/(M_SC+M_PEP)
4370 ! if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
4371 if (mnum.ge.5) mp(mnum)=0.0d0
4372 if (itype(i,mnum).eq.ntyp1_molec(mnum)&
4373 .or. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) cycle
4375 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
4377 Im(1,1)=Im(1,1)+mp(mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
4378 Im(1,2)=Im(1,2)-mp(mnum)*pr(1)*pr(2)
4379 Im(1,3)=Im(1,3)-mp(mnum)*pr(1)*pr(3)
4380 Im(2,3)=Im(2,3)-mp(mnum)*pr(2)*pr(3)
4381 Im(2,2)=Im(2,2)+mp(mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
4382 Im(3,3)=Im(3,3)+mp(mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
4387 iti=iabs(itype(i,mnum))
4388 if (iti.eq.ntyp1_molec(mnum)) cycle
4391 pr(j)=c(j,inres)-cm(j)
4393 Im(1,1)=Im(1,1)+msc(iabs(iti),mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
4394 Im(1,2)=Im(1,2)-msc(iabs(iti),mnum)*pr(1)*pr(2)
4395 Im(1,3)=Im(1,3)-msc(iabs(iti),mnum)*pr(1)*pr(3)
4396 Im(2,3)=Im(2,3)-msc(iabs(iti),mnum)*pr(2)*pr(3)
4397 Im(2,2)=Im(2,2)+msc(iabs(iti),mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
4398 Im(3,3)=Im(3,3)+msc(iabs(iti),mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
4403 if (itype(i,mnum).eq.ntyp1_molec(mnum)&
4404 .or. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) cycle
4405 Im(1,1)=Im(1,1)+Ip(mnum)*(1-dc_norm(1,i)*dc_norm(1,i))*&
4406 vbld(i+1)*vbld(i+1)*0.25d0
4407 Im(1,2)=Im(1,2)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(2,i))*&
4408 vbld(i+1)*vbld(i+1)*0.25d0
4409 Im(1,3)=Im(1,3)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(3,i))*&
4410 vbld(i+1)*vbld(i+1)*0.25d0
4411 Im(2,3)=Im(2,3)+Ip(mnum)*(-dc_norm(2,i)*dc_norm(3,i))*&
4412 vbld(i+1)*vbld(i+1)*0.25d0
4413 Im(2,2)=Im(2,2)+Ip(mnum)*(1-dc_norm(2,i)*dc_norm(2,i))*&
4414 vbld(i+1)*vbld(i+1)*0.25d0
4415 Im(3,3)=Im(3,3)+Ip(mnum)*(1-dc_norm(3,i)*dc_norm(3,i))*&
4416 vbld(i+1)*vbld(i+1)*0.25d0
4421 iti=iabs(itype(i,mnum))
4422 if (iti.ne.10 .and. iti.ne.ntyp1_molec(mnum).and.mnum.le.2) then
4424 Im(1,1)=Im(1,1)+Isc(iti,mnum)*(1-dc_norm(1,inres)*&
4425 dc_norm(1,inres))*vbld(inres)*vbld(inres)
4426 Im(1,2)=Im(1,2)-Isc(iti,mnum)*(dc_norm(1,inres)*&
4427 dc_norm(2,inres))*vbld(inres)*vbld(inres)
4428 Im(1,3)=Im(1,3)-Isc(iti,mnum)*(dc_norm(1,inres)*&
4429 dc_norm(3,inres))*vbld(inres)*vbld(inres)
4430 Im(2,3)=Im(2,3)-Isc(iti,mnum)*(dc_norm(2,inres)*&
4431 dc_norm(3,inres))*vbld(inres)*vbld(inres)
4432 Im(2,2)=Im(2,2)+Isc(iti,mnum)*(1-dc_norm(2,inres)*&
4433 dc_norm(2,inres))*vbld(inres)*vbld(inres)
4434 Im(3,3)=Im(3,3)+Isc(iti,mnum)*(1-dc_norm(3,inres)*&
4435 dc_norm(3,inres))*vbld(inres)*vbld(inres)
4444 !c Copng the Im matrix for the djacob subroutine
4451 !c Finding the eigenvectors and eignvalues of the inertia tensor
4452 call djacob(3,3,10000,1.0d-10,Imcp,eigvec,eigval)
4454 if (dabs(eigval(i)).gt.1.0d-15) then
4455 Id(i,i)=1.0d0/eigval(i)
4462 Imcp(i,j)=eigvec(j,i)
4468 pr1(i,j)=pr1(i,j)+Id(i,k)*Imcp(k,j)
4475 pr2(i,j)=pr2(i,j)+eigvec(i,k)*pr1(k,j)
4479 !c Calculating the total rotational velocity of the molecule
4482 vrot(i)=vrot(i)+pr2(i,j)*L(j)
4488 ! write (iout,*) itype(i+1,mnum1),itype(i,mnum)
4489 if (itype(i+1,mnum1).ne.ntyp1_molec(mnum1) &
4490 .and. itype(i,mnum).eq.ntyp1_molec(mnum) .or.&
4491 itype(i,mnum).ne.ntyp1_molec(mnum) &
4492 .and. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) cycle
4493 call vecpr(vrot(1),dc(1,i),vp)
4495 d_t(j,i)=d_t(j,i)-vp(j)
4500 if(itype(i,1).ne.10 .and.itype(i,mnum).ne.ntyp1_molec(mnum)&
4501 .and.mnum.le.2) then
4503 call vecpr(vrot(1),dc(1,inres),vp)
4505 d_t(j,inres)=d_t(j,inres)-vp(j)
4511 end subroutine inertia_tensor
4512 !------------------------------------------------------------
4513 subroutine angmom(cm,L)
4515 double precision L(3),cm(3),pr(3),vp(3),vrot(3),incr(3),v(3),&
4517 integer iti,inres,i,j,mnum,mnum1
4518 !c Calculate the angular momentum
4528 ! if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
4529 if (mnum.ge.5) mp(mnum)=0.0d0
4530 if (itype(i,mnum).eq.ntyp1_molec(mnum)&
4531 .or. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) cycle
4533 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
4536 v(j)=incr(j)+0.5d0*d_t(j,i)
4539 incr(j)=incr(j)+d_t(j,i)
4541 call vecpr(pr(1),v(1),vp)
4543 L(j)=L(j)+mp(mnum)*vp(j)
4547 pp(j)=0.5d0*d_t(j,i)
4549 call vecpr(pr(1),pp(1),vp)
4551 L(j)=L(j)+Ip(mnum)*vp(j)
4559 iti=iabs(itype(i,mnum))
4562 pr(j)=c(j,inres)-cm(j)
4564 if (itype(i,1).ne.10 .and.itype(i,mnum).ne.ntyp1_molec(mnum)&
4565 .and.mnum.le.2) then
4567 v(j)=incr(j)+d_t(j,inres)
4574 call vecpr(pr(1),v(1),vp)
4575 !c write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),
4576 !c & " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
4577 ! if (mnum.gt.4) then
4583 L(j)=L(j)+mscab*vp(j)
4585 !c write (iout,*) "L",(l(j),j=1,3)
4586 if (itype(i,1).ne.10 .and.itype(i,mnum).ne.ntyp1_molec(mnum)&
4587 .and.mnum.le.2) then
4589 v(j)=incr(j)+d_t(j,inres)
4591 call vecpr(dc(1,inres),d_t(1,inres),vp)
4593 L(j)=L(j)+Isc(iti,mnum)*vp(j)
4597 incr(j)=incr(j)+d_t(j,i)
4601 end subroutine angmom
4602 !---------------------------------------------------------------
4603 subroutine vcm_vel(vcm)
4604 double precision vcm(3),vv(3),summas,amas
4605 integer i,j,mnum,mnum1
4613 if ((mnum.ge.5).or.(mnum.eq.3))&
4614 mp(mnum)=msc(itype(i,mnum),mnum)
4616 summas=summas+mp(mnum)
4618 vcm(j)=vcm(j)+mp(mnum)*(vv(j)+0.5d0*d_t(j,i))
4622 amas=msc(iabs(itype(i,mnum)),mnum)
4626 ! amas=msc(iabs(itype(i)))
4628 ! if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
4629 if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4630 .and.(mnum.lt.4)) then
4632 vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres))
4636 vcm(j)=vcm(j)+amas*vv(j)
4640 vv(j)=vv(j)+d_t(j,i)
4643 !c write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas
4645 vcm(j)=vcm(j)/summas
4648 end subroutine vcm_vel
4650 subroutine inertia_tensor
4652 ! Calculating the intertia tensor for the entire protein in order to
4653 ! remove the perpendicular components of velocity matrix which cause
4654 ! the molecule to rotate.
4657 ! implicit real*8 (a-h,o-z)
4658 ! include 'DIMENSIONS'
4659 ! include 'COMMON.CONTROL'
4660 ! include 'COMMON.VAR'
4661 ! include 'COMMON.MD'
4662 ! include 'COMMON.CHAIN'
4663 ! include 'COMMON.DERIV'
4664 ! include 'COMMON.GEO'
4665 ! include 'COMMON.LOCAL'
4666 ! include 'COMMON.INTERACT'
4667 ! include 'COMMON.IOUNITS'
4668 ! include 'COMMON.NAMES'
4670 real(kind=8),dimension(3,3) :: Im,Imcp,eigvec,Id
4671 real(kind=8),dimension(3) :: pr,eigval,L,vp,vrot
4672 real(kind=8) :: M_SC,mag,mag2,M_PEP
4673 real(kind=8),dimension(3,0:nres) :: vpp !(3,0:MAXRES)
4674 real(kind=8),dimension(3) :: vs_p,pp,incr,v
4675 real(kind=8),dimension(3,3) :: pr1,pr2
4677 !el common /gucio/ cm
4678 integer :: iti,inres,i,j,k,mnum
4689 ! calculating the center of the mass of the protein
4693 if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
4694 write(iout,*) "WTF",itype(i,mnum),i,mnum,mp(mnum)
4695 ! if (itype(i,mnum).eq.ntyp1_molec(mnum)) cycle
4696 M_PEP=M_PEP+mp(mnum)
4699 cm(j)=cm(j)+(c(j,i)+0.5d0*dc(j,i))*mp(mnum)
4708 if (mnum.ge.5) cycle
4709 iti=iabs(itype(i,mnum))
4710 M_SC=M_SC+msc(iabs(iti),mnum)
4713 cm(j)=cm(j)+msc(iabs(iti),mnum)*c(j,inres)
4717 cm(j)=cm(j)/(M_SC+M_PEP)
4719 ! write(iout,*) "Center of mass:",cm
4722 if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
4724 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
4726 Im(1,1)=Im(1,1)+mp(mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
4727 Im(1,2)=Im(1,2)-mp(mnum)*pr(1)*pr(2)
4728 Im(1,3)=Im(1,3)-mp(mnum)*pr(1)*pr(3)
4729 Im(2,3)=Im(2,3)-mp(mnum)*pr(2)*pr(3)
4730 Im(2,2)=Im(2,2)+mp(mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
4731 Im(3,3)=Im(3,3)+mp(mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
4734 ! write(iout,*) "The angular momentum before msc add"
4736 ! write (iout,*) (Im(i,j),j=1,3)
4741 iti=iabs(itype(i,mnum))
4742 if (mnum.ge.5) cycle
4745 pr(j)=c(j,inres)-cm(j)
4747 Im(1,1)=Im(1,1)+msc(iabs(iti),mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
4748 Im(1,2)=Im(1,2)-msc(iabs(iti),mnum)*pr(1)*pr(2)
4749 Im(1,3)=Im(1,3)-msc(iabs(iti),mnum)*pr(1)*pr(3)
4750 Im(2,3)=Im(2,3)-msc(iabs(iti),mnum)*pr(2)*pr(3)
4751 Im(2,2)=Im(2,2)+msc(iabs(iti),mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
4752 Im(3,3)=Im(3,3)+msc(iabs(iti),mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
4754 ! write(iout,*) "The angular momentum before Ip add"
4756 ! write (iout,*) (Im(i,j),j=1,3)
4761 Im(1,1)=Im(1,1)+Ip(mnum)*(1-dc_norm(1,i)*dc_norm(1,i))* &
4762 vbld(i+1)*vbld(i+1)*0.25d0
4763 Im(1,2)=Im(1,2)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(2,i))* &
4764 vbld(i+1)*vbld(i+1)*0.25d0
4765 Im(1,3)=Im(1,3)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(3,i))* &
4766 vbld(i+1)*vbld(i+1)*0.25d0
4767 Im(2,3)=Im(2,3)+Ip(mnum)*(-dc_norm(2,i)*dc_norm(3,i))* &
4768 vbld(i+1)*vbld(i+1)*0.25d0
4769 Im(2,2)=Im(2,2)+Ip(mnum)*(1-dc_norm(2,i)*dc_norm(2,i))* &
4770 vbld(i+1)*vbld(i+1)*0.25d0
4771 Im(3,3)=Im(3,3)+Ip(mnum)*(1-dc_norm(3,i)*dc_norm(3,i))* &
4772 vbld(i+1)*vbld(i+1)*0.25d0
4774 ! write(iout,*) "The angular momentum before Isc add"
4776 ! write (iout,*) (Im(i,j),j=1,3)
4782 ! if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
4783 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4784 .and.(mnum.lt.4)) then
4785 iti=iabs(itype(i,mnum))
4787 Im(1,1)=Im(1,1)+Isc(iti,mnum)*(1-dc_norm(1,inres)* &
4788 dc_norm(1,inres))*vbld(inres)*vbld(inres)
4789 Im(1,2)=Im(1,2)-Isc(iti,mnum)*(dc_norm(1,inres)* &
4790 dc_norm(2,inres))*vbld(inres)*vbld(inres)
4791 Im(1,3)=Im(1,3)-Isc(iti,mnum)*(dc_norm(1,inres)* &
4792 dc_norm(3,inres))*vbld(inres)*vbld(inres)
4793 Im(2,3)=Im(2,3)-Isc(iti,mnum)*(dc_norm(2,inres)* &
4794 dc_norm(3,inres))*vbld(inres)*vbld(inres)
4795 Im(2,2)=Im(2,2)+Isc(iti,mnum)*(1-dc_norm(2,inres)* &
4796 dc_norm(2,inres))*vbld(inres)*vbld(inres)
4797 Im(3,3)=Im(3,3)+Isc(iti,mnum)*(1-dc_norm(3,inres)* &
4798 dc_norm(3,inres))*vbld(inres)*vbld(inres)
4802 ! write(iout,*) "The angular momentum before agnom:"
4804 ! write (iout,*) (Im(i,j),j=1,3)
4808 ! write(iout,*) "The angular momentum before adjustment:"
4809 ! write(iout,*) (L(j),j=1,3)
4811 ! write (iout,*) (Im(i,j),j=1,3)
4817 ! Copying the Im matrix for the djacob subroutine
4825 ! Finding the eigenvectors and eignvalues of the inertia tensor
4826 call djacob(3,3,10000,1.0d-10,Imcp,eigvec,eigval)
4827 ! write (iout,*) "Eigenvalues & Eigenvectors"
4828 ! write (iout,'(5x,3f10.5)') (eigval(i),i=1,3)
4831 ! write (iout,'(i5,3f10.5)') i,(eigvec(i,j),j=1,3)
4833 ! Constructing the diagonalized matrix
4835 if (dabs(eigval(i)).gt.1.0d-15) then
4836 Id(i,i)=1.0d0/eigval(i)
4843 Imcp(i,j)=eigvec(j,i)
4849 pr1(i,j)=pr1(i,j)+Id(i,k)*Imcp(k,j)
4856 pr2(i,j)=pr2(i,j)+eigvec(i,k)*pr1(k,j)
4860 ! Calculating the total rotational velocity of the molecule
4863 vrot(i)=vrot(i)+pr2(i,j)*L(j)
4866 ! Resetting the velocities
4868 call vecpr(vrot(1),dc(1,i),vp)
4870 ! print *,"HERE2",d_t(j,i),vp(j)
4871 d_t(j,i)=d_t(j,i)-vp(j)
4872 ! print *,"HERE2",d_t(j,i),vp(j)
4877 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4878 .and.(mnum.lt.4)) then
4879 ! if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
4880 ! if(itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
4882 call vecpr(vrot(1),dc(1,inres),vp)
4884 d_t(j,inres)=d_t(j,inres)-vp(j)
4889 ! write(iout,*) "The angular momentum after adjustment:"
4890 ! write(iout,*) (L(j),j=1,3)
4893 end subroutine inertia_tensor
4894 !-----------------------------------------------------------------------------
4895 subroutine angmom(cm,L)
4898 ! implicit real*8 (a-h,o-z)
4899 ! include 'DIMENSIONS'
4900 ! include 'COMMON.CONTROL'
4901 ! include 'COMMON.VAR'
4902 ! include 'COMMON.MD'
4903 ! include 'COMMON.CHAIN'
4904 ! include 'COMMON.DERIV'
4905 ! include 'COMMON.GEO'
4906 ! include 'COMMON.LOCAL'
4907 ! include 'COMMON.INTERACT'
4908 ! include 'COMMON.IOUNITS'
4909 ! include 'COMMON.NAMES'
4910 real(kind=8) :: mscab
4911 real(kind=8),dimension(3) :: L,cm,pr,vp,vrot,incr,v,pp
4912 integer :: iti,inres,i,j,mnum
4913 ! Calculate the angular momentum
4922 if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
4924 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
4927 v(j)=incr(j)+0.5d0*d_t(j,i)
4930 incr(j)=incr(j)+d_t(j,i)
4932 call vecpr(pr(1),v(1),vp)
4934 L(j)=L(j)+mp(mnum)*vp(j)
4935 ! print *,"HERE3",J,i,L(j),mp(mnum),Ip(mnum),mnum
4939 pp(j)=0.5d0*d_t(j,i)
4941 call vecpr(pr(1),pp(1),vp)
4944 L(j)=L(j)+Ip(mnum)*vp(j)
4952 iti=iabs(itype(i,mnum))
4960 pr(j)=c(j,inres)-cm(j)
4963 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4964 .and.(mnum.lt.4)) then
4966 v(j)=incr(j)+d_t(j,inres)
4973 ! print *,i,pr,"pr",v
4974 call vecpr(pr(1),v(1),vp)
4975 ! write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),&
4976 ! " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
4978 L(j)=L(j)+mscab*vp(j)
4980 ! write (iout,*) "L",(l(j),j=1,3)
4981 ! print *,"L",(l(j),j=1,3),i,vp(1)
4983 if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4984 .and.(mnum.lt.4)) then
4986 v(j)=incr(j)+d_t(j,inres)
4988 call vecpr(dc(1,inres),d_t(1,inres),vp)
4990 L(j)=L(j)+Isc(iti,mnum)*vp(j)
4994 incr(j)=incr(j)+d_t(j,i)
4998 end subroutine angmom
4999 !-----------------------------------------------------------------------------
5000 subroutine vcm_vel(vcm)
5003 ! implicit real*8 (a-h,o-z)
5004 ! include 'DIMENSIONS'
5005 ! include 'COMMON.VAR'
5006 ! include 'COMMON.MD'
5007 ! include 'COMMON.CHAIN'
5008 ! include 'COMMON.DERIV'
5009 ! include 'COMMON.GEO'
5010 ! include 'COMMON.LOCAL'
5011 ! include 'COMMON.INTERACT'
5012 ! include 'COMMON.IOUNITS'
5013 real(kind=8),dimension(3) :: vcm,vv
5014 real(kind=8) :: summas,amas
5024 if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
5026 summas=summas+mp(mnum)
5028 vcm(j)=vcm(j)+mp(mnum)*(vv(j)+0.5d0*d_t(j,i))
5029 ! print *,i,j,vv(j),d_t(j,i)
5033 amas=msc(iabs(itype(i,mnum)),mnum)
5038 if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5039 .and.(mnum.lt.4)) then
5041 vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres))
5045 vcm(j)=vcm(j)+amas*vv(j)
5049 vv(j)=vv(j)+d_t(j,i)
5052 ! write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas
5054 vcm(j)=vcm(j)/summas
5057 end subroutine vcm_vel
5059 !-----------------------------------------------------------------------------
5061 !-----------------------------------------------------------------------------
5063 ! RATTLE algorithm for velocity Verlet - step 1, UNRES
5067 ! implicit real*8 (a-h,o-z)
5068 ! include 'DIMENSIONS'
5070 ! include 'COMMON.CONTROL'
5071 ! include 'COMMON.VAR'
5072 ! include 'COMMON.MD'
5074 ! include 'COMMON.LANGEVIN'
5076 ! include 'COMMON.LANGEVIN.lang0'
5078 ! include 'COMMON.CHAIN'
5079 ! include 'COMMON.DERIV'
5080 ! include 'COMMON.GEO'
5081 ! include 'COMMON.LOCAL'
5082 ! include 'COMMON.INTERACT'
5083 ! include 'COMMON.IOUNITS'
5084 ! include 'COMMON.NAMES'
5085 ! include 'COMMON.TIME1'
5086 !el real(kind=8) :: gginv(2*nres,2*nres),&
5087 !el gdc(3,2*nres,2*nres)
5088 real(kind=8) :: dC_uncor(3,2*nres) !,&
5089 !el real(kind=8) :: Cmat(2*nres,2*nres)
5090 real(kind=8) :: x(2*nres),xcorr(3,2*nres) !maxres2=2*maxres
5091 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
5092 !el common /przechowalnia/ nbond
5093 integer :: max_rattle = 5
5094 logical :: lprn = .false., lprn1 = .false., not_done
5095 real(kind=8) :: tol_rattle = 1.0d-5
5097 integer :: ii,i,j,jj,l,ind,ind1,nres2
5100 !el /common/ przechowalnia
5102 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
5103 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
5104 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
5106 if (lprn) write (iout,*) "RATTLE1"
5110 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5111 .and.(mnum.lt.4)) nbond=nbond+1
5113 ! Make a folded form of the Ginv-matrix
5126 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
5131 if (itype(k,1).ne.10 .and. itype(k,mnum).ne.ntyp1_molec(mnum)&
5132 .and.(mnum.lt.4)) then
5136 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
5144 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5155 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
5159 if (itype(k,1).ne.10) then
5163 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
5171 write (iout,*) "Matrix GGinv"
5172 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
5178 if (iter.gt.max_rattle) then
5179 write (iout,*) "Error - too many iterations in RATTLE."
5182 ! Calculate the matrix C = GG**(-1) dC_old o dC
5187 dC_uncor(j,ind1)=dC(j,i)
5191 if (itype(i,1).ne.10) then
5194 dC_uncor(j,ind1)=dC(j,i+nres)
5203 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
5207 if (itype(k,1).ne.10) then
5210 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
5215 ! Calculate deviations from standard virtual-bond lengths
5219 x(ind)=vbld(i+1)**2-vbl**2
5222 if (itype(i,1).ne.10) then
5224 x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
5228 write (iout,*) "Coordinates and violations"
5230 write(iout,'(i5,3f10.5,5x,e15.5)') &
5231 i,(dC_uncor(j,i),j=1,3),x(i)
5233 write (iout,*) "Velocities and violations"
5237 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5238 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
5242 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5243 .and.(mnum.lt.4)) then
5246 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5247 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
5248 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
5251 ! write (iout,*) "gdc"
5253 ! write (iout,*) "i",i
5255 ! write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
5261 if (dabs(x(i)).gt.xmax) then
5265 if (xmax.lt.tol_rattle) then
5269 ! Calculate the matrix of the system of equations
5274 Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
5279 write (iout,*) "Matrix Cmat"
5280 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
5282 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
5283 ! Add constraint term to positions
5290 xx = xx+x(ii)*gdc(j,ind,ii)
5294 d_t_new(j,i)=d_t_new(j,i)-xx/d_time
5298 if (itype(i,1).ne.10) then
5303 xx = xx+x(ii)*gdc(j,ind,ii)
5306 dC(j,i+nres)=dC(j,i+nres)-xx
5307 d_t_new(j,i+nres)=d_t_new(j,i+nres)-xx/d_time
5311 ! Rebuild the chain using the new coordinates
5312 call chainbuild_cart
5314 write (iout,*) "New coordinates, Lagrange multipliers,",&
5315 " and differences between actual and standard bond lengths"
5319 xx=vbld(i+1)**2-vbl**2
5320 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
5321 i,(dC(j,i),j=1,3),x(ind),xx
5325 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5328 xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
5329 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
5330 i,(dC(j,i+nres),j=1,3),x(ind),xx
5333 write (iout,*) "Velocities and violations"
5337 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5338 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
5341 if (itype(i,1).ne.10) then
5343 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5344 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
5345 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
5352 10 write (iout,*) "Error - singularity in solving the system",&
5353 " of equations for Lagrange multipliers."
5357 "RATTLE inactive; use -DRATTLE switch at compile time."
5360 end subroutine rattle1
5361 !-----------------------------------------------------------------------------
5363 ! RATTLE algorithm for velocity Verlet - step 2, UNRES
5367 ! implicit real*8 (a-h,o-z)
5368 ! include 'DIMENSIONS'
5370 ! include 'COMMON.CONTROL'
5371 ! include 'COMMON.VAR'
5372 ! include 'COMMON.MD'
5374 ! include 'COMMON.LANGEVIN'
5376 ! include 'COMMON.LANGEVIN.lang0'
5378 ! include 'COMMON.CHAIN'
5379 ! include 'COMMON.DERIV'
5380 ! include 'COMMON.GEO'
5381 ! include 'COMMON.LOCAL'
5382 ! include 'COMMON.INTERACT'
5383 ! include 'COMMON.IOUNITS'
5384 ! include 'COMMON.NAMES'
5385 ! include 'COMMON.TIME1'
5386 !el real(kind=8) :: gginv(2*nres,2*nres),&
5387 !el gdc(3,2*nres,2*nres)
5388 real(kind=8) :: dC_uncor(3,2*nres) !,&
5389 !el Cmat(2*nres,2*nres)
5390 real(kind=8) :: x(2*nres) !maxres2=2*maxres
5391 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
5392 !el common /przechowalnia/ nbond
5393 integer :: max_rattle = 5
5394 logical :: lprn = .false., lprn1 = .false., not_done
5395 real(kind=8) :: tol_rattle = 1.0d-5
5399 !el /common/ przechowalnia
5400 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
5401 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
5402 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
5404 if (lprn) write (iout,*) "RATTLE2"
5405 if (lprn) write (iout,*) "Velocity correction"
5406 ! Calculate the matrix G dC
5412 gdc(j,i,ind)=GGinv(i,ind)*dC(j,k)
5417 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5418 .and.(mnum.lt.4)) then
5421 gdc(j,i,ind)=GGinv(i,ind)*dC(j,k+nres)
5427 ! write (iout,*) "gdc"
5429 ! write (iout,*) "i",i
5431 ! write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
5435 ! Calculate the matrix of the system of equations
5442 Cmat(ind,j)=Cmat(ind,j)+dC(k,i)*gdc(k,ind,j)
5448 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5449 .and.(mnum.lt.4)) then
5454 Cmat(ind,j)=Cmat(ind,j)+dC(k,i+nres)*gdc(k,ind,j)
5459 ! Calculate the scalar product dC o d_t_new
5463 x(ind)=scalar(d_t(1,i),dC(1,i))
5467 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5468 .and.(mnum.lt.4)) then
5470 x(ind)=scalar(d_t(1,i+nres),dC(1,i+nres))
5474 write (iout,*) "Velocities and violations"
5478 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5479 i,ind,(d_t(j,i),j=1,3),x(ind)
5483 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5484 .and.(mnum.lt.4)) then
5486 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5487 i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind)
5493 if (dabs(x(i)).gt.xmax) then
5497 if (xmax.lt.tol_rattle) then
5502 write (iout,*) "Matrix Cmat"
5503 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
5505 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
5506 ! Add constraint term to velocities
5513 xx = xx+x(ii)*gdc(j,ind,ii)
5515 d_t(j,i)=d_t(j,i)-xx
5520 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5521 .and.(mnum.lt.4)) then
5526 xx = xx+x(ii)*gdc(j,ind,ii)
5528 d_t(j,i+nres)=d_t(j,i+nres)-xx
5534 "New velocities, Lagrange multipliers violations"
5538 if (lprn) write (iout,'(2i5,3f10.5,5x,2e15.5)') &
5539 i,ind,(d_t(j,i),j=1,3),x(ind),scalar(d_t(1,i),dC(1,i))
5543 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5546 write (iout,'(2i5,3f10.5,5x,2e15.5)') &
5547 i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind),&
5548 scalar(d_t(1,i+nres),dC(1,i+nres))
5554 10 write (iout,*) "Error - singularity in solving the system",&
5555 " of equations for Lagrange multipliers."
5559 "RATTLE inactive; use -DRATTLE option at compile time."
5562 end subroutine rattle2
5563 !-----------------------------------------------------------------------------
5564 subroutine rattle_brown
5565 ! RATTLE/LINCS algorithm for Brownian dynamics, UNRES
5569 ! implicit real*8 (a-h,o-z)
5570 ! include 'DIMENSIONS'
5572 ! include 'COMMON.CONTROL'
5573 ! include 'COMMON.VAR'
5574 ! include 'COMMON.MD'
5576 ! include 'COMMON.LANGEVIN'
5578 ! include 'COMMON.LANGEVIN.lang0'
5580 ! include 'COMMON.CHAIN'
5581 ! include 'COMMON.DERIV'
5582 ! include 'COMMON.GEO'
5583 ! include 'COMMON.LOCAL'
5584 ! include 'COMMON.INTERACT'
5585 ! include 'COMMON.IOUNITS'
5586 ! include 'COMMON.NAMES'
5587 ! include 'COMMON.TIME1'
5588 !el real(kind=8) :: gginv(2*nres,2*nres),&
5589 !el gdc(3,2*nres,2*nres)
5590 real(kind=8) :: dC_uncor(3,2*nres) !,&
5591 !el real(kind=8) :: Cmat(2*nres,2*nres)
5592 real(kind=8) :: x(2*nres) !maxres2=2*maxres
5593 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
5594 !el common /przechowalnia/ nbond
5595 integer :: max_rattle = 5
5596 logical :: lprn = .false., lprn1 = .false., not_done
5597 real(kind=8) :: tol_rattle = 1.0d-5
5601 !el /common/ przechowalnia
5602 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
5603 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
5604 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
5607 if (lprn) write (iout,*) "RATTLE_BROWN"
5610 if (itype(i,1).ne.10) nbond=nbond+1
5612 ! Make a folded form of the Ginv-matrix
5625 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
5629 if (itype(k,1).ne.10) then
5633 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
5640 if (itype(i,1).ne.10) then
5650 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
5654 if (itype(k,1).ne.10) then
5658 if (j.eq.1 .and. l.eq.1)GGinv(ii,jj)=fricmat(ind,ind1)
5666 write (iout,*) "Matrix GGinv"
5667 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
5673 if (iter.gt.max_rattle) then
5674 write (iout,*) "Error - too many iterations in RATTLE."
5677 ! Calculate the matrix C = GG**(-1) dC_old o dC
5682 dC_uncor(j,ind1)=dC(j,i)
5686 if (itype(i,1).ne.10) then
5689 dC_uncor(j,ind1)=dC(j,i+nres)
5698 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
5702 if (itype(k,1).ne.10) then
5705 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
5710 ! Calculate deviations from standard virtual-bond lengths
5714 x(ind)=vbld(i+1)**2-vbl**2
5717 if (itype(i,1).ne.10) then
5719 x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
5723 write (iout,*) "Coordinates and violations"
5725 write(iout,'(i5,3f10.5,5x,e15.5)') &
5726 i,(dC_uncor(j,i),j=1,3),x(i)
5728 write (iout,*) "Velocities and violations"
5732 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5733 i,ind,(d_t(j,i),j=1,3),scalar(d_t(1,i),dC_old(1,i))
5736 if (itype(i,1).ne.10) then
5738 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5739 i+nres,ind,(d_t(j,i+nres),j=1,3),&
5740 scalar(d_t(1,i+nres),dC_old(1,i+nres))
5743 write (iout,*) "gdc"
5745 write (iout,*) "i",i
5747 write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
5753 if (dabs(x(i)).gt.xmax) then
5757 if (xmax.lt.tol_rattle) then
5761 ! Calculate the matrix of the system of equations
5766 Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
5771 write (iout,*) "Matrix Cmat"
5772 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
5774 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
5775 ! Add constraint term to positions
5782 xx = xx+x(ii)*gdc(j,ind,ii)
5785 d_t(j,i)=d_t(j,i)+xx/d_time
5790 if (itype(i,1).ne.10) then
5795 xx = xx+x(ii)*gdc(j,ind,ii)
5798 d_t(j,i+nres)=d_t(j,i+nres)+xx/d_time
5799 dC(j,i+nres)=dC(j,i+nres)+xx
5803 ! Rebuild the chain using the new coordinates
5804 call chainbuild_cart
5806 write (iout,*) "New coordinates, Lagrange multipliers,",&
5807 " and differences between actual and standard bond lengths"
5811 xx=vbld(i+1)**2-vbl**2
5812 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
5813 i,(dC(j,i),j=1,3),x(ind),xx
5816 if (itype(i,1).ne.10) then
5818 xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
5819 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
5820 i,(dC(j,i+nres),j=1,3),x(ind),xx
5823 write (iout,*) "Velocities and violations"
5827 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5828 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
5831 if (itype(i,1).ne.10) then
5833 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5834 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
5835 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
5842 10 write (iout,*) "Error - singularity in solving the system",&
5843 " of equations for Lagrange multipliers."
5847 "RATTLE inactive; use -DRATTLE option at compile time"
5850 end subroutine rattle_brown
5851 !-----------------------------------------------------------------------------
5853 !-----------------------------------------------------------------------------
5854 subroutine friction_force
5859 ! implicit real*8 (a-h,o-z)
5860 ! include 'DIMENSIONS'
5861 ! include 'COMMON.VAR'
5862 ! include 'COMMON.CHAIN'
5863 ! include 'COMMON.DERIV'
5864 ! include 'COMMON.GEO'
5865 ! include 'COMMON.LOCAL'
5866 ! include 'COMMON.INTERACT'
5867 ! include 'COMMON.MD'
5869 ! include 'COMMON.LANGEVIN'
5871 ! include 'COMMON.LANGEVIN.lang0'
5873 ! include 'COMMON.IOUNITS'
5874 !el real(kind=8),dimension(6*nres) :: gamvec !(MAXRES6) maxres6=6*maxres
5875 !el common /syfek/ gamvec
5877 integer iposc,ichain,n,innt,inct
5878 double precision rs(nres*2)
5879 real(kind=8) ::v_work(3,6*nres),vvec(2*nres)
5881 real(kind=8) :: v_work(6*nres)
5884 real(kind=8) :: vv(3),vvtot(3,nres)!,v_work(6*nres) !,&
5885 !el ginvfric(2*nres,2*nres) !maxres2=2*maxres
5886 !el common /przechowalnia/ ginvfric
5888 logical :: lprn, checkmode
5889 integer :: i,j,ind,k,nres2,nres6,mnum
5894 ! if (large) lprn=.true.
5895 ! if (large) checkmode=.true.
5897 !c Here accelerations due to friction forces are computed right after forces.
5898 d_t_work(:6*nres)=0.0d0
5900 v_work(j,1)=d_t(j,0)
5901 v_work(j,nnt)=d_t(j,0)
5905 v_work(j,i)=v_work(j,i-1)+d_t(j,i-1)
5910 if (iabs(itype(i,1)).ne.10 .and. iabs(itype(i,mnum)).ne.ntyp1_molec(mnum).and.mnum.lt.3) then
5912 v_work(j,i+nres)=v_work(j,i)+d_t(j,i+nres)
5917 write (iout,*) "v_work"
5919 write (iout,'(i5,3f10.5)') i,(v_work(j,i),j=1,3)
5925 n=dimen_chain(ichain)
5926 iposc=iposd_chain(ichain)
5927 !c write (iout,*) "friction_force j",j," ichain",ichain,
5928 !c & " n",n," iposc",iposc,iposc+n-1
5929 innt=chain_border(1,ichain)
5930 inct=chain_border(2,ichain)
5932 !c innt=chain_border(1,1)
5933 !c inct=chain_border(2,1)
5935 vvec(ind+1)=v_work(j,i)
5937 ! if (iabs(itype(i)).ne.10) then
5938 if (iabs(itype(i,1)).ne.10 .and. iabs(itype(i,mnum)).ne.ntyp1_molec(mnum).and.mnum.lt.3) then
5939 vvec(ind+1)=v_work(j,i+nres)
5944 write (iout,*) "vvec ind",ind," n",n
5945 write (iout,'(f10.5)') (vvec(i),i=iposc,ind)
5947 !c write (iout,*) "chain",i," ind",ind," n",n
5955 ! if (large) print *,"before fivediagmult"
5956 call fivediagmult(n,DMfric(iposc),DU1fric(iposc),&
5957 DU2fric(iposc),vvec(iposc),rs)
5958 ! if (large) print *,"after fivediagmult"
5962 time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
5964 time_fricmatmult=time_fricmatmult+tcpu()-time01
5969 write (iout,'(f10.5)') (rs(i),i=1,n)
5971 do i=iposc,iposc+n-1
5972 ! write (iout,*) "ichain",ichain," i",i," j",j,&
5973 ! "index",3*(i-1)+j,"rs",rs(i-iposc+1)
5974 fric_work(3*(i-1)+j)=-rs(i-iposc+1)
5979 write (iout,*) "Vector fric_work dimen3",dimen3
5980 write (iout,'(3f10.5)') (fric_work(j),j=1,dimen3)
5983 if(.not.allocated(gamvec)) allocate(gamvec(nres6)) !(MAXRES6)
5984 if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
5992 d_t_work(j)=d_t(j,0)
5997 d_t_work(ind+j)=d_t(j,i)
6003 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
6004 .and.(mnum.lt.4)) then
6006 d_t_work(ind+j)=d_t(j,i+nres)
6012 call fricmat_mult(d_t_work,fric_work)
6014 if (.not.checkmode) return
6017 write (iout,*) "d_t_work and fric_work"
6019 write (iout,'(i3,2e15.5)') i,d_t_work(i),fric_work(i)
6023 friction(j,0)=fric_work(j)
6028 friction(j,i)=fric_work(ind+j)
6034 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
6035 .and.(mnum.lt.4)) then
6036 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
6038 friction(j,i+nres)=fric_work(ind+j)
6044 write(iout,*) "Friction backbone"
6046 write(iout,'(i5,3e15.5,5x,3e15.5)') &
6047 i,(friction(j,i),j=1,3),(d_t(j,i),j=1,3)
6049 write(iout,*) "Friction side chain"
6051 write(iout,'(i5,3e15.5,5x,3e15.5)') &
6052 i,(friction(j,i+nres),j=1,3),(d_t(j,i+nres),j=1,3)
6061 vvtot(j,i)=vv(j)+0.5d0*d_t(j,i)
6062 vvtot(j,i+nres)=vv(j)+d_t(j,i+nres)
6063 vv(j)=vv(j)+d_t(j,i)
6066 write (iout,*) "vvtot backbone and sidechain"
6068 write (iout,'(i5,3e15.5,5x,3e15.5)') i,(vvtot(j,i),j=1,3),&
6069 (vvtot(j,i+nres),j=1,3)
6074 v_work(ind+j)=vvtot(j,i)
6080 v_work(ind+j)=vvtot(j,i+nres)
6084 write (iout,*) "v_work gamvec and site-based friction forces"
6086 write (iout,'(i5,3e15.5)') i,v_work(i),gamvec(i),&
6090 ! fric_work1(i)=0.0d0
6092 ! fric_work1(i)=fric_work1(i)-A(j,i)*gamvec(j)*v_work(j)
6095 ! write (iout,*) "fric_work and fric_work1"
6097 ! write (iout,'(i5,2e15.5)') i,fric_work(i),fric_work1(i)
6103 ginvfric(i,j)=ginvfric(i,j)+ginv(i,k)*fricmat(k,j)
6107 write (iout,*) "ginvfric"
6109 write (iout,'(i5,100f8.3)') i,(ginvfric(i,j),j=1,dimen)
6111 write (iout,*) "symmetry check"
6114 write (iout,*) i,j,ginvfric(i,j)-ginvfric(j,i)
6120 end subroutine friction_force
6121 !-----------------------------------------------------------------------------
6122 subroutine setup_fricmat
6126 use control_data, only:time_Bcast
6127 use control, only:tcpu
6129 ! implicit real*8 (a-h,o-z)
6133 real(kind=8) :: time00
6135 ! include 'DIMENSIONS'
6136 ! include 'COMMON.VAR'
6137 ! include 'COMMON.CHAIN'
6138 ! include 'COMMON.DERIV'
6139 ! include 'COMMON.GEO'
6140 ! include 'COMMON.LOCAL'
6141 ! include 'COMMON.INTERACT'
6142 ! include 'COMMON.MD'
6143 ! include 'COMMON.SETUP'
6144 ! include 'COMMON.TIME1'
6145 ! integer licznik /0/
6148 ! include 'COMMON.LANGEVIN'
6150 ! include 'COMMON.LANGEVIN.lang0'
6152 ! include 'COMMON.IOUNITS'
6154 integer :: i,j,ind,ind1,m,ichain,innt,inct
6155 logical :: lprn = .false.
6156 real(kind=8) :: dtdi !el ,gamvec(2*nres)
6157 !el real(kind=8),dimension(2*nres,2*nres) :: ginvfric,fcopy
6158 ! real(kind=8),allocatable,dimension(:,:) :: fcopy
6159 !el real(kind=8),dimension(2*nres*(2*nres+1)/2) :: Ghalf !(mmaxres2) (mmaxres2=(maxres2*(maxres2+1)/2))
6160 !el common /syfek/ gamvec
6161 real(kind=8) :: work(8*2*nres)
6162 integer :: iwork(2*nres)
6163 !el common /przechowalnia/ ginvfric,Ghalf,fcopy
6164 integer :: ii,iti,k,l,nzero,nres2,nres6,ierr,mnum
6169 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
6170 if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
6171 if (fg_rank.ne.king) goto 10
6177 if(.not.allocated(gamvec)) allocate(gamvec(nres2)) !(MAXRES2)
6179 if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
6180 if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
6181 !el allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
6182 if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !maxres2=2*maxres
6184 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
6187 if (.not.allocated(DMfric)) allocate(DMfric(nres2))
6188 if (.not.allocated(DU1fric)) allocate(DU1fric(nres2))
6189 if (.not.allocated(DU2fric)) allocate(DU2fric(nres2))
6190 ! Load the friction coefficients corresponding to peptide groups
6195 gamvec(ind1)=gamp(mnum)
6198 if (molnum(nct).eq.5) then
6201 gamvec(ind1)=gamp(mnum)
6203 ! Load the friction coefficients corresponding to side chains
6207 gamsc(ntyp1_molec(j),j)=1.0d0
6214 gamvec(ii)=gamsc(iabs(iti),mnum)
6216 if (surfarea) call sdarea(gamvec)
6222 innt=chain_border(1,ichain)
6223 inct=chain_border(2,ichain)
6224 !c write (iout,*) "ichain",ichain," innt",innt," inct",inct
6227 DMfric(ind)=gamvec(innt-nnt+1)/4
6228 if (iabs(itype(innt,1)).eq.10.or.mnum.gt.2) then
6229 DMfric(ind)=DMfric(ind)+gamvec(m+innt-nnt+1)
6232 DMfric(ind+1)=gamvec(m+innt-nnt+1)
6235 !c write (iout,*) "DMfric init ind",ind
6239 DMfric(ind)=gamvec(i-nnt+1)/2
6240 if (iabs(itype(i,1)).eq.10.or.mnum.gt.2) then
6241 DMfric(ind)=DMfric(ind)+gamvec(m+i-nnt+1)
6244 DMfric(ind+1)=gamvec(m+i-nnt+1)
6248 !c write (iout,*) "DMfric endloop ind",ind
6249 if (inct.gt.innt) then
6250 DMfric(ind)=gamvec(inct-1-nnt+1)/4
6252 if (iabs(itype(inct,1)).eq.10.or.mnum.gt.2) then
6253 DMfric(ind)=DMfric(ind)+gamvec(inct+m-nnt+1)
6256 DMfric(ind+1)=gamvec(inct+m-nnt+1)
6260 !c write (iout,*) "DMfric end ind",ind
6266 ind=iposd_chain(ichain)
6267 innt=chain_border(1,ichain)
6268 inct=chain_border(2,ichain)
6270 if (iabs(itype(i,1)).ne.10.and.mnum.le.2) then
6273 DU1fric(ind)=gamvec(i-nnt+1)/4
6281 ind=iposd_chain(ichain)
6282 innt=chain_border(1,ichain)
6283 inct=chain_border(2,ichain)
6285 if (iabs(itype(i,1)).ne.10.and.mnum.le.2) then
6286 DU2fric(ind)=gamvec(i-nnt+1)/4
6287 DU2fric(ind+1)=0.0d0
6296 write(iout,*)"The upper part of the five-diagonal friction matrix"
6298 write (iout,'(a,i5)') 'Chain',ichain
6299 innt=iposd_chain(ichain)
6300 inct=iposd_chain(ichain)+dimen_chain(ichain)-1
6302 if (i.lt.inct-1) then
6303 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMfric(i),DU1fric(i),&
6305 else if (i.eq.inct-1) then
6306 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMfric(i),DU1fric(i)
6308 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMfric(i)
6317 ! Zeroing out fricmat
6323 ! Load the friction coefficients corresponding to peptide groups
6328 gamvec(ind1)=gamp(mnum)
6331 if (molnum(nct).eq.5) then
6334 gamvec(ind1)=gamp(mnum)
6336 ! Load the friction coefficients corresponding to side chains
6340 gamsc(ntyp1_molec(j),j)=1.0d0
6347 gamvec(ii)=gamsc(iabs(iti),mnum)
6349 if (surfarea) call sdarea(gamvec)
6351 ! write (iout,*) "Matrix A and vector gamma"
6353 ! write (iout,'(i2,$)') i
6355 ! write (iout,'(f4.1,$)') A(i,j)
6357 ! write (iout,'(f8.3)') gamvec(i)
6361 write (iout,*) "Vector gamvec"
6363 write (iout,'(i5,f10.5)') i, gamvec(i)
6367 ! The friction matrix
6372 dtdi=dtdi+A(j,k)*A(j,i)*gamvec(j)
6379 write (iout,'(//a)') "Matrix fricmat"
6380 call matout2(dimen,dimen,nres2,nres2,fricmat)
6382 if (lang.eq.2 .or. lang.eq.3) then
6383 ! Mass-scale the friction matrix if non-direct integration will be performed
6389 Ginvfric(i,j)=Ginvfric(i,j)+ &
6390 Gsqrm(i,k)*Gsqrm(l,j)*fricmat(k,l)
6395 ! Diagonalize the friction matrix
6400 Ghalf(ind)=Ginvfric(i,j)
6403 call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
6406 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
6407 " mass-scaled friction matrix"
6408 call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
6410 ! Precompute matrices for tinker stochastic integrator
6417 mt1(i,j)=mt1(i,j)+fricvec(k,i)*gsqrm(k,j)
6418 mt2(i,j)=mt2(i,j)+fricvec(k,i)*gsqrp(k,j)
6424 else if (lang.eq.4) then
6425 ! Diagonalize the friction matrix
6430 Ghalf(ind)=fricmat(i,j)
6433 call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
6436 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
6438 call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
6440 ! Determine the number of zero eigenvalues of the friction matrix
6441 nzero=max0(dimen-dimen1,0)
6442 ! do while (fricgam(nzero+1).le.1.0d-5 .and. nzero.lt.dimen)
6445 write (iout,*) "Number of zero eigenvalues:",nzero
6450 fricmat(i,j)=fricmat(i,j) &
6451 +fricvec(i,k)*fricvec(j,k)/fricgam(k)
6456 write (iout,'(//a)') "Generalized inverse of fricmat"
6457 call matout(dimen,dimen,nres6,nres6,fricmat)
6462 if (nfgtasks.gt.1) then
6463 if (fg_rank.eq.0) then
6464 ! The matching BROADCAST for fg processors is called in ERGASTULUM
6470 call MPI_Bcast(10,1,MPI_INTEGER,king,FG_COMM,IERROR)
6472 time_Bcast=time_Bcast+MPI_Wtime()-time00
6474 time_Bcast=time_Bcast+tcpu()-time00
6476 ! print *,"Processor",myrank,
6477 ! & " BROADCAST iorder in SETUP_FRICMAT"
6480 write (iout,*) "setup_fricmat licznik"!,licznik !sp
6486 ! Scatter the friction matrix
6487 call MPI_Scatterv(fricmat(1,1),nginv_counts(0),&
6488 nginv_start(0),MPI_DOUBLE_PRECISION,fcopy(1,1),&
6489 myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
6492 time_scatter=time_scatter+MPI_Wtime()-time00
6493 time_scatter_fmat=time_scatter_fmat+MPI_Wtime()-time00
6495 time_scatter=time_scatter+tcpu()-time00
6496 time_scatter_fmat=time_scatter_fmat+tcpu()-time00
6500 do j=1,2*my_ng_count
6501 fricmat(j,i)=fcopy(i,j)
6504 ! write (iout,*) "My chunk of fricmat"
6505 ! call MATOUT2(my_ng_count,dimen,maxres2,maxres2,fcopy)
6510 end subroutine setup_fricmat
6511 !-----------------------------------------------------------------------------
6512 subroutine stochastic_force(stochforcvec)
6515 use random, only:anorm_distr
6516 ! implicit real*8 (a-h,o-z)
6517 ! include 'DIMENSIONS'
6518 use control, only: tcpu
6523 ! include 'COMMON.VAR'
6524 ! include 'COMMON.CHAIN'
6525 ! include 'COMMON.DERIV'
6526 ! include 'COMMON.GEO'
6527 ! include 'COMMON.LOCAL'
6528 ! include 'COMMON.INTERACT'
6529 ! include 'COMMON.MD'
6530 ! include 'COMMON.TIME1'
6532 ! include 'COMMON.LANGEVIN'
6534 ! include 'COMMON.LANGEVIN.lang0'
6536 ! include 'COMMON.IOUNITS'
6538 real(kind=8) :: x,sig,lowb,highb
6539 real(kind=8) :: ff(3),force(3,0:2*nres),zeta2,lowb2
6540 real(kind=8) :: highb2,sig2,forcvec(6*nres),stochforcvec(6*nres)
6541 real(kind=8) :: time00
6542 logical :: lprn = .false.
6543 integer :: i,j,ind,mnum
6545 integer ichain,innt,inct,iposc
6550 stochforc(j,i)=0.0d0
6560 ! Compute the stochastic forces acting on bodies. Store in force.
6566 force(j,i)=anorm_distr(x,sig,lowb,highb)
6574 force(j,i+nres)=anorm_distr(x,sig2,lowb2,highb2)
6578 time_fsample=time_fsample+MPI_Wtime()-time00
6580 time_fsample=time_fsample+tcpu()-time00
6585 innt=chain_border(1,ichain)
6586 inct=chain_border(2,ichain)
6587 iposc=iposd_chain(ichain)
6588 !c for debugging only
6589 !c innt=chain_border(1,1)
6590 !c inct=chain_border(2,1)
6591 !c iposc=iposd_chain(1)
6592 !c write (iout,*)"stochastic_force ichain=",ichain," innt",innt,
6593 !c & " inct",inct," iposc",iposc
6595 stochforcvec(ind+j)=0.5d0*force(j,innt)
6597 if (iabs(itype(innt,molnum(innt))).eq.10.or.molnum(innt).ge.3) then
6599 stochforcvec(ind+j)=stochforcvec(ind+j)+force(j,innt+nres)
6605 stochforcvec(ind+j)=force(j,innt+nres)
6611 stochforcvec(ind+j)=0.5d0*(force(j,i)+force(j,i-1))
6613 if (iabs(itype(i,1).eq.10).or.molnum(i).ge.3) then
6615 stochforcvec(ind+j)=stochforcvec(ind+j)+force(j,i+nres)
6621 stochforcvec(ind+j)=force(j,i+nres)
6627 stochforcvec(ind+j)=0.5d0*force(j,inct-1)
6629 if (iabs(itype(inct,molnum(inct))).eq.10.or.molnum(inct).ge.3) then
6631 stochforcvec(ind+j)=stochforcvec(ind+j)+force(j,inct+nres)
6637 stochforcvec(ind+j)=force(j,inct+nres)
6641 !c write (iout,*) "chain",ichain," ind",ind
6644 write (iout,*) "stochforcvec"
6645 write (iout,'(3f10.5)') (stochforcvec(j),j=1,ind)
6648 ! Compute the stochastic forces acting on virtual-bond vectors.
6654 stochforc(j,i)=ff(j)+0.5d0*force(j,i)
6657 ff(j)=ff(j)+force(j,i)
6659 ! if (itype(i+1,1).ne.ntyp1) then
6661 if (itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
6663 stochforc(j,i)=stochforc(j,i)+force(j,i+nres+1)
6664 ff(j)=ff(j)+force(j,i+nres+1)
6669 stochforc(j,0)=ff(j)+force(j,nnt+nres)
6673 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
6674 .and.(mnum.lt.4)) then
6675 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
6677 stochforc(j,i+nres)=force(j,i+nres)
6683 stochforcvec(j)=stochforc(j,0)
6688 stochforcvec(ind+j)=stochforc(j,i)
6694 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
6695 .and.(mnum.lt.4)) then
6696 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
6698 stochforcvec(ind+j)=stochforc(j,i+nres)
6704 write (iout,*) "stochforcvec"
6706 write(iout,'(i5,e15.5)') i,stochforcvec(i)
6708 write(iout,*) "Stochastic forces backbone"
6710 write(iout,'(i5,3e15.5)') i,(stochforc(j,i),j=1,3)
6712 write(iout,*) "Stochastic forces side chain"
6714 write(iout,'(i5,3e15.5)') &
6715 i,(stochforc(j,i+nres),j=1,3)
6723 write (iout,*) i,ind
6725 forcvec(ind+j)=force(j,i)
6730 write (iout,*) i,ind
6732 forcvec(j+ind)=force(j,i+nres)
6737 write (iout,*) "forcvec"
6741 write (iout,'(2i3,2f10.5)') i,j,force(j,i),&
6748 write (iout,'(2i3,2f10.5)') i,j,force(j,i+nres),&
6757 end subroutine stochastic_force
6758 !-----------------------------------------------------------------------------
6759 subroutine sdarea(gamvec)
6761 ! Scale the friction coefficients according to solvent accessible surface areas
6762 ! Code adapted from TINKER
6766 ! implicit real*8 (a-h,o-z)
6767 ! include 'DIMENSIONS'
6768 ! include 'COMMON.CONTROL'
6769 ! include 'COMMON.VAR'
6770 ! include 'COMMON.MD'
6772 ! include 'COMMON.LANGEVIN'
6774 ! include 'COMMON.LANGEVIN.lang0'
6776 ! include 'COMMON.CHAIN'
6777 ! include 'COMMON.DERIV'
6778 ! include 'COMMON.GEO'
6779 ! include 'COMMON.LOCAL'
6780 ! include 'COMMON.INTERACT'
6781 ! include 'COMMON.IOUNITS'
6782 ! include 'COMMON.NAMES'
6783 real(kind=8),dimension(2*nres) :: radius,gamvec !(maxres2)
6784 real(kind=8),parameter :: twosix = 1.122462048309372981d0
6785 logical :: lprn = .false.
6786 real(kind=8) :: probe,area,ratio
6787 integer :: i,j,ind,iti,mnum
6789 ! determine new friction coefficients every few SD steps
6791 ! set the atomic radii to estimates of sigma values
6793 ! print *,"Entered sdarea"
6799 ! Load peptide group radii
6802 radius(i)=pstok(mnum)
6804 ! Load side chain radii
6808 radius(i+nres)=restok(iti,mnum)
6811 ! write (iout,*) "i",i," radius",radius(i)
6814 radius(i) = radius(i) / twosix
6815 if (radius(i) .ne. 0.0d0) radius(i) = radius(i) + probe
6818 ! scale atomic friction coefficients by accessible area
6820 if (lprn) write (iout,*) &
6821 "Original gammas, surface areas, scaling factors, new gammas, ",&
6822 "std's of stochastic forces"
6825 if (radius(i).gt.0.0d0) then
6826 call surfatom (i,area,radius)
6827 ratio = dmax1(area/(4.0d0*pi*radius(i)**2),1.0d-1)
6828 if (lprn) write (iout,'(i5,3f10.5,$)') &
6829 i,gamvec(ind+1),area,ratio
6832 gamvec(ind) = ratio * gamvec(ind)
6835 stdforcp(i)=stdfp(mnum)*dsqrt(gamvec(ind))
6836 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcp(i)
6840 if (radius(i+nres).gt.0.0d0) then
6841 call surfatom (i+nres,area,radius)
6842 ratio = dmax1(area/(4.0d0*pi*radius(i+nres)**2),1.0d-1)
6843 if (lprn) write (iout,'(i5,3f10.5,$)') &
6844 i,gamvec(ind+1),area,ratio
6847 gamvec(ind) = ratio * gamvec(ind)
6850 stdforcsc(i)=stdfsc(itype(i,mnum),mnum)*dsqrt(gamvec(ind))
6851 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcsc(i)
6856 end subroutine sdarea
6857 !-----------------------------------------------------------------------------
6859 !-----------------------------------------------------------------------------
6862 ! ###################################################
6863 ! ## COPYRIGHT (C) 1996 by Jay William Ponder ##
6864 ! ## All Rights Reserved ##
6865 ! ###################################################
6867 ! ################################################################
6869 ! ## subroutine surfatom -- exposed surface area of an atom ##
6871 ! ################################################################
6874 ! "surfatom" performs an analytical computation of the surface
6875 ! area of a specified atom; a simplified version of "surface"
6877 ! literature references:
6879 ! T. J. Richmond, "Solvent Accessible Surface Area and
6880 ! Excluded Volume in Proteins", Journal of Molecular Biology,
6883 ! L. Wesson and D. Eisenberg, "Atomic Solvation Parameters
6884 ! Applied to Molecular Dynamics of Proteins in Solution",
6885 ! Protein Science, 1, 227-235 (1992)
6887 ! variables and parameters:
6889 ! ir number of atom for which area is desired
6890 ! area accessible surface area of the atom
6891 ! radius radii of each of the individual atoms
6894 subroutine surfatom(ir,area,radius)
6896 ! implicit real*8 (a-h,o-z)
6897 ! include 'DIMENSIONS'
6899 ! include 'COMMON.GEO'
6900 ! include 'COMMON.IOUNITS'
6902 integer :: nsup,nstart_sup
6903 ! double precision c,dc,dc_old,d_c_work,xloc,xrot,dc_norm
6904 ! common /chain/ c(3,maxres2+2),dc(3,0:maxres2),dc_old(3,0:maxres2),
6905 ! & xloc(3,maxres),xrot(3,maxres),dc_norm(3,0:maxres2),
6906 ! & dc_work(MAXRES6),nres,nres0
6907 integer,parameter :: maxarc=300
6911 integer :: mi,ni,narc
6912 integer :: key(maxarc)
6913 integer :: intag(maxarc)
6914 integer :: intag1(maxarc)
6915 real(kind=8) :: area,arcsum
6916 real(kind=8) :: arclen,exang
6917 real(kind=8) :: delta,delta2
6918 real(kind=8) :: eps,rmove
6919 real(kind=8) :: xr,yr,zr
6920 real(kind=8) :: rr,rrsq
6921 real(kind=8) :: rplus,rminus
6922 real(kind=8) :: axx,axy,axz
6923 real(kind=8) :: ayx,ayy
6924 real(kind=8) :: azx,azy,azz
6925 real(kind=8) :: uxj,uyj,uzj
6926 real(kind=8) :: tx,ty,tz
6927 real(kind=8) :: txb,tyb,td
6928 real(kind=8) :: tr2,tr,txr,tyr
6929 real(kind=8) :: tk1,tk2
6930 real(kind=8) :: thec,the,t,tb
6931 real(kind=8) :: txk,tyk,tzk
6932 real(kind=8) :: t1,ti,tf,tt
6933 real(kind=8) :: txj,tyj,tzj
6934 real(kind=8) :: ccsq,cc,xysq
6935 real(kind=8) :: bsqk,bk,cosine
6936 real(kind=8) :: dsqj,gi,pix2
6937 real(kind=8) :: therk,dk,gk
6938 real(kind=8) :: risqk,rik
6939 real(kind=8) :: radius(2*nres) !(maxatm) (maxatm=maxres2)
6940 real(kind=8) :: ri(maxarc),risq(maxarc)
6941 real(kind=8) :: ux(maxarc),uy(maxarc),uz(maxarc)
6942 real(kind=8) :: xc(maxarc),yc(maxarc),zc(maxarc)
6943 real(kind=8) :: xc1(maxarc),yc1(maxarc),zc1(maxarc)
6944 real(kind=8) :: dsq(maxarc),bsq(maxarc)
6945 real(kind=8) :: dsq1(maxarc),bsq1(maxarc)
6946 real(kind=8) :: arci(maxarc),arcf(maxarc)
6947 real(kind=8) :: ex(maxarc),lt(maxarc),gr(maxarc)
6948 real(kind=8) :: b(maxarc),b1(maxarc),bg(maxarc)
6949 real(kind=8) :: kent(maxarc),kout(maxarc)
6950 real(kind=8) :: ther(maxarc)
6951 logical :: moved,top
6952 logical :: omit(maxarc)
6955 ! maxatm = 2*nres !maxres2 maxres2=2*maxres
6956 ! maxlight = 8*maxatm
6959 ! maxtors = 4*maxatm
6961 ! zero out the surface area for the sphere of interest
6964 ! write (2,*) "ir",ir," radius",radius(ir)
6965 if (radius(ir) .eq. 0.0d0) return
6967 ! set the overlap significance and connectivity shift
6971 delta2 = delta * delta
6976 ! store coordinates and radius of the sphere of interest
6984 ! initialize values of some counters and summations
6993 ! test each sphere to see if it overlaps the sphere of interest
6996 if (i.eq.ir .or. radius(i).eq.0.0d0) goto 30
6997 rplus = rr + radius(i)
6999 if (abs(tx) .ge. rplus) goto 30
7001 if (abs(ty) .ge. rplus) goto 30
7003 if (abs(tz) .ge. rplus) goto 30
7005 ! check for sphere overlap by testing distance against radii
7007 xysq = tx*tx + ty*ty
7008 if (xysq .lt. delta2) then
7015 if (rplus-cc .le. delta) goto 30
7016 rminus = rr - radius(i)
7018 ! check to see if sphere of interest is completely buried
7020 if (cc-abs(rminus) .le. delta) then
7021 if (rminus .le. 0.0d0) goto 170
7025 ! check for too many overlaps with sphere of interest
7027 if (io .ge. maxarc) then
7029 20 format (/,' SURFATOM -- Increase the Value of MAXARC')
7033 ! get overlap between current sphere and sphere of interest
7042 gr(io) = (ccsq+rplus*rminus) / (2.0d0*rr*b1(io))
7048 ! case where no other spheres overlap the sphere of interest
7051 area = 4.0d0 * pi * rrsq
7055 ! case where only one sphere overlaps the sphere of interest
7058 area = pix2 * (1.0d0 + gr(1))
7059 area = mod(area,4.0d0*pi) * rrsq
7063 ! case where many spheres intersect the sphere of interest;
7064 ! sort the intersecting spheres by their degree of overlap
7066 call sort2 (io,gr,key)
7069 intag(i) = intag1(k)
7078 ! get radius of each overlap circle on surface of the sphere
7083 risq(i) = rrsq - gi*gi
7084 ri(i) = sqrt(risq(i))
7085 ther(i) = 0.5d0*pi - asin(min(1.0d0,max(-1.0d0,gr(i))))
7088 ! find boundary of inaccessible area on sphere of interest
7091 if (.not. omit(k)) then
7098 ! check to see if J circle is intersecting K circle;
7099 ! get distance between circle centers and sum of radii
7102 if (omit(j)) goto 60
7103 cc = (txk*xc(j)+tyk*yc(j)+tzk*zc(j))/(bk*b(j))
7104 cc = acos(min(1.0d0,max(-1.0d0,cc)))
7105 td = therk + ther(j)
7107 ! check to see if circles enclose separate regions
7109 if (cc .ge. td) goto 60
7111 ! check for circle J completely inside circle K
7113 if (cc+ther(j) .lt. therk) goto 40
7115 ! check for circles that are essentially parallel
7117 if (cc .gt. delta) goto 50
7122 ! check to see if sphere of interest is completely buried
7125 if (pix2-cc .le. td) goto 170
7131 ! find T value of circle intersections
7134 if (omit(k)) goto 110
7149 ! rotation matrix elements
7161 if (.not. omit(j)) then
7166 ! rotate spheres so K vector colinear with z-axis
7168 uxj = txj*axx + tyj*axy - tzj*axz
7169 uyj = tyj*ayy - txj*ayx
7170 uzj = txj*azx + tyj*azy + tzj*azz
7171 cosine = min(1.0d0,max(-1.0d0,uzj/b(j)))
7172 if (acos(cosine) .lt. therk+ther(j)) then
7173 dsqj = uxj*uxj + uyj*uyj
7178 tr2 = risqk*dsqj - tb*tb
7184 ! get T values of intersection for K circle
7187 tb = min(1.0d0,max(-1.0d0,tb))
7189 if (tyb-txr .lt. 0.0d0) tk1 = pix2 - tk1
7191 tb = min(1.0d0,max(-1.0d0,tb))
7193 if (tyb+txr .lt. 0.0d0) tk2 = pix2 - tk2
7194 thec = (rrsq*uzj-gk*bg(j)) / (rik*ri(j)*b(j))
7195 if (abs(thec) .lt. 1.0d0) then
7197 else if (thec .ge. 1.0d0) then
7199 else if (thec .le. -1.0d0) then
7203 ! see if "tk1" is entry or exit point; check t=0 point;
7204 ! "ti" is exit point, "tf" is entry point
7206 cosine = min(1.0d0,max(-1.0d0, &
7207 (uzj*gk-uxj*rik)/(b(j)*rr)))
7208 if ((acos(cosine)-ther(j))*(tk2-tk1) .le. 0.0d0) then
7216 if (narc .ge. maxarc) then
7218 70 format (/,' SURFATOM -- Increase the Value',&
7222 if (tf .le. ti) then
7243 ! special case; K circle without intersections
7245 if (narc .le. 0) goto 90
7247 ! general case; sum up arclength and set connectivity code
7249 call sort2 (narc,arci,key)
7254 if (narc .gt. 1) then
7257 if (t .lt. arci(j)) then
7258 arcsum = arcsum + arci(j) - t
7259 exang = exang + ex(ni)
7261 if (jb .ge. maxarc) then
7263 80 format (/,' SURFATOM -- Increase the Value',&
7268 kent(jb) = maxarc*i + k
7270 kout(jb) = maxarc*k + i
7279 arcsum = arcsum + pix2 - t
7281 exang = exang + ex(ni)
7284 kent(jb) = maxarc*i + k
7286 kout(jb) = maxarc*k + i
7293 arclen = arclen + gr(k)*arcsum
7296 if (arclen .eq. 0.0d0) goto 170
7297 if (jb .eq. 0) goto 150
7299 ! find number of independent boundaries and check connectivity
7303 if (kout(k) .ne. 0) then
7310 if (m .eq. kent(ii)) then
7313 if (j .eq. jb) goto 150
7325 ! attempt to fix connectivity error by moving atom slightly
7329 140 format (/,' SURFATOM -- Connectivity Error at Atom',i6)
7338 ! compute the exposed surface area for the sphere of interest
7341 area = ib*pix2 + exang + arclen
7342 area = mod(area,4.0d0*pi) * rrsq
7344 ! attempt to fix negative area by moving atom slightly
7346 if (area .lt. 0.0d0) then
7349 160 format (/,' SURFATOM -- Negative Area at Atom',i6)
7360 end subroutine surfatom
7361 !----------------------------------------------------------------
7362 !----------------------------------------------------------------
7363 subroutine alloc_MD_arrays
7364 !EL Allocation of arrays used by MD module
7366 integer :: nres2,nres6
7369 !----------------------
7373 allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
7374 allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
7375 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
7376 allocate(fricvec(nres2,nres2))
7377 allocate(pfric_mat(nres2,nres2),vfric_mat(nres2,nres2))
7378 allocate(afric_mat(nres2,nres2),prand_mat(nres2,nres2))
7379 allocate(vrand_mat1(nres2,nres2),vrand_mat2(nres2,nres2)) !(MAXRES2,MAXRES2)
7380 allocate(pfric0_mat(nres2,nres2,0:maxflag_stoch))
7381 allocate(afric0_mat(nres2,nres2,0:maxflag_stoch))
7382 allocate(vfric0_mat(nres2,nres2,0:maxflag_stoch))
7383 allocate(prand0_mat(nres2,nres2,0:maxflag_stoch))
7384 allocate(vrand0_mat1(nres2,nres2,0:maxflag_stoch))
7385 allocate(vrand0_mat2(nres2,nres2,0:maxflag_stoch)) !(MAXRES2,MAXRES2,0:maxflag_stoch)
7386 allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
7388 allocate(mt1(nres2,nres2),mt2(nres2,nres2),mt3(nres2,nres2)) !(maxres2,maxres2)
7389 !----------------------
7391 ! commom.langevin.lang0
7393 allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
7395 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
7396 allocate(fricvec(nres2,nres2)) !(MAXRES2,MAXRES2)
7398 allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
7399 allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
7402 if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2))
7404 !----------------------
7405 ! commom.hairpin in CSA module
7406 !----------------------
7407 ! common.mce in MCM_MD module
7408 !----------------------
7410 ! common /mdgrad/ in module.energy
7411 ! common /back_constr/ in module.energy
7412 ! common /qmeas/ in module.energy
7415 allocate(potEcomp(0:n_ene+4)) !(0:n_ene+4)
7417 allocate(d_t(3,0:nres2),d_a(3,0:nres2),d_t_old(3,0:nres2)) !(3,0:MAXRES2)
7418 allocate(d_a_work(nres6)) !(6*MAXRES)
7420 allocate(DM(nres2),DU1(nres2),DU2(nres2))
7421 allocate(DMorig(nres2),DU1orig(nres2),DU2orig(nres2))
7423 allocate(Gvec(1300,1300))!maximum number of protein in one chain
7426 write (iout,*) "Before A Allocation",nfgtasks-1
7428 allocate(Gmat(nres2,nres2),A(nres2,nres2))
7429 if(.not.allocated(Ginv)) allocate(Ginv(nres2,nres2)) !in control: ergastulum
7430 allocate(Gsqrp(nres2,nres2),Gsqrm(nres2,nres2),Gvec(nres2,nres2)) !(maxres2,maxres2)
7432 allocate(Geigen(nres2)) !(maxres2)
7433 if(.not.allocated(vtot)) allocate(vtot(nres2)) !(maxres2)
7434 ! common /inertia/ in io_conf: parmread
7435 ! real(kind=8),dimension(:),allocatable :: ISC,msc !(ntyp+1)
7436 ! common /langevin/in io read_MDpar
7437 ! real(kind=8),dimension(:),allocatable :: gamsc !(ntyp1)
7438 ! real(kind=8),dimension(:),allocatable :: stdfsc !(ntyp)
7439 ! in io_conf: parmread
7440 ! real(kind=8),dimension(:),allocatable :: restok !(ntyp+1)
7441 ! common /mdpmpi/ in control: ergastulum
7442 if(.not.allocated(ng_start)) allocate(ng_start(0:nfgtasks-1))
7443 if(.not.allocated(ng_counts)) allocate(ng_counts(0:nfgtasks-1))
7444 if(.not.allocated(nginv_counts)) allocate(nginv_counts(0:nfgtasks-1)) !(0:MaxProcs-1)
7445 if(.not.allocated(nginv_start)) allocate(nginv_start(0:nfgtasks)) !(0:MaxProcs)
7446 !----------------------
7447 ! common.muca in read_muca
7448 ! common /double_muca/
7449 ! real(kind=8) :: elow,ehigh,factor,hbin,factor_min
7450 ! real(kind=8),dimension(:),allocatable :: emuca,nemuca,&
7451 ! nemuca2,hist !(4*maxres)
7452 ! common /integer_muca/
7453 ! integer :: nmuca,imtime,muca_smooth
7455 ! real(kind=8),dimension(:),allocatable :: elowi,ehighi !(maxprocs)
7456 !----------------------
7458 ! common /mdgrad/ in module.energy
7459 ! common /back_constr/ in module.energy
7460 ! common /qmeas/ in module.energy
7464 allocate(d_t_work(nres6),d_t_work_new(nres6),d_af_work(nres6))
7465 allocate(d_as_work(nres6),kinetic_force(nres6)) !(MAXRES6)
7466 allocate(d_t_new(3,0:nres2),d_a_old(3,0:nres2),d_a_short(3,0:nres2)) !,d_a !(3,0:MAXRES2)
7467 allocate(stdforcp(nres),stdforcsc(nres)) !(MAXRES)
7468 !----------------------
7470 allocate(D_ban(nres6)) !(MAXRES6) maxres6=6*maxres
7471 ! common /stochcalc/ stochforcvec
7472 allocate(stochforcvec(nres6)) !(MAXRES6) maxres6=6*maxres
7473 !----------------------
7475 end subroutine alloc_MD_arrays
7476 !-----------------------------------------------------------------------------
7477 !-----------------------------------------------------------------------------