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,mnum).eq.ntyp1_molec(mnum)) then
707 if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
714 v(j)=incr(j)+d_t(j,nres+i)
717 ! if (mnum.ne.5) then
718 ! write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3)
719 ! write (iout,*) "i",i," msc",msc(iti,mnum)," v",(v(j),j=1,3)
720 KEt_sc=KEt_sc+msc(iti,mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
721 vtot(i+nres)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
724 incr(j)=incr(j)+d_t(j,i)
728 ! write(iout,*) 'KEt_sc', KEt_sc
729 ! The part due to stretching and rotation of the peptide groups
732 if (itype(i,mnum).ne.ntyp1_molec(mnum)&
733 .and.itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
734 if (mnum.eq.5) Ip(mnum)=0.0
735 ! write (iout,*) "i",i
736 ! write (iout,*) "i",i," mag1",mag1," mag2",mag2
740 !c write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3)
741 KEr_p=KEr_p+Ip(mnum)*(incr(1)*incr(1)+incr(2)*incr(2) &
746 !c write(iout,*) 'KEr_p', KEr_p
747 !c The rotational part of the side chain virtual bond
750 iti=iabs(itype(i,mnum))
751 if (itype(i,1).ne.10.and.itype(i,mnum).ne.ntyp1_molec(mnum)&
754 incr(j)=d_t(j,nres+i)
756 ! write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3)
757 KEr_sc=KEr_sc+Isc(iti,mnum)*(incr(1)*incr(1)+incr(2)*incr(2)+&
761 !c The total kinetic energy
763 ! write(iout,*) ' KEt_p',KEt_p,' KEt_sc',KEt_sc,' KEr_p',KEr_p, &
765 KE_total=0.5d0*(KEt_p+KEt_sc+0.25d0*KEr_p+KEr_sc)
766 !c write (iout,*) "KE_total",KE_total
768 end subroutine kinetic
771 !-----------------------------------------------------------------------------
772 subroutine kinetic(KE_total)
773 !----------------------------------------------------------------
774 ! This subroutine calculates the total kinetic energy of the chain
775 !-----------------------------------------------------------------
777 ! implicit real*8 (a-h,o-z)
778 ! include 'DIMENSIONS'
779 ! include 'COMMON.VAR'
780 ! include 'COMMON.CHAIN'
781 ! include 'COMMON.DERIV'
782 ! include 'COMMON.GEO'
783 ! include 'COMMON.LOCAL'
784 ! include 'COMMON.INTERACT'
785 ! include 'COMMON.MD'
786 ! include 'COMMON.IOUNITS'
787 real(kind=8) :: KE_total,mscab
789 integer :: i,j,k,iti,mnum,term
790 real(kind=8) :: KEt_p,KEt_sc,KEr_p,KEr_sc,incr(3),&
793 write (iout,*) "Velocities, kietic"
795 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
796 (d_t(j,i+nres),j=1,3)
801 ! write (iout,*) "ISC",(isc(itype(i,1)),i=1,nres)
802 ! The translational part for peptide virtual bonds
807 ! if (molnum(nct).gt.3) term=nct
810 if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
811 ! write (iout,*) "Kinetic trp:",i,(incr(j),j=1,3),mp(mnum)
814 v(j)=incr(j)+0.5d0*d_t(j,i)
818 v(j)=incr(j)+0.5d0*d_t(j,i)
821 vtot(i)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
822 KEt_p=KEt_p+mp(mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
824 incr(j)=incr(j)+d_t(j,i)
827 ! write(iout,*) 'KEt_p', KEt_p
828 ! The translational part for the side chain virtual bond
829 ! Only now we can initialize incr with zeros. It must be equal
830 ! to the velocities of the first Calpha.
836 iti=iabs(itype(i,mnum))
837 ! if (mnum.ge.4) then
842 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
843 if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
844 .or.(mnum.ge.4)) then
850 v(j)=incr(j)+d_t(j,nres+i)
853 ! write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3)
854 ! write (iout,*) "i",i," msc",msc(iti,mnum)," v",(v(j),j=1,3)
855 KEt_sc=KEt_sc+mscab*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
856 vtot(i+nres)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
858 incr(j)=incr(j)+d_t(j,i)
862 ! write(iout,*) 'KEt_sc', KEt_sc
863 ! The part due to stretching and rotation of the peptide groups
867 ! write (iout,*) "i",i
868 ! write (iout,*) "i",i," mag1",mag1," mag2",mag2
872 ! write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3)
873 KEr_p=KEr_p+Ip(mnum)*(incr(1)*incr(1)+incr(2)*incr(2) &
877 ! write(iout,*) 'KEr_p', KEr_p
878 ! The rotational part of the side chain virtual bond
882 iti=iabs(itype(i,mnum))
883 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
884 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
885 .and.(mnum.lt.4)) then
887 incr(j)=d_t(j,nres+i)
889 ! write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3)
890 KEr_sc=KEr_sc+Isc(iti,mnum)*(incr(1)*incr(1)+incr(2)*incr(2)+ &
894 ! The total kinetic energy
896 ! write(iout,*) 'KEr_sc', KEr_sc
897 KE_total=0.5d0*(KEt_p+KEt_sc+0.25d0*KEr_p+KEr_sc)
898 ! write (iout,*) "KE_total",KE_total
900 end subroutine kinetic
901 !-----------------------------------------------------------------------------
903 subroutine kinetic_CASC(KE_total)
904 !c----------------------------------------------------------------
905 !c Compute the kinetic energy of the system using the Calpha-SC
907 !c-----------------------------------------------------------------
909 double precision KE_total
911 integer i,j,k,iti,ichain,innt,inct,mnum
912 double precision KEt_p,KEt_sc,KEr_p,KEr_sc,incr(3),&
919 !c write (iout,*) "ISC",(isc(itype(i)),i=1,nres)
920 !c The translational part for peptide virtual bonds
923 innt=chain_border(1,ichain)
924 inct=chain_border(2,ichain)
925 !c write (iout,*) "Kinetic_CASC chain",ichain," innt",innt,
930 if (mnum.eq.5) mp(mnum)=0.0d0
931 ! if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
932 !c write (iout,*) i,(d_t(j,i),j=1,3),(d_t(j,i+1),j=1,3)
934 v(j)=0.5d0*(d_t(j,i)+d_t(j,i+1))
936 !c write (iout,*) "Kinetic trp i",i," v",(v(j),j=1,3)
937 KEt_p=KEt_p+mp(mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
939 !c write(iout,*) 'KEt_p', KEt_p
940 !c The translational part for the side chain virtual bond
941 !c Only now we can initialize incr with zeros. It must be equal
942 !c to the velocities of the first Calpha.
948 iti=iabs(itype(i,mnum))
950 if (itype(i,1).eq.10.or.mnum.ge.3.or. itype(i,mnum).eq.ntyp1_molec(mnum)) then
951 !c write (iout,*) i,iti,(d_t(j,i),j=1,3)
956 !c write (iout,*) i,iti,(d_t(j,nres+i),j=1,3)
961 !c write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3)
962 !c write (iout,*) "i",i," msc",msc(iti)," v",(v(j),j=1,3)
963 KEt_sc=KEt_sc+msc(iti,mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
966 !c write(iout,*) 'KEt_sc', KEt_sc
967 !c The part due to stretching and rotation of the peptide groups
971 incr(j)=d_t(j,i+1)-d_t(j,i)
973 if (mnum.eq.5) Ip(mnum)=0.0d0
974 !c write (iout,*) i,(incr(j),j=1,3)
975 !c write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3)
976 KEr_p=KEr_p+Ip(mnum)*(incr(1)*incr(1)+incr(2)*incr(2)&
980 !c write(iout,*) 'KEr_p', KEr_p
981 !c The rotational part of the side chain virtual bond
984 ! iti=iabs(itype(i,mnum))
988 iti=iabs(itype(i,mnum))
991 ! if (iti.ne.10.and.mnum.lt.3) then
992 if (itype(i,1).ne.10.and.mnum.lt.3.and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
994 incr(j)=d_t(j,nres+i)-d_t(j,i)
996 !c write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3)
997 KEr_sc=KEr_sc+Isc(iti,mnum)*(incr(1)*incr(1)+incr(2)*incr(2)+&
1003 !c The total kinetic energy
1005 !c write(iout,*) ' KEt_p',KEt_p,' KEt_sc',KEt_sc,' KEr_p',KEr_p,
1006 !c & ' KEr_sc', KEr_sc
1007 KE_total=0.5d0*(KEt_p+KEt_sc+0.25d0*KEr_p+KEr_sc)
1008 !c write (iout,*) "KE_total",KE_tota
1010 write (iout,*) "Need to compile with -DFIVEDIAG to use this sub!"
1014 end subroutine kinetic_CASC
1017 !-----------------------------------------------------------------------------
1019 !------------------------------------------------
1020 ! The driver for molecular dynamics subroutines
1021 !------------------------------------------------
1024 use control, only:tcpu,ovrtim
1025 ! use io_comm, only:ilen
1027 use compare, only:secondary2,hairpin
1028 use io, only:cartout,statout
1029 ! implicit real*8 (a-h,o-z)
1030 ! include 'DIMENSIONS'
1033 integer :: IERROR,ERRCODE
1035 ! include 'COMMON.SETUP'
1036 ! include 'COMMON.CONTROL'
1037 ! include 'COMMON.VAR'
1038 ! include 'COMMON.MD'
1040 ! include 'COMMON.LANGEVIN'
1042 ! include 'COMMON.LANGEVIN.lang0'
1044 ! include 'COMMON.CHAIN'
1045 ! include 'COMMON.DERIV'
1046 ! include 'COMMON.GEO'
1047 ! include 'COMMON.LOCAL'
1048 ! include 'COMMON.INTERACT'
1049 ! include 'COMMON.IOUNITS'
1050 ! include 'COMMON.NAMES'
1051 ! include 'COMMON.TIME1'
1052 ! include 'COMMON.HAIRPIN'
1053 real(kind=8),dimension(3) :: L,vcm,boxx
1055 real(kind=8),dimension(6*nres) :: v_work,v_transf !(maxres6) maxres6=6*maxres
1057 integer :: rstcount !ilen,
1059 character(len=50) :: tytul
1060 !el common /gucio/ cm
1061 integer :: i,j,nharp
1062 integer,dimension(4,nres) :: iharp !(4,nres/3)(4,maxres/3)
1065 real(kind=8) :: tt0,scalfac
1066 integer :: nres2,itime
1075 print *,"MY tmpdir",tmpdir,ilen(tmpdir)
1076 if (ilen(tmpdir).gt.0) &
1077 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_" &
1078 //liczba(:ilen(liczba))//'.rst')
1080 if (ilen(tmpdir).gt.0) &
1081 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"//'.rst')
1088 write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
1094 print *,"just befor setup matix",nres
1095 ! Determine the inverse of the inertia matrix.
1096 call setup_MD_matrices
1098 print *,"AFTER SETUP MATRICES"
1100 print *,"AFTER INIT MD"
1103 t_MDsetup = MPI_Wtime()-tt0
1105 t_MDsetup = tcpu()-tt0
1108 ! Entering the MD loop
1114 if (lang.eq.2 .or. lang.eq.3) then
1118 call sd_verlet_p_setup
1120 call sd_verlet_ciccotti_setup
1124 pfric0_mat(i,j,0)=pfric_mat(i,j)
1125 afric0_mat(i,j,0)=afric_mat(i,j)
1126 vfric0_mat(i,j,0)=vfric_mat(i,j)
1127 prand0_mat(i,j,0)=prand_mat(i,j)
1128 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
1129 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
1132 flag_stoch(0)=.true.
1133 do i=1,maxflag_stoch
1134 flag_stoch(i)=.false.
1138 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
1140 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1144 else if (lang.eq.1 .or. lang.eq.4) then
1145 print *,"before setup_fricmat"
1147 print *,"after setup_fricmat"
1150 t_langsetup=MPI_Wtime()-tt0
1153 t_langsetup=tcpu()-tt0
1156 do itime=1,n_timestep
1157 if (large) print *,itime,ntwe
1159 if (large.and. mod(itime,ntwe).eq.0) &
1160 write (iout,*) "itime",itime
1162 if (lang.gt.0 .and. surfarea .and. &
1163 mod(itime,reset_fricmat).eq.0) then
1164 if (lang.eq.2 .or. lang.eq.3) then
1168 call sd_verlet_p_setup
1170 call sd_verlet_ciccotti_setup
1174 pfric0_mat(i,j,0)=pfric_mat(i,j)
1175 afric0_mat(i,j,0)=afric_mat(i,j)
1176 vfric0_mat(i,j,0)=vfric_mat(i,j)
1177 prand0_mat(i,j,0)=prand_mat(i,j)
1178 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
1179 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
1182 flag_stoch(0)=.true.
1183 do i=1,maxflag_stoch
1184 flag_stoch(i)=.false.
1187 else if (lang.eq.1 .or. lang.eq.4) then
1188 print *,"before setup_fricmat"
1190 print *,"after setup_fricmat"
1192 write (iout,'(a,i10)') &
1193 "Friction matrix reset based on surface area, itime",itime
1195 if (reset_vel .and. tbf .and. lang.eq.0 &
1196 .and. mod(itime,count_reset_vel).eq.0) then
1199 if (molnum(i).eq.5) then
1200 call to_box(c(1,i),c(2,i),c(3,i))
1202 if (c(j,i).le.0) c(j,i)=c(j,i)+boxx(j)
1205 ! write(iout,*) "COORD",c(1,i),c(2,i),c(3,i)
1209 write(iout,'(a,f20.2)') &
1210 "Velocities reset to random values, time",totT
1213 d_t_old(j,i)=d_t(j,i)
1217 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
1221 d_t(j,0)=d_t(j,0)-vcm(j)
1224 kinetic_T=2.0d0/(dimen3*Rb)*EK
1225 scalfac=dsqrt(T_bath/kinetic_T)
1226 write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
1229 d_t_old(j,i)=scalfac*d_t(j,i)
1235 ! Time-reversible RESPA algorithm
1236 ! (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
1237 call RESPA_step(itime)
1239 ! Variable time step algorithm.
1240 if (large) print *,"before verlet_step"
1241 call velverlet_step(itime)
1242 if (large) print *,"after verlet_step"
1246 call brown_step(itime)
1248 print *,"Brown dynamics not here!"
1250 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1257 if (mod(itime,ntwe).eq.0) then
1261 ! call check_ecartint
1271 v_work(ind)=d_t(j,i)
1276 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum).and.mnum.lt.4) then
1279 v_work(ind)=d_t(j,i+nres)
1284 write (66,'(80f10.5)') &
1285 ((d_t(j,i),j=1,3),i=0,nres-1),((d_t(j,i+nres),j=1,3),i=1,nres)
1289 v_transf(i)=v_transf(i)+gvec(j,i)*v_work(j)
1291 v_transf(i)= v_transf(i)*dsqrt(geigen(i))
1293 write (67,'(80f10.5)') (v_transf(i),i=1,ind)
1296 if (mod(itime,ntwx).eq.0) then
1298 call enerprint(potEcomp)
1300 write (tytul,'("time",f8.2)') totT
1302 write(iout,*) "before hairpin"
1303 call hairpin(.true.,nharp,iharp)
1304 write(iout,*) "before secondary"
1305 call secondary2(.true.)
1306 write(iout,*) "before pdbout"
1307 call pdbout(potE,tytul,ipdb)
1308 ! call enerprint(potEcomp)
1313 write(iout,*) "starting fodstep"
1314 call fodstep(nfodstep)
1315 write(iout,*) "after fodstep"
1318 call hairpin(.true.,nharp,iharp)
1319 call secondary2(.true.)
1320 call pdbout(potE,tytul,ipdb)
1327 if (rstcount.eq.1000.or.itime.eq.n_timestep) then
1328 open(irest2,file=rest2name,status='unknown')
1329 write(irest2,*) totT,EK,potE,totE,t_bath
1331 ! AL 4/17/17: Now writing d_t(0,:) too
1333 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
1335 ! AL 4/17/17: Now writing d_c(0,:) too
1337 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
1345 t_MD=MPI_Wtime()-tt0
1349 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))') &
1351 'MD calculations setup:',t_MDsetup,&
1352 'Energy & gradient evaluation:',t_enegrad,&
1353 'Stochastic MD setup:',t_langsetup,&
1354 'Stochastic MD step setup:',t_sdsetup,&
1356 write (iout,'(/28(1h=),a25,27(1h=))') &
1357 ' End of MD calculation '
1359 write (iout,*) "time for etotal",t_etotal," elong",t_elong,&
1361 write (iout,*) "time_fric",time_fric," time_stoch",time_stoch,&
1362 " time_fricmatmult",time_fricmatmult," time_fsample ",&
1367 !-----------------------------------------------------------------------------
1368 subroutine velverlet_step(itime)
1369 !-------------------------------------------------------------------------------
1370 ! Perform a single velocity Verlet step; the time step can be rescaled if
1371 ! increments in accelerations exceed the threshold
1372 !-------------------------------------------------------------------------------
1373 ! implicit real*8 (a-h,o-z)
1374 ! include 'DIMENSIONS'
1376 use control, only:tcpu
1378 use minimm, only:minim_dc
1381 integer :: ierror,ierrcode
1382 real(kind=8) :: errcode
1384 ! include 'COMMON.SETUP'
1385 ! include 'COMMON.VAR'
1386 ! include 'COMMON.MD'
1388 ! include 'COMMON.LANGEVIN'
1390 ! include 'COMMON.LANGEVIN.lang0'
1392 ! include 'COMMON.CHAIN'
1393 ! include 'COMMON.DERIV'
1394 ! include 'COMMON.GEO'
1395 ! include 'COMMON.LOCAL'
1396 ! include 'COMMON.INTERACT'
1397 ! include 'COMMON.IOUNITS'
1398 ! include 'COMMON.NAMES'
1399 ! include 'COMMON.TIME1'
1400 ! include 'COMMON.MUCA'
1401 real(kind=8),dimension(3) :: vcm,incr
1402 real(kind=8),dimension(3) :: L
1403 integer :: count,rstcount !ilen,
1405 character(len=50) :: tytul
1406 integer :: maxcount_scale = 30
1407 !el common /gucio/ cm
1408 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1409 !el common /stochcalc/ stochforcvec
1410 integer :: icount_scale,itime_scal,i,j,ifac_time,iretcode,itime
1412 real(kind=8) :: epdrift,tt0,fac_time
1414 if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres)) !(MAXRES6) maxres6=6*maxres
1420 if (large) print *,"after sddir_precalc"
1421 else if (lang.eq.2 .or. lang.eq.3) then
1423 call stochastic_force(stochforcvec)
1426 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
1428 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1435 icount_scale=icount_scale+1
1436 ! write(iout,*) "icount_scale",icount_scale,scalel
1437 if (icount_scale.gt.maxcount_scale) then
1439 "ERROR: too many attempts at scaling down the time step. ",&
1440 "amax=",amax,"epdrift=",epdrift,&
1441 "damax=",damax,"edriftmax=",edriftmax,&
1445 call MPI_Abort(MPI_COMM_WORLD,IERROR,IERRCODE)
1449 ! First step of the velocity Verlet algorithm
1454 else if (lang.eq.3) then
1456 call sd_verlet1_ciccotti
1458 else if (lang.eq.1) then
1463 ! Build the chain from the newly calculated coordinates
1464 call chainbuild_cart
1465 if (rattle) call rattle1
1467 if (large) then !.and. mod(itime,ntwe).eq.0) then
1468 write (iout,*) "Cartesian and internal coordinates: step 1"
1473 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),&
1474 (dc(j,i+nres),j=1,3)
1476 write (iout,*) "Accelerations"
1478 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1479 (d_a(j,i+nres),j=1,3)
1481 write (iout,*) "Velocities, step 1"
1483 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1484 (d_t(j,i+nres),j=1,3)
1493 ! Calculate energy and forces
1495 call etotal(potEcomp)
1496 ! AL 4/17/17: Reduce the steps if NaNs occurred.
1497 if (potEcomp(0).gt.0.99e18 .or. isnan(potEcomp(0)).gt.0) then
1498 call enerprint(potEcomp)
1500 if (icount_scale.gt.15) then
1501 write (iout,*) "Tu jest problem",potEcomp(0),d_time
1502 ! call gen_rand_conf(1,*335)
1503 ! call minim_dc(potEcomp(0),iretcode,100)
1506 ! call etotal(potEcomp)
1507 ! write(iout,*) "needed to repara,",potEcomp
1510 ! 335 write(iout,*) "Failed genrand"
1514 if (large.and. mod(itime,ntwe).eq.0) &
1515 call enerprint(potEcomp)
1518 t_etotal=t_etotal+MPI_Wtime()-tt0
1520 t_etotal=t_etotal+tcpu()-tt0
1523 potE=potEcomp(0)-potEcomp(51)
1525 ! Get the new accelerations
1528 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1530 t_enegrad=t_enegrad+tcpu()-tt0
1532 ! Determine maximum acceleration and scale down the timestep if needed
1534 amax=amax/(itime_scal+1)**2
1535 call predict_edrift(epdrift)
1536 ! write(iout,*) "amax=",amax,damax,epdrift,edriftmax,amax/(itime_scal+1)
1538 ! write (iout,*) "before enter if",scalel,icount_scale
1539 if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then
1540 ! write(iout,*) "I enter if"
1541 ! Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
1543 ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax)) &
1545 itime_scal=itime_scal+ifac_time
1546 ! fac_time=dmin1(damax/amax,0.5d0)
1547 fac_time=0.5d0**ifac_time
1548 d_time=d_time*fac_time
1549 if (lang.eq.2 .or. lang.eq.3) then
1551 ! write (iout,*) "Calling sd_verlet_setup: 1"
1552 ! Rescale the stochastic forces and recalculate or restore
1553 ! the matrices of tinker integrator
1554 if (itime_scal.gt.maxflag_stoch) then
1555 if (large) write (iout,'(a,i5,a)') &
1556 "Calculate matrices for stochastic step;",&
1557 " itime_scal ",itime_scal
1559 call sd_verlet_p_setup
1561 call sd_verlet_ciccotti_setup
1563 write (iout,'(2a,i3,a,i3,1h.)') &
1564 "Warning: cannot store matrices for stochastic",&
1565 " integration because the index",itime_scal,&
1566 " is greater than",maxflag_stoch
1567 write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",&
1568 " integration Langevin algorithm for better efficiency."
1569 else if (flag_stoch(itime_scal)) then
1570 if (large) write (iout,'(a,i5,a,l1)') &
1571 "Restore matrices for stochastic step; itime_scal ",&
1572 itime_scal," flag ",flag_stoch(itime_scal)
1575 pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
1576 afric_mat(i,j)=afric0_mat(i,j,itime_scal)
1577 vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
1578 prand_mat(i,j)=prand0_mat(i,j,itime_scal)
1579 vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
1580 vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
1584 if (large) write (iout,'(2a,i5,a,l1)') &
1585 "Calculate & store matrices for stochastic step;",&
1586 " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
1588 call sd_verlet_p_setup
1590 call sd_verlet_ciccotti_setup
1592 flag_stoch(ifac_time)=.true.
1595 pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
1596 afric0_mat(i,j,itime_scal)=afric_mat(i,j)
1597 vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
1598 prand0_mat(i,j,itime_scal)=prand_mat(i,j)
1599 vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
1600 vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
1604 fac_time=1.0d0/dsqrt(fac_time)
1606 stochforcvec(i)=fac_time*stochforcvec(i)
1609 else if (lang.eq.1) then
1610 ! Rescale the accelerations due to stochastic forces
1611 fac_time=1.0d0/dsqrt(fac_time)
1613 d_as_work(i)=d_as_work(i)*fac_time
1616 if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)') &
1617 "itime",itime," Timestep scaled down to ",&
1618 d_time," ifac_time",ifac_time," itime_scal",itime_scal
1620 ! Second step of the velocity Verlet algorithm
1625 else if (lang.eq.3) then
1627 call sd_verlet2_ciccotti
1629 else if (lang.eq.1) then
1634 if (rattle) call rattle2
1637 if (d_time.ne.d_time0) then
1640 if (lang.eq.2 .or. lang.eq.3) then
1641 if (large) write (iout,'(a)') &
1642 "Restore original matrices for stochastic step"
1643 ! write (iout,*) "Calling sd_verlet_setup: 2"
1644 ! Restore the matrices of tinker integrator if the time step has been restored
1647 pfric_mat(i,j)=pfric0_mat(i,j,0)
1648 afric_mat(i,j)=afric0_mat(i,j,0)
1649 vfric_mat(i,j)=vfric0_mat(i,j,0)
1650 prand_mat(i,j)=prand0_mat(i,j,0)
1651 vrand_mat1(i,j)=vrand0_mat1(i,j,0)
1652 vrand_mat2(i,j)=vrand0_mat2(i,j,0)
1660 ! Calculate the kinetic and the total energy and the kinetic temperature
1664 ! call kinetic1(EK1)
1665 ! write (iout,*) "step",itime," EK",EK," EK1",EK1
1667 ! Couple the system to Berendsen bath if needed
1668 if (tbf .and. lang.eq.0) then
1671 kinetic_T=2.0d0/(dimen3*Rb)*EK
1672 ! Backup the coordinates, velocities, and accelerations
1676 d_t_old(j,i)=d_t(j,i)
1677 d_a_old(j,i)=d_a(j,i)
1681 if (mod(itime,ntwe).eq.0 .and. large) then
1682 write (iout,*) "Velocities, step 2"
1684 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1685 (d_t(j,i+nres),j=1,3)
1690 end subroutine velverlet_step
1691 !-----------------------------------------------------------------------------
1692 subroutine RESPA_step(itime)
1693 !-------------------------------------------------------------------------------
1694 ! Perform a single RESPA step.
1695 !-------------------------------------------------------------------------------
1696 ! implicit real*8 (a-h,o-z)
1697 ! include 'DIMENSIONS'
1701 use control, only:tcpu
1703 ! use io_conf, only:cartprint
1706 integer :: IERROR,ERRCODE
1708 ! include 'COMMON.SETUP'
1709 ! include 'COMMON.CONTROL'
1710 ! include 'COMMON.VAR'
1711 ! include 'COMMON.MD'
1713 ! include 'COMMON.LANGEVIN'
1715 ! include 'COMMON.LANGEVIN.lang0'
1717 ! include 'COMMON.CHAIN'
1718 ! include 'COMMON.DERIV'
1719 ! include 'COMMON.GEO'
1720 ! include 'COMMON.LOCAL'
1721 ! include 'COMMON.INTERACT'
1722 ! include 'COMMON.IOUNITS'
1723 ! include 'COMMON.NAMES'
1724 ! include 'COMMON.TIME1'
1725 real(kind=8),dimension(0:n_ene) :: energia_short,energia_long
1726 real(kind=8),dimension(3) :: L,vcm,incr
1727 real(kind=8),dimension(3,0:2*nres) :: dc_old0,d_t_old0,d_a_old0 !(3,0:maxres2) maxres2=2*maxres
1728 logical :: PRINT_AMTS_MSG = .false.
1729 integer :: count,rstcount !ilen,
1731 character(len=50) :: tytul
1732 integer :: maxcount_scale = 10
1733 !el common /gucio/ cm
1734 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1735 !el common /stochcalc/ stochforcvec
1736 integer :: itt,i,j,itsplit,itime
1738 !el common /cipiszcze/ itt
1740 real(kind=8) :: epdrift,tt0,epdriftmax
1743 if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres)) !(MAXRES6) maxres6=6*maxres
1747 if (large.and. mod(itime,ntwe).eq.0) then
1748 write (iout,*) "***************** RESPA itime",itime
1749 write (iout,*) "Cartesian and internal coordinates: step 0"
1751 call pdbout(0.0d0,"cipiszcze",iout)
1753 write (iout,*) "Accelerations from long-range forces"
1755 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1756 (d_a(j,i+nres),j=1,3)
1758 write (iout,*) "Velocities, step 0"
1760 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1761 (d_t(j,i+nres),j=1,3)
1766 ! Perform the initial RESPA step (increment velocities)
1767 ! write (iout,*) "*********************** RESPA ini"
1770 if (mod(itime,ntwe).eq.0 .and. large) then
1771 write (iout,*) "Velocities, end"
1773 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1774 (d_t(j,i+nres),j=1,3)
1778 ! Compute the short-range forces
1784 ! 7/2/2009 commented out
1786 ! call etotal_short(energia_short)
1789 ! 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
1792 d_a(j,i)=d_a_short(j,i)
1796 if (large.and. mod(itime,ntwe).eq.0) then
1797 write (iout,*) "energia_short",energia_short(0)
1798 write (iout,*) "Accelerations from short-range forces"
1800 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1801 (d_a(j,i+nres),j=1,3)
1806 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1808 t_enegrad=t_enegrad+tcpu()-tt0
1813 d_t_old(j,i)=d_t(j,i)
1814 d_a_old(j,i)=d_a(j,i)
1817 ! 6/30/08 A-MTS: attempt at increasing the split number
1820 dc_old0(j,i)=dc_old(j,i)
1821 d_t_old0(j,i)=d_t_old(j,i)
1822 d_a_old0(j,i)=d_a_old(j,i)
1825 if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
1826 if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
1833 ! write (iout,*) "itime",itime," ntime_split",ntime_split
1834 ! Split the time step
1835 d_time=d_time0/ntime_split
1836 ! Perform the short-range RESPA steps (velocity Verlet increments of
1837 ! positions and velocities using short-range forces)
1838 ! write (iout,*) "*********************** RESPA split"
1839 do itsplit=1,ntime_split
1842 else if (lang.eq.2 .or. lang.eq.3) then
1844 call stochastic_force(stochforcvec)
1847 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
1849 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1854 ! First step of the velocity Verlet algorithm
1859 else if (lang.eq.3) then
1861 call sd_verlet1_ciccotti
1863 else if (lang.eq.1) then
1868 ! Build the chain from the newly calculated coordinates
1869 call chainbuild_cart
1870 if (rattle) call rattle1
1872 if (large.and. mod(itime,ntwe).eq.0) then
1873 write (iout,*) "***** ITSPLIT",itsplit
1874 write (iout,*) "Cartesian and internal coordinates: step 1"
1875 call pdbout(0.0d0,"cipiszcze",iout)
1878 write (iout,*) "Velocities, step 1"
1880 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1881 (d_t(j,i+nres),j=1,3)
1890 ! Calculate energy and forces
1892 call etotal_short(energia_short)
1893 ! AL 4/17/17: Exit itime_split loop when energy goes infinite
1894 if (energia_short(0).gt.0.99e20 .or. isnan(energia_short(0)) ) then
1895 if (PRINT_AMTS_MSG) &
1896 write (iout,*) "Infinities/NaNs in energia_short",energia_short(0),"; increasing ntime_split to",ntime_split
1897 ntime_split=ntime_split*2
1898 if (ntime_split.gt.maxtime_split) then
1901 "Cannot rescue the run; aborting job. Retry with a smaller time step"
1903 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1906 "Cannot rescue the run; terminating. Retry with a smaller time step"
1912 if (large.and. mod(itime,ntwe).eq.0) &
1913 call enerprint(energia_short)
1916 t_eshort=t_eshort+MPI_Wtime()-tt0
1918 t_eshort=t_eshort+tcpu()-tt0
1922 ! Get the new accelerations
1924 ! 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
1927 d_a_short(j,i)=d_a(j,i)
1931 if (large.and. mod(itime,ntwe).eq.0) then
1932 write (iout,*)"energia_short",energia_short(0)
1933 write (iout,*) "Accelerations from short-range forces"
1935 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1936 (d_a(j,i+nres),j=1,3)
1941 ! Determine maximum acceleration and scale down the timestep if needed
1943 amax=amax/ntime_split**2
1944 call predict_edrift(epdrift)
1945 if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0) &
1946 write (iout,*) "amax",amax," damax",damax,&
1947 " epdrift",epdrift," epdriftmax",epdriftmax
1948 ! Exit loop and try with increased split number if the change of
1949 ! acceleration is too big
1950 if (amax.gt.damax .or. epdrift.gt.edriftmax) then
1951 if (ntime_split.lt.maxtime_split) then
1953 ntime_split=ntime_split*2
1954 ! AL 4/17/17: We should exit the itime_split loop when acceleration change is too big
1958 dc_old(j,i)=dc_old0(j,i)
1959 d_t_old(j,i)=d_t_old0(j,i)
1960 d_a_old(j,i)=d_a_old0(j,i)
1963 if (PRINT_AMTS_MSG) then
1964 write (iout,*) "acceleration/energy drift too large",amax,&
1965 epdrift," split increased to ",ntime_split," itime",itime,&
1971 "Uh-hu. Bumpy landscape. Maximum splitting number",&
1973 " already reached!!! Trying to carry on!"
1977 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1979 t_enegrad=t_enegrad+tcpu()-tt0
1981 ! Second step of the velocity Verlet algorithm
1986 else if (lang.eq.3) then
1988 call sd_verlet2_ciccotti
1990 else if (lang.eq.1) then
1995 if (rattle) call rattle2
1996 ! Backup the coordinates, velocities, and accelerations
2000 d_t_old(j,i)=d_t(j,i)
2001 d_a_old(j,i)=d_a(j,i)
2008 ! Restore the time step
2010 ! Compute long-range forces
2017 call etotal_long(energia_long)
2018 if (energia_long(0).gt.0.99e20 .or. isnan(energia_long(0))) then
2021 "Infinitied/NaNs in energia_long, Aborting MPI job."
2023 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
2025 write (iout,*) "Infinitied/NaNs in energia_long, terminating."
2029 if (large.and. mod(itime,ntwe).eq.0) &
2030 call enerprint(energia_long)
2033 t_elong=t_elong+MPI_Wtime()-tt0
2035 t_elong=t_elong+tcpu()-tt0
2038 potE=potEcomp(0)-potEcomp(51)
2042 t_enegrad=t_enegrad+MPI_Wtime()-tt0
2044 t_enegrad=t_enegrad+tcpu()-tt0
2046 ! Compute accelerations from long-range forces
2048 if (large.and. mod(itime,ntwe).eq.0) then
2049 write (iout,*) "energia_long",energia_long(0)
2050 write (iout,*) "Cartesian and internal coordinates: step 2"
2052 call pdbout(0.0d0,"cipiszcze",iout)
2054 write (iout,*) "Accelerations from long-range forces"
2056 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
2057 (d_a(j,i+nres),j=1,3)
2059 write (iout,*) "Velocities, step 2"
2061 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2062 (d_t(j,i+nres),j=1,3)
2066 ! Compute the final RESPA step (increment velocities)
2067 ! write (iout,*) "*********************** RESPA fin"
2069 ! Compute the complete potential energy
2071 potEcomp(i)=energia_short(i)+energia_long(i)
2073 potE=potEcomp(0)-potEcomp(51)
2074 ! potE=energia_short(0)+energia_long(0)
2077 ! Calculate the kinetic and the total energy and the kinetic temperature
2080 ! Couple the system to Berendsen bath if needed
2081 if (tbf .and. lang.eq.0) then
2084 kinetic_T=2.0d0/(dimen3*Rb)*EK
2085 ! Backup the coordinates, velocities, and accelerations
2087 if (mod(itime,ntwe).eq.0 .and. large) then
2088 write (iout,*) "Velocities, end"
2090 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2091 (d_t(j,i+nres),j=1,3)
2096 end subroutine RESPA_step
2097 !-----------------------------------------------------------------------------
2098 subroutine RESPA_vel
2099 ! First and last RESPA step (incrementing velocities using long-range
2102 ! implicit real*8 (a-h,o-z)
2103 ! include 'DIMENSIONS'
2104 ! include 'COMMON.CONTROL'
2105 ! include 'COMMON.VAR'
2106 ! include 'COMMON.MD'
2107 ! include 'COMMON.CHAIN'
2108 ! include 'COMMON.DERIV'
2109 ! include 'COMMON.GEO'
2110 ! include 'COMMON.LOCAL'
2111 ! include 'COMMON.INTERACT'
2112 ! include 'COMMON.IOUNITS'
2113 ! include 'COMMON.NAMES'
2114 integer :: i,j,inres,mnum
2117 d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
2121 d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
2126 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2127 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2128 .and.(mnum.lt.4)) then
2131 d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
2136 end subroutine RESPA_vel
2137 !-----------------------------------------------------------------------------
2139 ! Applying velocity Verlet algorithm - step 1 to coordinates
2141 ! implicit real*8 (a-h,o-z)
2142 ! include 'DIMENSIONS'
2143 ! include 'COMMON.CONTROL'
2144 ! include 'COMMON.VAR'
2145 ! include 'COMMON.MD'
2146 ! include 'COMMON.CHAIN'
2147 ! include 'COMMON.DERIV'
2148 ! include 'COMMON.GEO'
2149 ! include 'COMMON.LOCAL'
2150 ! include 'COMMON.INTERACT'
2151 ! include 'COMMON.IOUNITS'
2152 ! include 'COMMON.NAMES'
2153 real(kind=8) :: adt,adt2
2154 integer :: i,j,inres,mnum
2157 write (iout,*) "VELVERLET1 START: DC"
2159 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),&
2160 (dc(j,i+nres),j=1,3)
2164 adt=d_a_old(j,0)*d_time
2166 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
2167 d_t_new(j,0)=d_t_old(j,0)+adt2
2168 d_t(j,0)=d_t_old(j,0)+adt
2172 adt=d_a_old(j,i)*d_time
2174 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
2175 d_t_new(j,i)=d_t_old(j,i)+adt2
2176 d_t(j,i)=d_t_old(j,i)+adt
2181 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2182 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2183 .and.(mnum.lt.4)) then
2186 adt=d_a_old(j,inres)*d_time
2188 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
2189 d_t_new(j,inres)=d_t_old(j,inres)+adt2
2190 d_t(j,inres)=d_t_old(j,inres)+adt
2195 write (iout,*) "VELVERLET1 END: DC"
2197 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),&
2198 (dc(j,i+nres),j=1,3)
2202 end subroutine verlet1
2203 !-----------------------------------------------------------------------------
2205 ! Step 2 of the velocity Verlet algorithm: update velocities
2207 ! implicit real*8 (a-h,o-z)
2208 ! include 'DIMENSIONS'
2209 ! include 'COMMON.CONTROL'
2210 ! include 'COMMON.VAR'
2211 ! include 'COMMON.MD'
2212 ! include 'COMMON.CHAIN'
2213 ! include 'COMMON.DERIV'
2214 ! include 'COMMON.GEO'
2215 ! include 'COMMON.LOCAL'
2216 ! include 'COMMON.INTERACT'
2217 ! include 'COMMON.IOUNITS'
2218 ! include 'COMMON.NAMES'
2219 integer :: i,j,inres,mnum
2222 d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
2226 d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
2231 ! iti=iabs(itype(i,mnum))
2232 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2233 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2234 .and.(mnum.lt.4)) then
2237 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
2242 end subroutine verlet2
2243 !-----------------------------------------------------------------------------
2244 subroutine sddir_precalc
2245 ! Applying velocity Verlet algorithm - step 1 to coordinates
2246 ! implicit real*8 (a-h,o-z)
2247 ! include 'DIMENSIONS'
2253 ! include 'COMMON.CONTROL'
2254 ! include 'COMMON.VAR'
2255 ! include 'COMMON.MD'
2257 ! include 'COMMON.LANGEVIN'
2259 ! include 'COMMON.LANGEVIN.lang0'
2261 ! include 'COMMON.CHAIN'
2262 ! include 'COMMON.DERIV'
2263 ! include 'COMMON.GEO'
2264 ! include 'COMMON.LOCAL'
2265 ! include 'COMMON.INTERACT'
2266 ! include 'COMMON.IOUNITS'
2267 ! include 'COMMON.NAMES'
2268 ! include 'COMMON.TIME1'
2269 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
2270 !el common /stochcalc/ stochforcvec
2271 real(kind=8) :: time00
2274 ! Compute friction and stochastic forces
2278 if (large) print *,"before friction_force"
2280 if (large) print *,"after friction_force"
2281 time_fric=time_fric+MPI_Wtime()-time00
2283 call stochastic_force(stochforcvec)
2284 time_stoch=time_stoch+MPI_Wtime()-time00
2287 ! Compute the acceleration due to friction forces (d_af_work) and stochastic
2288 ! forces (d_as_work)
2290 ! call ginv_mult(fric_work, d_af_work)
2291 ! call ginv_mult(stochforcvec, d_as_work)
2293 write(iout,*) "forces before fivediaginv"
2295 write(iout,*) "fricwork",i,fric_work(i)
2297 call fivediaginv_mult(dimen,fric_work, d_af_work)
2298 call fivediaginv_mult(dimen,stochforcvec, d_as_work)
2300 write(iout,*),"dimen",dimen
2302 write(iout,*) "fricwork",fric_work(i), d_af_work(i)
2303 write(iout,*) "stochforcevec", stochforcvec(i), d_as_work(i)
2307 call ginv_mult(fric_work, d_af_work)
2308 call ginv_mult(stochforcvec, d_as_work)
2312 end subroutine sddir_precalc
2313 !-----------------------------------------------------------------------------
2314 subroutine sddir_verlet1
2315 ! Applying velocity Verlet algorithm - step 1 to velocities
2318 ! implicit real*8 (a-h,o-z)
2319 ! include 'DIMENSIONS'
2320 ! include 'COMMON.CONTROL'
2321 ! include 'COMMON.VAR'
2322 ! include 'COMMON.MD'
2324 ! include 'COMMON.LANGEVIN'
2326 ! include 'COMMON.LANGEVIN.lang0'
2328 ! include 'COMMON.CHAIN'
2329 ! include 'COMMON.DERIV'
2330 ! include 'COMMON.GEO'
2331 ! include 'COMMON.LOCAL'
2332 ! include 'COMMON.INTERACT'
2333 ! include 'COMMON.IOUNITS'
2334 ! include 'COMMON.NAMES'
2335 ! Revised 3/31/05 AL: correlation between random contributions to
2336 ! position and velocity increments included.
2337 real(kind=8) :: sqrt13 = 0.57735026918962576451d0 ! 1/sqrt(3)
2338 real(kind=8) :: adt,adt2
2339 integer :: i,j,ind,inres,mnum
2341 ! Add the contribution from BOTH friction and stochastic force to the
2342 ! coordinates, but ONLY the contribution from the friction forces to velocities
2345 adt=(d_a_old(j,0)+d_af_work(j))*d_time
2346 adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
2347 ! write(iout,*) i,"adt",adt,"ads",adt2,d_a_old(j,0),d_af_work(j),d_time
2348 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
2349 d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
2350 d_t(j,0)=d_t_old(j,0)+adt
2355 adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
2356 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
2357 ! write(iout,*) i,"adt",adt,"ads",adt2,d_a_old(j,i),d_af_work(ind+j)
2358 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
2359 d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
2360 d_t(j,i)=d_t_old(j,i)+adt
2366 ! iti=iabs(itype(i,mnum))
2367 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2368 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2369 .and.(mnum.lt.4)) then
2372 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
2373 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
2374 ! write(iout,*) i,"adt",adt,"ads",adt2,d_a_old(j,inres),d_af_work(ind+j)
2375 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
2376 d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
2377 d_t(j,inres)=d_t_old(j,inres)+adt
2384 end subroutine sddir_verlet1
2385 !-----------------------------------------------------------------------------
2386 subroutine sddir_verlet2
2387 ! Calculating the adjusted velocities for accelerations
2390 ! implicit real*8 (a-h,o-z)
2391 ! include 'DIMENSIONS'
2392 ! include 'COMMON.CONTROL'
2393 ! include 'COMMON.VAR'
2394 ! include 'COMMON.MD'
2396 ! include 'COMMON.LANGEVIN'
2398 ! include 'COMMON.LANGEVIN.lang0'
2400 ! include 'COMMON.CHAIN'
2401 ! include 'COMMON.DERIV'
2402 ! include 'COMMON.GEO'
2403 ! include 'COMMON.LOCAL'
2404 ! include 'COMMON.INTERACT'
2405 ! include 'COMMON.IOUNITS'
2406 ! include 'COMMON.NAMES'
2407 real(kind=8),dimension(:),allocatable :: stochforcvec,d_as_work1 !(MAXRES6) maxres6=6*maxres
2408 real(kind=8) :: cos60 = 0.5d0, sin60 = 0.86602540378443864676d0
2409 integer :: i,j,ind,inres,mnum
2410 if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres))
2411 if (.not.allocated(d_as_work1)) allocate(d_as_work1(6*nres))
2412 ! Revised 3/31/05 AL: correlation between random contributions to
2413 ! position and velocity increments included.
2414 ! The correlation coefficients are calculated at low-friction limit.
2415 ! Also, friction forces are now not calculated with new velocities.
2417 ! call friction_force
2418 call stochastic_force(stochforcvec)
2420 ! Compute the acceleration due to friction forces (d_af_work) and stochastic
2421 ! forces (d_as_work)
2424 call fivediaginv_mult(6*nres,stochforcvec, d_as_work1)
2426 call ginv_mult(stochforcvec, d_as_work1)
2433 d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j)) &
2434 +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
2439 d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j)) &
2440 +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
2446 ! iti=iabs(itype(i,mnum))
2447 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2448 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2449 .and.(mnum.lt.4)) then
2452 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres) &
2453 +d_af_work(ind+j))+sin60*d_as_work(ind+j) &
2454 +cos60*d_as_work1(ind+j))*d_time
2460 end subroutine sddir_verlet2
2461 !-----------------------------------------------------------------------------
2462 subroutine max_accel
2464 ! Find the maximum difference in the accelerations of the the sites
2465 ! at the beginning and the end of the time step.
2468 ! implicit real*8 (a-h,o-z)
2469 ! include 'DIMENSIONS'
2470 ! include 'COMMON.CONTROL'
2471 ! include 'COMMON.VAR'
2472 ! include 'COMMON.MD'
2473 ! include 'COMMON.CHAIN'
2474 ! include 'COMMON.DERIV'
2475 ! include 'COMMON.GEO'
2476 ! include 'COMMON.LOCAL'
2477 ! include 'COMMON.INTERACT'
2478 ! include 'COMMON.IOUNITS'
2479 real(kind=8),dimension(3) :: aux,accel,accel_old
2480 real(kind=8) :: dacc
2484 ! aux(j)=d_a(j,0)-d_a_old(j,0)
2485 accel_old(j)=d_a_old(j,0)
2492 ! 7/3/08 changed to asymmetric difference
2494 ! accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
2495 accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
2496 accel(j)=accel(j)+0.5d0*d_a(j,i)
2497 ! if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
2498 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
2499 dacc=dabs(accel(j)-accel_old(j))
2500 ! write (iout,*) i,dacc
2501 if (dacc.gt.amax) amax=dacc
2509 accel_old(j)=d_a_old(j,0)
2514 accel_old(j)=accel_old(j)+d_a_old(j,1)
2515 accel(j)=accel(j)+d_a(j,1)
2520 ! iti=iabs(itype(i,mnum))
2521 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2522 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2523 .and.(mnum.lt.4)) then
2525 ! accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
2526 accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
2527 accel(j)=accel(j)+d_a(j,i+nres)
2531 ! if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
2532 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
2533 dacc=dabs(accel(j)-accel_old(j))
2534 ! write (iout,*) "side-chain",i,dacc
2535 if (dacc.gt.amax) amax=dacc
2539 accel_old(j)=accel_old(j)+d_a_old(j,i)
2540 accel(j)=accel(j)+d_a(j,i)
2541 ! aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
2545 end subroutine max_accel
2546 !-----------------------------------------------------------------------------
2547 subroutine predict_edrift(epdrift)
2549 ! Predict the drift of the potential energy
2552 use control_data, only: lmuca
2553 ! implicit real*8 (a-h,o-z)
2554 ! include 'DIMENSIONS'
2555 ! include 'COMMON.CONTROL'
2556 ! include 'COMMON.VAR'
2557 ! include 'COMMON.MD'
2558 ! include 'COMMON.CHAIN'
2559 ! include 'COMMON.DERIV'
2560 ! include 'COMMON.GEO'
2561 ! include 'COMMON.LOCAL'
2562 ! include 'COMMON.INTERACT'
2563 ! include 'COMMON.IOUNITS'
2564 ! include 'COMMON.MUCA'
2565 real(kind=8) :: epdrift,epdriftij
2567 ! Drift of the potential energy
2573 epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
2574 if (lmuca) epdriftij=epdriftij*factor
2575 ! write (iout,*) "back",i,j,epdriftij
2576 if (epdriftij.gt.epdrift) epdrift=epdriftij
2580 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1.and.&
2581 molnum(i).lt.4) then
2584 dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
2585 if (lmuca) epdriftij=epdriftij*factor
2586 ! write (iout,*) "side",i,j,epdriftij
2587 if (epdriftij.gt.epdrift) epdrift=epdriftij
2591 epdrift=0.5d0*epdrift*d_time*d_time
2592 ! write (iout,*) "epdrift",epdrift
2594 end subroutine predict_edrift
2595 !-----------------------------------------------------------------------------
2596 subroutine verlet_bath
2598 ! Coupling to the thermostat by using the Berendsen algorithm
2601 ! implicit real*8 (a-h,o-z)
2602 ! include 'DIMENSIONS'
2603 ! include 'COMMON.CONTROL'
2604 ! include 'COMMON.VAR'
2605 ! include 'COMMON.MD'
2606 ! include 'COMMON.CHAIN'
2607 ! include 'COMMON.DERIV'
2608 ! include 'COMMON.GEO'
2609 ! include 'COMMON.LOCAL'
2610 ! include 'COMMON.INTERACT'
2611 ! include 'COMMON.IOUNITS'
2612 ! include 'COMMON.NAMES'
2613 real(kind=8) :: T_half,fact
2614 integer :: i,j,inres,mnum
2616 T_half=2.0d0/(dimen3*Rb)*EK
2617 fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
2618 ! write(iout,*) "T_half", T_half
2619 ! write(iout,*) "EK", EK
2620 ! write(iout,*) "fact", fact
2622 d_t(j,0)=fact*d_t(j,0)
2626 d_t(j,i)=fact*d_t(j,i)
2631 ! iti=iabs(itype(i,mnum))
2632 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2633 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2634 .and.(mnum.lt.4)) then
2637 d_t(j,inres)=fact*d_t(j,inres)
2642 end subroutine verlet_bath
2643 !-----------------------------------------------------------------------------
2645 ! Set up the initial conditions of a MD simulation
2648 use control, only:tcpu
2649 !el use io_basic, only:ilen
2652 use minimm, only:minim_dc,minimize,sc_move
2653 use io_config, only:readrst
2654 use io, only:statout
2655 use random, only: iran_num
2656 ! implicit real*8 (a-h,o-z)
2657 ! include 'DIMENSIONS'
2660 character(len=16) :: form
2661 integer :: IERROR,ERRCODE
2663 integer :: iranmin,itrial,itmp,n_model_try,k, &
2665 integer, dimension(:),allocatable :: list_model_try
2666 integer, dimension(0:nodes-1) :: i_start_models
2667 ! include 'COMMON.SETUP'
2668 ! include 'COMMON.CONTROL'
2669 ! include 'COMMON.VAR'
2670 ! include 'COMMON.MD'
2672 ! include 'COMMON.LANGEVIN'
2674 ! include 'COMMON.LANGEVIN.lang0'
2676 ! include 'COMMON.CHAIN'
2677 ! include 'COMMON.DERIV'
2678 ! include 'COMMON.GEO'
2679 ! include 'COMMON.LOCAL'
2680 ! include 'COMMON.INTERACT'
2681 ! include 'COMMON.IOUNITS'
2682 ! include 'COMMON.NAMES'
2683 ! include 'COMMON.REMD'
2684 real(kind=8),dimension(0:n_ene) :: energia_long,energia_short,energia
2685 real(kind=8),dimension(3) :: vcm,incr,L
2686 real(kind=8) :: xv,sigv,lowb,highb
2687 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
2688 character(len=256) :: qstr
2691 character(len=50) :: tytul
2692 logical :: file_exist
2693 !el common /gucio/ cm
2694 integer :: i,j,ipos,iq,iw,nft_sc,iretcode,ierr,mnum,itime
2698 real(kind=8) :: etot,tt0
2702 ! write(iout,*) "d_time", d_time
2703 ! Compute the standard deviations of stochastic forces for Langevin dynamics
2704 ! if the friction coefficients do not depend on surface area
2705 if (lang.gt.0 .and. .not.surfarea) then
2708 stdforcp(i)=stdfp(mnum)*dsqrt(gamp(mnum))
2712 stdforcsc(i)=stdfsc(iabs(itype(i,mnum)),mnum) &
2713 *dsqrt(gamsc(iabs(itype(i,mnum)),mnum))
2717 ! Open the pdb file for snapshotshots
2720 if (ilen(tmpdir).gt.0) &
2721 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2722 liczba(:ilen(liczba))//".pdb")
2724 file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2728 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) &
2729 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2730 liczba(:ilen(liczba))//".x")
2731 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2734 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) &
2735 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2736 liczba(:ilen(liczba))//".cx")
2737 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2743 if (ilen(tmpdir).gt.0) &
2744 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
2745 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
2747 if (ilen(tmpdir).gt.0) &
2748 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
2749 cartname=prefix(:ilen(prefix))//"_MD.cx"
2753 write (qstr,'(256(1h ))')
2756 iq = qinfrag(i,iset)*10
2757 iw = wfrag(i,iset)/100
2759 if(me.eq.king.or..not.out1file) &
2760 write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
2761 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
2766 iq = qinpair(i,iset)*10
2767 iw = wpair(i,iset)/100
2769 if(me.eq.king.or..not.out1file) &
2770 write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
2771 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
2775 ! pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
2777 ! cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
2779 ! cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
2781 ! statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
2785 if (restart1file) then
2787 inquire(file=mremd_rst_name,exist=file_exist)
2788 write (*,*) me," Before broadcast: file_exist",file_exist
2790 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,&
2793 write (*,*) me," After broadcast: file_exist",file_exist
2794 ! inquire(file=mremd_rst_name,exist=file_exist)
2795 if(me.eq.king.or..not.out1file) &
2796 write(iout,*) "Initial state read by master and distributed"
2798 if (ilen(tmpdir).gt.0) &
2799 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_' &
2800 //liczba(:ilen(liczba))//'.rst')
2801 inquire(file=rest2name,exist=file_exist)
2804 if(.not.restart1file) then
2805 if(me.eq.king.or..not.out1file) &
2806 write(iout,*) "Initial state will be read from file ",&
2807 rest2name(:ilen(rest2name))
2810 call rescale_weights(t_bath)
2812 if(me.eq.king.or..not.out1file)then
2813 if (restart1file) then
2814 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),&
2817 write(iout,*) "File ",rest2name(:ilen(rest2name)),&
2820 write(iout,*) "Initial velocities randomly generated"
2827 ! Generate initial velocities
2828 if(me.eq.king.or..not.out1file) &
2829 write(iout,*) "Initial velocities randomly generated"
2834 ! rest2name = prefix(:ilen(prefix))//'.rst'
2835 if(me.eq.king.or..not.out1file)then
2836 write (iout,*) "Initial velocities"
2838 write (iout,'(i6,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2839 (d_t(j,i+nres),j=1,3)
2841 ! Zeroing the total angular momentum of the system
2842 write(iout,*) "Calling the zero-angular momentum subroutine"
2845 ! Getting the potential energy and forces and velocities and accelerations
2847 ! write (iout,*) "velocity of the center of the mass:"
2848 ! write (iout,*) (vcm(j),j=1,3)
2850 d_t(j,0)=d_t(j,0)-vcm(j)
2852 ! Removing the velocity of the center of mass
2854 if(me.eq.king.or..not.out1file)then
2855 write (iout,*) "vcm right after adjustment:"
2856 write (iout,*) (vcm(j),j=1,3)
2863 if ((.not.rest).or.(forceminim)) then
2864 if (forceminim) call chainbuild_cart
2866 if(iranconf.ne.0 .or.indpdb.gt.0.and..not.unres_pdb .or.preminim) then
2868 print *, 'Calling OVERLAP_SC'
2869 call overlap_sc(fail)
2870 print *,'after OVERLAP'
2873 print *,'call SC_MOVE'
2874 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
2875 print *,'SC_move',nft_sc,etot
2876 if(me.eq.king.or..not.out1file) &
2877 write(iout,*) 'SC_move',nft_sc,etot
2881 print *, 'Calling MINIM_DC'
2882 call minim_dc(etot,iretcode,nfun)
2884 call geom_to_var(nvar,varia)
2885 print *,'Calling MINIMIZE.'
2886 call minimize(etot,varia,iretcode,nfun)
2887 call var_to_geom(nvar,varia)
2889 write(iout,*) "just before minimin"
2891 if(me.eq.king.or..not.out1file) &
2892 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2894 write(iout,*) "just after minimin"
2896 if(iranconf.ne.0) then
2897 !c 8/22/17 AL Loop to produce a low-energy random conformation
2900 if(me.eq.king.or..not.out1file) &
2901 write (iout,*) 'Calling OVERLAP_SC'
2902 call overlap_sc(fail)
2903 endif !endif overlap
2906 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
2907 print *,'SC_move',nft_sc,etot
2908 if(me.eq.king.or..not.out1file) &
2909 write(iout,*) 'SC_move',nft_sc,etot
2913 print *, 'Calling MINIM_DC'
2914 call minim_dc(etot,iretcode,nfun)
2915 call int_from_cart1(.false.)
2917 call geom_to_var(nvar,varia)
2918 print *,'Calling MINIMIZE.'
2919 call minimize(etot,varia,iretcode,nfun)
2920 call var_to_geom(nvar,varia)
2922 if(me.eq.king.or..not.out1file) &
2923 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2924 write(iout,*) "just after minimin"
2926 if (isnan(etot) .or. etot.gt.4.0d6) then
2927 write (iout,*) "Energy too large",etot, &
2928 " trying another random conformation"
2931 call gen_rand_conf(itmp,*30)
2933 30 write (iout,*) 'Failed to generate random conformation', &
2935 write (*,*) 'Processor:',me, &
2936 ' Failed to generate random conformation',&
2945 write (iout,'(a,i3,a)') 'Processor:',me, &
2946 ' error in generating random conformation.'
2947 write (*,'(a,i3,a)') 'Processor:',me, &
2948 ' error in generating random conformation.'
2951 ! call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
2952 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
2962 write (iout,'(a,i3,a)') 'Processor:',me, &
2963 ' failed to generate a low-energy random conformation.'
2964 write (*,'(a,i3,a,f10.3)') 'Processor:',me, &
2965 ' failed to generate a low-energy random conformation.',etot
2969 ! call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
2970 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
2975 else if (preminim) then
2976 if (start_from_model) then
2980 do while (fail .and. n_model_try.lt.nmodel_start)
2981 write (iout,*) "n_model_try",n_model_try
2983 i_model=iran_num(1,nmodel_start)
2985 if (i_model.eq.list_model_try(k)) exit
2987 if (k.gt.n_model_try) exit
2989 n_model_try=n_model_try+1
2990 list_model_try(n_model_try)=i_model
2991 if (me.eq.king .or. .not. out1file) &
2992 write (iout,*) 'Trying to start from model ',&
2993 pdbfiles_chomo(i_model)(:ilen(pdbfiles_chomo(i_model)))
2996 c(j,i)=chomo(j,i,i_model)
2999 call int_from_cart(.true.,.false.)
3000 call sc_loc_geom(.false.)
3004 dc(j,i)=c(j,i+1)-c(j,i)
3005 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
3010 dc(j,i+nres)=c(j,i+nres)-c(j,i)
3011 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
3014 if (me.eq.king.or..not.out1file) then
3015 write (iout,*) "Energies before removing overlaps"
3016 call etotal(energia(0))
3017 call enerprint(energia(0))
3019 ! Remove SC overlaps if requested
3021 write (iout,*) 'Calling OVERLAP_SC'
3022 call overlap_sc(fail)
3025 "Failed to remove overlap from model",i_model
3029 if (me.eq.king.or..not.out1file) then
3030 write (iout,*) "Energies after removing overlaps"
3031 call etotal(energia(0))
3032 call enerprint(energia(0))
3035 ! Search for better SC rotamers if requested
3037 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
3038 print *,'SC_move',nft_sc,etot
3039 if (me.eq.king.or..not.out1file)&
3040 write(iout,*) 'SC_move',nft_sc,etot
3042 call etotal(energia(0))
3045 call MPI_Gather(i_model,1,MPI_INTEGER,i_start_models(0),&
3046 1,MPI_INTEGER,king,CG_COMM,IERROR)
3047 if (n_model_try.gt.nmodel_start .and.&
3048 (me.eq.king .or. out1file)) then
3050 "All models have irreparable overlaps. Trying randoms starts."
3052 i_model=nmodel_start+1
3056 ! Remove SC overlaps if requested
3058 write (iout,*) 'Calling OVERLAP_SC'
3059 call overlap_sc(fail)
3062 "Failed to remove overlap"
3065 if (me.eq.king.or..not.out1file) then
3066 write (iout,*) "Energies after removing overlaps"
3067 call etotal(energia(0))
3068 call enerprint(energia(0))
3071 ! 8/22/17 AL Minimize initial structure
3073 if (me.eq.king.or..not.out1file) write(iout,*)&
3074 'Minimizing initial PDB structure: Calling MINIM_DC'
3075 call minim_dc(etot,iretcode,nfun)
3077 if (me.eq.king.or..not.out1file)&
3078 write(iout,*) 'LBFGS return code is ',statusbf,' eval ',nfun
3081 call geom_to_var(nvar,varia)
3082 if(me.eq.king.or..not.out1file) write (iout,*)&
3083 'Minimizing initial PDB structure: Calling MINIMIZE.'
3084 call minimize(etot,varia,iretcode,nfun)
3085 call var_to_geom(nvar,varia)
3087 if (me.eq.king.or..not.out1file)&
3088 write(iout,*) 'LBFGS return code is ',statusbf,' eval ',nfun
3089 if(me.eq.king.or..not.out1file)&
3090 write(iout,*) 'LBFGS return code is ',statusbf,' eval ',nfun
3092 if (me.eq.king.or..not.out1file)&
3093 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
3094 if(me.eq.king.or..not.out1file)&
3095 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
3099 if (nmodel_start.gt.0 .and. me.eq.king) then
3100 write (iout,'(a)') "Task Starting model"
3102 if (i_start_models(i).gt.nmodel_start) then
3103 write (iout,'(i4,2x,a)') i,"RANDOM STRUCTURE"
3105 write(iout,'(i4,2x,a)')i,pdbfiles_chomo(i_start_models(i)) &
3106 (:ilen(pdbfiles_chomo(i_start_models(i))))
3111 call chainbuild_cart
3113 write(iout,*) "just after kinetic"
3118 kinetic_T=2.0d0/(dimen3*Rb)*EK
3119 if(me.eq.king.or..not.out1file)then
3120 write(iout,*) "just after verlet_bath"
3130 write(iout,*) "before ETOTAL"
3131 call etotal(potEcomp)
3132 if (large) call enerprint(potEcomp)
3135 t_etotal=t_etotal+MPI_Wtime()-tt0
3137 t_etotal=t_etotal+tcpu()-tt0
3142 write(iout,*) "before lagrangian"
3144 write(iout,*) "before max_accel"
3146 if (amax*d_time .gt. dvmax) then
3147 d_time=d_time*dvmax/amax
3148 if(me.eq.king.or..not.out1file) write (iout,*) &
3149 "Time step reduced to",d_time,&
3150 " because of too large initial acceleration."
3152 if(me.eq.king.or..not.out1file)then
3153 write(iout,*) "Potential energy and its components"
3154 call enerprint(potEcomp)
3155 ! write(iout,*) (potEcomp(i),i=0,n_ene)
3157 potE=potEcomp(0)-potEcomp(51)
3161 if (ntwe.ne.0) call statout(itime)
3162 if(me.eq.king.or..not.out1file) &
3163 write (iout,'(/a/3(a25,1pe14.5/))') "Initial:", &
3164 " Kinetic energy",EK," Potential energy",potE, &
3165 " Total energy",totE," Maximum acceleration ", &
3168 write (iout,*) "Initial coordinates"
3170 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),&
3173 write (iout,*) "Initial dC"
3175 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),&
3176 (dc(j,i+nres),j=1,3)
3178 write (iout,*) "Initial velocities"
3179 write (iout,"(13x,' backbone ',23x,' side chain')")
3181 write (iout,'(i6,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
3182 (d_t(j,i+nres),j=1,3)
3184 write (iout,*) "Initial accelerations"
3186 ! write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
3187 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),&
3188 (d_a(j,i+nres),j=1,3)
3194 d_t_old(j,i)=d_t(j,i)
3195 d_a_old(j,i)=d_a(j,i)
3197 ! write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
3206 call etotal_short(energia_short)
3207 if (large) call enerprint(potEcomp)
3210 t_eshort=t_eshort+MPI_Wtime()-tt0
3212 t_eshort=t_eshort+tcpu()-tt0
3217 if(.not.out1file .and. large) then
3218 write (iout,*) "energia_long",energia_long(0),&
3219 " energia_short",energia_short(0),&
3220 " total",energia_long(0)+energia_short(0)
3221 write (iout,*) "Initial fast-force accelerations"
3223 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
3224 (d_a(j,i+nres),j=1,3)
3227 ! 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
3230 d_a_short(j,i)=d_a(j,i)
3239 call etotal_long(energia_long)
3240 if (large) call enerprint(potEcomp)
3243 t_elong=t_elong+MPI_Wtime()-tt0
3245 t_elong=t_elong+tcpu()-tt0
3250 if(.not.out1file .and. large) then
3251 write (iout,*) "energia_long",energia_long(0)
3252 write (iout,*) "Initial slow-force accelerations"
3254 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
3255 (d_a(j,i+nres),j=1,3)
3259 t_enegrad=t_enegrad+MPI_Wtime()-tt0
3261 t_enegrad=t_enegrad+tcpu()-tt0
3265 end subroutine init_MD
3266 !-----------------------------------------------------------------------------
3267 subroutine random_vel
3269 ! implicit real*8 (a-h,o-z)
3271 use random, only:anorm_distr
3273 ! include 'DIMENSIONS'
3274 ! include 'COMMON.CONTROL'
3275 ! include 'COMMON.VAR'
3276 ! include 'COMMON.MD'
3278 ! include 'COMMON.LANGEVIN'
3280 ! include 'COMMON.LANGEVIN.lang0'
3282 ! include 'COMMON.CHAIN'
3283 ! include 'COMMON.DERIV'
3284 ! include 'COMMON.GEO'
3285 ! include 'COMMON.LOCAL'
3286 ! include 'COMMON.INTERACT'
3287 ! include 'COMMON.IOUNITS'
3288 ! include 'COMMON.NAMES'
3289 ! include 'COMMON.TIME1'
3290 real(kind=8) :: xv,sigv,lowb,highb ,Ek1
3292 integer ichain,n,innt,inct,ibeg,ierr,innt_org
3293 real(kind=8) ,allocatable, dimension(:):: work
3294 integer,allocatable,dimension(:) :: iwork
3295 ! double precision Ghalf(mmaxres2_chain),Geigen(maxres2_chain),&
3296 ! Gvec(maxres2_chain,maxres2_chain)
3297 ! common /przechowalnia/Ghalf,Geigen,Gvec
3299 ! double precision inertia(maxres2_chain,maxres2_chain)
3304 real(kind=8) ,allocatable, dimension(:) :: xsolv,DML,rs
3305 real(kind=8) :: sumx,Ek2,Ek3,aux,masinv
3307 real(kind=8) ,allocatable, dimension(:) :: rsold
3308 real (kind=8),allocatable,dimension(:,:) :: matold,inertia
3312 integer :: i,j,ii,k,mark,imark,mnum,nres2
3313 integer(kind=8) :: ind
3314 ! Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
3316 ! First generate velocities in the eigenspace of the G matrix
3317 ! write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
3320 if(.not.allocated(work)) then
3321 allocate(work(48*nres))
3322 allocate(iwork(6*nres))
3324 print *,"IN RANDOM VEL"
3326 ! print *,size(ghalf)
3329 write (iout,*) "Random_vel, fivediag"
3331 allocate(inertia(2*nres,2*nres))
3338 write(iout,*), "nchain",nchain
3342 ! if(.not.allocated(ghalf)) print *,"COCO"
3343 ! if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2))
3345 n=dimen_chain(ichain)
3346 innt=iposd_chain(ichain)
3348 innt_org=chain_border(1,ichain)
3349 if ((molnum(innt_org).eq.5).or.(molnum(innt_org).eq.4)) go to 137
3350 if(.not.allocated(ghalf)) print *,"COCO"
3351 if(.not.allocated(Ghalf)) allocate(Ghalf(1300*(1300+1)/2))
3355 write (iout,*) "Chain",ichain," n",n," start",innt
3357 if (i.lt.inct-1) then
3358 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i),DU1orig(i),&
3360 else if (i.eq.inct-1) then
3361 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i),DU1orig(i)
3363 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i)
3368 ghalf(ind+1)=dmorig(innt)
3369 ghalf(ind+2)=du1orig(innt)
3370 ghalf(ind+3)=dmorig(innt+1)
3374 write (iout,*) "i",i," ind",ind," indu2",innt+i-2,&
3375 " indu1",innt+i-1," indm",innt+i
3376 ghalf(ind+1)=du2orig(innt-1+i-2)
3377 ghalf(ind+2)=du1orig(innt-1+i-1)
3378 ghalf(ind+3)=dmorig(innt-1+i)
3379 !c write (iout,'(3(a,i2,1x))') "DU2",innt-1+i-2,
3380 !c & "DU1",innt-1+i-1,"DM ",innt-1+i
3388 inertia(i,j)=ghalf(ind)
3389 inertia(j,i)=ghalf(ind)
3394 write (iout,*) "Chain ",ichain," ind",ind," dim",n*(n+1)/2
3395 write (iout,*) "Five-diagonal inertia matrix, lower triangle"
3396 ! call matoutr(n,ghalf)
3398 call gldiag(nres*2,n,n,Ghalf,work,Geigen,Gvec,ierr,iwork)
3400 write (iout,'(//a,i3)')&
3401 "Eigenvectors and eigenvalues of the G matrix chain",ichain
3402 call eigout(n,n,nres*2,nres*2,Gvec,Geigen)
3405 !c check diagonalization
3411 aux=aux+gvec(k,i)*gvec(l,j)*inertia(k,l)
3415 write (iout,*) i,j,aux,geigen(i)
3417 write (iout,*) i,j,aux
3423 write(iout,*) "HERE,",n,innt
3424 innt_org=chain_border(1,ichain)
3430 mnum=molnum(innt_org)
3431 if (molnum(innt_org).ge.4) geigen(i)=3.0/msc(itype(innt_org+i-1,mnum),mnum)
3432 ! if (molnum(innt).eq.5) write(iout,*) "typ",i,innt-1+i,itype(innt+i-1,5)
3433 sigv=dsqrt((Rb*t_bath)/geigen(i))
3436 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
3437 EK=EK+0.5d0*geigen(i)*d_t_work_new(ii)**2
3438 write (iout,*) "i",i," ii",ii," geigen",geigen(i), &
3439 " d_t_work_new",d_t_work_new(ii),innt_org+i-1
3442 if (molnum(innt_org).ge.4) then
3443 mnum=molnum(innt_org)
3448 masinv=1.0d0/msc(itype(innt_org+i-1,mnum),mnum)
3449 d_t_work(ind)=d_t_work(ind)&
3450 +masinv*d_t_work_new((i-1)*3+k)
3460 d_t_work(ind)=d_t_work(ind)&
3461 +Gvec(i,j)*d_t_work_new((j-1)*3+k)
3471 aux=aux+inertia(i,j)*d_t_work(3*(i-1)+k)*d_t_work(3*(j-1)+k)
3477 !c Transfer to the d_t vector
3478 innt=chain_border(1,ichain)
3479 inct=chain_border(2,ichain)
3481 !c write (iout,*) "ichain",ichain," innt",innt," inct",inct
3485 d_t(j,i)=d_t_work(ind)
3488 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum).and.mnum.le.2) then
3491 d_t(j,i+nres)=d_t_work(ind)
3498 write (iout,*) "Random velocities in the Calpha,SC space"
3501 write (iout,'(a3,1h(,i5,1h),3f10.5,3x,3f10.5)')&
3502 restyp(itype(i,mnum),mnum),i,(d_t(j,i),j=1,3),(d_t(j,i+nres),j=1,3)
3505 call kinetic_CASC(Ek1)
3507 ! Transform the velocities to virtual-bond space
3516 if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum).or.mnum.ge.3) then
3518 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
3522 d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
3523 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
3536 write (iout,*) "Random vel after 1st transf the Calpha,SC space"
3537 write (iout,'(3hORG,1h(,i5,1h),3f10.5)') 0,(d_t(j,0),j=1,3)
3540 write (iout,'(a3,1h(,i5,1h),3f10.5,3x,3f10.5)')&
3541 restyp(itype(i,mnum),mnum),i,(d_t(j,i),j=1,3),(d_t(j,i+nres),j=1,3)
3545 !c d_a(:,0)=d_a(:,1)
3547 !c write (iout,*) "Shifting accelerations"
3549 write(iout,*) "nchain",ichain,chain_border1(1,ichain),molnum(chain_border1(1,ichain))
3550 if (molnum(chain_border1(1,ichain)+1).eq.5) cycle
3551 !c write (iout,*) "ichain",chain_border1(1,ichain)-1,
3552 !c & chain_border1(1,ichain)
3553 d_t(:,chain_border1(1,ichain)-1)=d_t(:,chain_border1(1,ichain))
3554 d_t(:,chain_border1(1,ichain))=0.0d0
3556 !c write (iout,*) "Adding accelerations"
3558 if (molnum(chain_border1(1,ichain)+1).eq.5) cycle
3559 !c write (iout,*) "chain",ichain,chain_border1(1,ichain)-1,
3560 !c & chain_border(2,ichain-1)
3561 d_t(:,chain_border1(1,ichain)-1)=&
3562 d_t(:,chain_border1(1,ichain)-1)+d_t(:,chain_border(2,ichain-1))
3563 d_t(:,chain_border(2,ichain-1))=0.0d0
3566 write (iout,*) "chain",ichain,chain_border1(1,ichain)-1,&
3567 chain_border(2,ichain-1)
3568 if (molnum(chain_border1(1,ichain)+1).eq.5) cycle
3570 d_t(:,chain_border1(1,ichain)-1)=&
3571 d_t(:,chain_border1(1,ichain)-1)+d_t(:,chain_border(2,ichain-1))
3572 d_t(:,chain_border(2,ichain-1))=0.0d0
3576 write (iout,*) "Random vel after 2nd transf the Calpha,SC space"
3577 write (iout,'(3hORG,1h(,i5,1h),3f10.5)') 0,(d_t(j,0),j=1,3)
3580 write (iout,'(a3,1h(,i5,1h),3f10.5,3x,3f10.5)')&
3581 restyp(itype(i,mnum),mnum),i,(d_t(j,i),j=1,3),(d_t(j,i+nres),j=1,3)
3588 !c d_t(j,0)=d_t(j,nnt)
3591 innt=chain_border(1,ichain)
3592 inct=chain_border(2,ichain)
3593 !c write (iout,*) "ichain",ichain," innt",innt," inct",inct
3594 !c write (iout,*) "ibeg",ibeg
3596 d_t(j,ibeg)=d_t(j,innt)
3601 if (iabs(itype(i,1).eq.10).or.mnum.ge.3) then
3602 !c write (iout,*) "i",i,(d_t(j,i),j=1,3),(d_t(j,i+1),j=1,3)
3604 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
3608 d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
3609 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
3618 "Random velocities in the virtual-bond-vector space"
3619 write (iout,'(3hORG,1h(,i5,1h),3f10.5)') 0,(d_t(j,0),j=1,3)
3621 write (iout,'(a3,1h(,i5,1h),3f10.5,3x,3f10.5)')&
3622 restyp(itype(i,mnum),mnum),i,(d_t(j,i),j=1,3),(d_t(j,i+nres),j=1,3)
3625 write (iout,*) "Kinetic energy from inertia matrix eigenvalues",&
3628 "Kinetic temperatures from inertia matrix eigenvalues",&
3631 write (iout,*) "Kinetic energy from inertia matrix",Ek3
3632 write (iout,*) "Kinetic temperatures from inertia",&
3635 write (iout,*) "Kinetic energy from velocities in CA-SC space",&
3638 "Kinetic temperatures from velovities in CA-SC space",&
3642 "Kinetic energy from virtual-bond-vector velocities",Ek1
3644 "Kinetic temperature from virtual-bond-vector velocities ",&
3653 sigv=dsqrt((Rb*t_bath)/geigen(i))
3656 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
3658 write (iout,*) "i",i," ii",ii," geigen",geigen(i),&
3659 " d_t_work_new",d_t_work_new(ii)
3670 Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
3673 write (iout,*) "Ek from eigenvectors",Ek1
3674 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
3683 d_t_work(ind)=d_t_work(ind) &
3684 +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
3686 ! write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
3690 ! Transfer to the d_t vector
3692 d_t(j,0)=d_t_work(j)
3698 d_t(j,i)=d_t_work(ind)
3703 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3704 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3705 .and.(mnum.lt.4)) then
3708 d_t(j,i+nres)=d_t_work(ind)
3714 ! write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",&
3715 ! 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
3717 ! write(iout,*) "end init MD"
3720 end subroutine random_vel
3721 !-----------------------------------------------------------------------------
3723 subroutine sd_verlet_p_setup
3724 ! Sets up the parameters of stochastic Verlet algorithm
3725 ! implicit real*8 (a-h,o-z)
3726 ! include 'DIMENSIONS'
3727 use control, only: tcpu
3732 ! include 'COMMON.CONTROL'
3733 ! include 'COMMON.VAR'
3734 ! include 'COMMON.MD'
3736 ! include 'COMMON.LANGEVIN'
3738 ! include 'COMMON.LANGEVIN.lang0'
3740 ! include 'COMMON.CHAIN'
3741 ! include 'COMMON.DERIV'
3742 ! include 'COMMON.GEO'
3743 ! include 'COMMON.LOCAL'
3744 ! include 'COMMON.INTERACT'
3745 ! include 'COMMON.IOUNITS'
3746 ! include 'COMMON.NAMES'
3747 ! include 'COMMON.TIME1'
3748 real(kind=8),dimension(6*nres) :: emgdt !(MAXRES6) maxres6=6*maxres
3749 real(kind=8) :: pterm,vterm,rho,rhoc,vsig
3750 real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
3751 prand_vec,vrand_vec1,vrand_vec2 !(MAXRES6) maxres6=6*maxres
3752 logical :: lprn = .false.
3753 real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
3754 real(kind=8) :: ktm,gdt,egdt,gdt2,gdt3,gdt4,gdt5,gdt6,gdt7,gdt8,&
3756 integer :: i,maxres2
3763 ! AL 8/17/04 Code adapted from tinker
3765 ! Get the frictional and random terms for stochastic dynamics in the
3766 ! eigenspace of mass-scaled UNRES friction matrix
3770 gdt = fricgam(i) * d_time
3772 ! Stochastic dynamics reduces to simple MD for zero friction
3774 if (gdt .le. zero) then
3775 pfric_vec(i) = 1.0d0
3776 vfric_vec(i) = d_time
3777 afric_vec(i) = 0.5d0 * d_time * d_time
3778 prand_vec(i) = 0.0d0
3779 vrand_vec1(i) = 0.0d0
3780 vrand_vec2(i) = 0.0d0
3782 ! Analytical expressions when friction coefficient is large
3785 if (gdt .ge. gdt_radius) then
3788 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
3789 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
3790 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
3791 vterm = 1.0d0 - egdt**2
3792 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
3794 ! Use series expansions when friction coefficient is small
3805 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0 &
3806 - gdt5/120.0d0 + gdt6/720.0d0 &
3807 - gdt7/5040.0d0 + gdt8/40320.0d0 &
3808 - gdt9/362880.0d0) / fricgam(i)**2
3809 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
3810 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
3811 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0 &
3812 + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0 &
3813 + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0 &
3814 + 127.0d0*gdt9/90720.0d0
3815 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0 &
3816 - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0 &
3817 - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0 &
3818 - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
3819 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0 &
3820 - 17.0d0*gdt2/1280.0d0 &
3821 + 17.0d0*gdt3/6144.0d0 &
3822 + 40967.0d0*gdt4/34406400.0d0 &
3823 - 57203.0d0*gdt5/275251200.0d0 &
3824 - 1429487.0d0*gdt6/13212057600.0d0)
3827 ! Compute the scaling factors of random terms for the nonzero friction case
3829 ktm = 0.5d0*d_time/fricgam(i)
3830 psig = dsqrt(ktm*pterm) / fricgam(i)
3831 vsig = dsqrt(ktm*vterm)
3832 rhoc = dsqrt(1.0d0 - rho*rho)
3834 vrand_vec1(i) = vsig * rho
3835 vrand_vec2(i) = vsig * rhoc
3840 "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",&
3843 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),&
3844 afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
3848 ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
3851 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
3852 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
3853 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
3854 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
3855 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
3856 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
3859 t_sdsetup=t_sdsetup+MPI_Wtime()
3861 t_sdsetup=t_sdsetup+tcpu()-tt0
3864 end subroutine sd_verlet_p_setup
3865 !-----------------------------------------------------------------------------
3866 subroutine eigtransf1(n,ndim,ab,d,c)
3870 real(kind=8) :: ab(ndim,ndim,n),c(ndim,n),d(ndim)
3876 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
3881 end subroutine eigtransf1
3882 !-----------------------------------------------------------------------------
3883 subroutine eigtransf(n,ndim,a,b,d,c)
3887 real(kind=8) :: a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
3893 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
3898 end subroutine eigtransf
3899 !-----------------------------------------------------------------------------
3900 subroutine sd_verlet1
3902 ! Applying stochastic velocity Verlet algorithm - step 1 to velocities
3904 ! implicit real*8 (a-h,o-z)
3905 ! include 'DIMENSIONS'
3906 ! include 'COMMON.CONTROL'
3907 ! include 'COMMON.VAR'
3908 ! include 'COMMON.MD'
3910 ! include 'COMMON.LANGEVIN'
3912 ! include 'COMMON.LANGEVIN.lang0'
3914 ! include 'COMMON.CHAIN'
3915 ! include 'COMMON.DERIV'
3916 ! include 'COMMON.GEO'
3917 ! include 'COMMON.LOCAL'
3918 ! include 'COMMON.INTERACT'
3919 ! include 'COMMON.IOUNITS'
3920 ! include 'COMMON.NAMES'
3921 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
3922 !el common /stochcalc/ stochforcvec
3923 logical :: lprn = .false.
3924 real(kind=8) :: ddt1,ddt2
3925 integer :: i,j,ind,inres
3927 ! write (iout,*) "dc_old"
3929 ! write (iout,'(i5,3f10.5,5x,3f10.5)')
3930 ! & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3933 dc_work(j)=dc_old(j,0)
3934 d_t_work(j)=d_t_old(j,0)
3935 d_a_work(j)=d_a_old(j,0)
3940 dc_work(ind+j)=dc_old(j,i)
3941 d_t_work(ind+j)=d_t_old(j,i)
3942 d_a_work(ind+j)=d_a_old(j,i)
3948 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3949 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3950 .and.(mnum.lt.4)) then
3952 dc_work(ind+j)=dc_old(j,i+nres)
3953 d_t_work(ind+j)=d_t_old(j,i+nres)
3954 d_a_work(ind+j)=d_a_old(j,i+nres)
3962 "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",&
3966 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
3967 vfric_mat(i,j),afric_mat(i,j),&
3968 prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3976 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) &
3977 +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3978 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3979 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3981 d_t_work_new(i)=ddt1+0.5d0*ddt2
3982 d_t_work(i)=ddt1+ddt2
3987 d_t(j,0)=d_t_work(j)
3992 dc(j,i)=dc_work(ind+j)
3993 d_t(j,i)=d_t_work(ind+j)
3999 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
4000 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4001 .and.(mnum.lt.4)) then
4004 dc(j,inres)=dc_work(ind+j)
4005 d_t(j,inres)=d_t_work(ind+j)
4011 end subroutine sd_verlet1
4012 !-----------------------------------------------------------------------------
4013 subroutine sd_verlet2
4015 ! Calculating the adjusted velocities for accelerations
4017 ! implicit real*8 (a-h,o-z)
4018 ! include 'DIMENSIONS'
4019 ! include 'COMMON.CONTROL'
4020 ! include 'COMMON.VAR'
4021 ! include 'COMMON.MD'
4023 ! include 'COMMON.LANGEVIN'
4025 ! include 'COMMON.LANGEVIN.lang0'
4027 ! include 'COMMON.CHAIN'
4028 ! include 'COMMON.DERIV'
4029 ! include 'COMMON.GEO'
4030 ! include 'COMMON.LOCAL'
4031 ! include 'COMMON.INTERACT'
4032 ! include 'COMMON.IOUNITS'
4033 ! include 'COMMON.NAMES'
4034 !el real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV !(MAXRES6) maxres6=6*maxres
4035 real(kind=8),dimension(6*nres) :: stochforcvecV !(MAXRES6) maxres6=6*maxres
4036 !el common /stochcalc/ stochforcvec
4038 real(kind=8) :: ddt1,ddt2
4039 integer :: i,j,ind,inres
4040 ! Compute the stochastic forces which contribute to velocity change
4042 call stochastic_force(stochforcvecV)
4049 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
4050 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+ &
4051 vrand_mat2(i,j)*stochforcvecV(j)
4053 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
4057 d_t(j,0)=d_t_work(j)
4062 d_t(j,i)=d_t_work(ind+j)
4068 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
4069 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4070 .and.(mnum.lt.4)) then
4073 d_t(j,inres)=d_t_work(ind+j)
4079 end subroutine sd_verlet2
4080 !-----------------------------------------------------------------------------
4081 subroutine sd_verlet_ciccotti_setup
4083 ! Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
4085 ! implicit real*8 (a-h,o-z)
4086 ! include 'DIMENSIONS'
4087 use control, only: tcpu
4092 ! include 'COMMON.CONTROL'
4093 ! include 'COMMON.VAR'
4094 ! include 'COMMON.MD'
4096 ! include 'COMMON.LANGEVIN'
4098 ! include 'COMMON.LANGEVIN.lang0'
4100 ! include 'COMMON.CHAIN'
4101 ! include 'COMMON.DERIV'
4102 ! include 'COMMON.GEO'
4103 ! include 'COMMON.LOCAL'
4104 ! include 'COMMON.INTERACT'
4105 ! include 'COMMON.IOUNITS'
4106 ! include 'COMMON.NAMES'
4107 ! include 'COMMON.TIME1'
4108 real(kind=8),dimension(6*nres) :: emgdt !(MAXRES6) maxres6=6*maxres
4109 real(kind=8) :: pterm,vterm,rho,rhoc,vsig
4110 real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
4111 prand_vec,vrand_vec1,vrand_vec2 !(MAXRES6) maxres6=6*maxres
4112 logical :: lprn = .false.
4113 real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
4114 real(kind=8) :: ktm,gdt,egdt,tt0
4115 integer :: i,maxres2
4122 ! AL 8/17/04 Code adapted from tinker
4124 ! Get the frictional and random terms for stochastic dynamics in the
4125 ! eigenspace of mass-scaled UNRES friction matrix
4129 write (iout,*) "i",i," fricgam",fricgam(i)
4130 gdt = fricgam(i) * d_time
4132 ! Stochastic dynamics reduces to simple MD for zero friction
4134 if (gdt .le. zero) then
4135 pfric_vec(i) = 1.0d0
4136 vfric_vec(i) = d_time
4137 afric_vec(i) = 0.5d0*d_time*d_time
4138 prand_vec(i) = afric_vec(i)
4139 vrand_vec2(i) = vfric_vec(i)
4141 ! Analytical expressions when friction coefficient is large
4146 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
4147 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
4148 prand_vec(i) = afric_vec(i)
4149 vrand_vec2(i) = vfric_vec(i)
4151 ! Compute the scaling factors of random terms for the nonzero friction case
4153 ! ktm = 0.5d0*d_time/fricgam(i)
4154 ! psig = dsqrt(ktm*pterm) / fricgam(i)
4155 ! vsig = dsqrt(ktm*vterm)
4156 ! prand_vec(i) = psig*afric_vec(i)
4157 ! vrand_vec2(i) = vsig*vfric_vec(i)
4162 "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",&
4165 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),&
4166 afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
4170 ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
4172 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
4173 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
4174 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
4175 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
4176 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
4178 t_sdsetup=t_sdsetup+MPI_Wtime()
4180 t_sdsetup=t_sdsetup+tcpu()-tt0
4183 end subroutine sd_verlet_ciccotti_setup
4184 !-----------------------------------------------------------------------------
4185 subroutine sd_verlet1_ciccotti
4187 ! Applying stochastic velocity Verlet algorithm - step 1 to velocities
4188 ! implicit real*8 (a-h,o-z)
4190 ! include 'DIMENSIONS'
4194 ! include 'COMMON.CONTROL'
4195 ! include 'COMMON.VAR'
4196 ! include 'COMMON.MD'
4198 ! include 'COMMON.LANGEVIN'
4200 ! include 'COMMON.LANGEVIN.lang0'
4202 ! include 'COMMON.CHAIN'
4203 ! include 'COMMON.DERIV'
4204 ! include 'COMMON.GEO'
4205 ! include 'COMMON.LOCAL'
4206 ! include 'COMMON.INTERACT'
4207 ! include 'COMMON.IOUNITS'
4208 ! include 'COMMON.NAMES'
4209 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
4210 !el common /stochcalc/ stochforcvec
4211 logical :: lprn = .false.
4212 real(kind=8) :: ddt1,ddt2
4213 integer :: i,j,ind,inres
4214 ! write (iout,*) "dc_old"
4216 ! write (iout,'(i5,3f10.5,5x,3f10.5)')
4217 ! & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
4220 dc_work(j)=dc_old(j,0)
4221 d_t_work(j)=d_t_old(j,0)
4222 d_a_work(j)=d_a_old(j,0)
4227 dc_work(ind+j)=dc_old(j,i)
4228 d_t_work(ind+j)=d_t_old(j,i)
4229 d_a_work(ind+j)=d_a_old(j,i)
4234 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
4236 dc_work(ind+j)=dc_old(j,i+nres)
4237 d_t_work(ind+j)=d_t_old(j,i+nres)
4238 d_a_work(ind+j)=d_a_old(j,i+nres)
4247 "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",&
4251 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
4252 vfric_mat(i,j),afric_mat(i,j),&
4253 prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
4261 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) &
4262 +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
4263 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
4264 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
4266 d_t_work_new(i)=ddt1+0.5d0*ddt2
4267 d_t_work(i)=ddt1+ddt2
4272 d_t(j,0)=d_t_work(j)
4277 dc(j,i)=dc_work(ind+j)
4278 d_t(j,i)=d_t_work(ind+j)
4284 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
4285 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4286 .and.(mnum.lt.4)) then
4289 dc(j,inres)=dc_work(ind+j)
4290 d_t(j,inres)=d_t_work(ind+j)
4296 end subroutine sd_verlet1_ciccotti
4297 !-----------------------------------------------------------------------------
4298 subroutine sd_verlet2_ciccotti
4300 ! Calculating the adjusted velocities for accelerations
4302 ! implicit real*8 (a-h,o-z)
4303 ! include 'DIMENSIONS'
4304 ! include 'COMMON.CONTROL'
4305 ! include 'COMMON.VAR'
4306 ! include 'COMMON.MD'
4308 ! include 'COMMON.LANGEVIN'
4310 ! include 'COMMON.LANGEVIN.lang0'
4312 ! include 'COMMON.CHAIN'
4313 ! include 'COMMON.DERIV'
4314 ! include 'COMMON.GEO'
4315 ! include 'COMMON.LOCAL'
4316 ! include 'COMMON.INTERACT'
4317 ! include 'COMMON.IOUNITS'
4318 ! include 'COMMON.NAMES'
4319 !el real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV !(MAXRES6) maxres6=6*maxres
4320 real(kind=8),dimension(6*nres) :: stochforcvecV !(MAXRES6) maxres6=6*maxres
4321 !el common /stochcalc/ stochforcvec
4322 real(kind=8) :: ddt1,ddt2
4323 integer :: i,j,ind,inres
4325 ! Compute the stochastic forces which contribute to velocity change
4327 call stochastic_force(stochforcvecV)
4334 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
4335 ! ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
4336 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
4338 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
4342 d_t(j,0)=d_t_work(j)
4347 d_t(j,i)=d_t_work(ind+j)
4353 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4355 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
4358 d_t(j,inres)=d_t_work(ind+j)
4364 end subroutine sd_verlet2_ciccotti
4366 !-----------------------------------------------------------------------------
4368 !-----------------------------------------------------------------------------
4370 subroutine inertia_tensor
4373 real(kind=8) Im(3,3),Imcp(3,3),pr(3),M_SC,&
4374 eigvec(3,3),Id(3,3),eigval(3),L(3),vp(3),vrot(3),&
4375 vpp(3,0:MAXRES),vs_p(3),pr1(3,3),&
4376 pr2(3,3),pp(3),incr(3),v(3),mag,mag2,M_PEP
4377 integer iti,inres,i,j,k,mnum,mnum1
4390 !c caulating the center of the mass of the protein
4394 if (itype(i,mnum).eq.ntyp1_molec(mnum)&
4395 .or. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) cycle
4396 ! if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
4397 if (mnum.ge.5) mp(mnum)=0.0d0
4398 M_PEP=M_PEP+mp(mnum)
4401 cm(j)=cm(j)+(c(j,i)+0.5d0*dc(j,i))*mp(mnum)
4411 iti=iabs(itype(i,mnum))
4412 if (iti.eq.ntyp1_molec(mnum)) cycle
4413 M_SC=M_SC+msc(iabs(iti),mnum)
4416 cm(j)=cm(j)+msc(iabs(iti),mnum)*c(j,inres)
4420 cm(j)=cm(j)/(M_SC+M_PEP)
4425 ! if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
4426 if (mnum.ge.5) mp(mnum)=0.0d0
4427 if (itype(i,mnum).eq.ntyp1_molec(mnum)&
4428 .or. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) cycle
4430 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
4432 Im(1,1)=Im(1,1)+mp(mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
4433 Im(1,2)=Im(1,2)-mp(mnum)*pr(1)*pr(2)
4434 Im(1,3)=Im(1,3)-mp(mnum)*pr(1)*pr(3)
4435 Im(2,3)=Im(2,3)-mp(mnum)*pr(2)*pr(3)
4436 Im(2,2)=Im(2,2)+mp(mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
4437 Im(3,3)=Im(3,3)+mp(mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
4442 iti=iabs(itype(i,mnum))
4443 if (iti.eq.ntyp1_molec(mnum)) cycle
4446 pr(j)=c(j,inres)-cm(j)
4448 Im(1,1)=Im(1,1)+msc(iabs(iti),mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
4449 Im(1,2)=Im(1,2)-msc(iabs(iti),mnum)*pr(1)*pr(2)
4450 Im(1,3)=Im(1,3)-msc(iabs(iti),mnum)*pr(1)*pr(3)
4451 Im(2,3)=Im(2,3)-msc(iabs(iti),mnum)*pr(2)*pr(3)
4452 Im(2,2)=Im(2,2)+msc(iabs(iti),mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
4453 Im(3,3)=Im(3,3)+msc(iabs(iti),mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
4458 if (itype(i,mnum).eq.ntyp1_molec(mnum)&
4459 .or. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) cycle
4460 Im(1,1)=Im(1,1)+Ip(mnum)*(1-dc_norm(1,i)*dc_norm(1,i))*&
4461 vbld(i+1)*vbld(i+1)*0.25d0
4462 Im(1,2)=Im(1,2)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(2,i))*&
4463 vbld(i+1)*vbld(i+1)*0.25d0
4464 Im(1,3)=Im(1,3)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(3,i))*&
4465 vbld(i+1)*vbld(i+1)*0.25d0
4466 Im(2,3)=Im(2,3)+Ip(mnum)*(-dc_norm(2,i)*dc_norm(3,i))*&
4467 vbld(i+1)*vbld(i+1)*0.25d0
4468 Im(2,2)=Im(2,2)+Ip(mnum)*(1-dc_norm(2,i)*dc_norm(2,i))*&
4469 vbld(i+1)*vbld(i+1)*0.25d0
4470 Im(3,3)=Im(3,3)+Ip(mnum)*(1-dc_norm(3,i)*dc_norm(3,i))*&
4471 vbld(i+1)*vbld(i+1)*0.25d0
4476 iti=iabs(itype(i,mnum))
4477 if (iti.ne.10 .and. iti.ne.ntyp1_molec(mnum).and.mnum.le.2) then
4479 Im(1,1)=Im(1,1)+Isc(iti,mnum)*(1-dc_norm(1,inres)*&
4480 dc_norm(1,inres))*vbld(inres)*vbld(inres)
4481 Im(1,2)=Im(1,2)-Isc(iti,mnum)*(dc_norm(1,inres)*&
4482 dc_norm(2,inres))*vbld(inres)*vbld(inres)
4483 Im(1,3)=Im(1,3)-Isc(iti,mnum)*(dc_norm(1,inres)*&
4484 dc_norm(3,inres))*vbld(inres)*vbld(inres)
4485 Im(2,3)=Im(2,3)-Isc(iti,mnum)*(dc_norm(2,inres)*&
4486 dc_norm(3,inres))*vbld(inres)*vbld(inres)
4487 Im(2,2)=Im(2,2)+Isc(iti,mnum)*(1-dc_norm(2,inres)*&
4488 dc_norm(2,inres))*vbld(inres)*vbld(inres)
4489 Im(3,3)=Im(3,3)+Isc(iti,mnum)*(1-dc_norm(3,inres)*&
4490 dc_norm(3,inres))*vbld(inres)*vbld(inres)
4499 !c Copng the Im matrix for the djacob subroutine
4506 !c Finding the eigenvectors and eignvalues of the inertia tensor
4507 call djacob(3,3,10000,1.0d-10,Imcp,eigvec,eigval)
4509 if (dabs(eigval(i)).gt.1.0d-15) then
4510 Id(i,i)=1.0d0/eigval(i)
4517 Imcp(i,j)=eigvec(j,i)
4523 pr1(i,j)=pr1(i,j)+Id(i,k)*Imcp(k,j)
4530 pr2(i,j)=pr2(i,j)+eigvec(i,k)*pr1(k,j)
4534 !c Calculating the total rotational velocity of the molecule
4537 vrot(i)=vrot(i)+pr2(i,j)*L(j)
4543 ! write (iout,*) itype(i+1,mnum1),itype(i,mnum)
4544 if (itype(i+1,mnum1).ne.ntyp1_molec(mnum1) &
4545 .and. itype(i,mnum).eq.ntyp1_molec(mnum) .or.&
4546 itype(i,mnum).ne.ntyp1_molec(mnum) &
4547 .and. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) cycle
4548 call vecpr(vrot(1),dc(1,i),vp)
4550 d_t(j,i)=d_t(j,i)-vp(j)
4555 if(itype(i,1).ne.10 .and.itype(i,mnum).ne.ntyp1_molec(mnum)&
4556 .and.mnum.le.2) then
4558 call vecpr(vrot(1),dc(1,inres),vp)
4560 d_t(j,inres)=d_t(j,inres)-vp(j)
4566 end subroutine inertia_tensor
4567 !------------------------------------------------------------
4568 subroutine angmom(cm,L)
4570 double precision L(3),cm(3),pr(3),vp(3),vrot(3),incr(3),v(3),&
4572 integer iti,inres,i,j,mnum,mnum1
4573 !c Calculate the angular momentum
4583 ! if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
4584 if (mnum.ge.5) mp(mnum)=0.0d0
4585 if (itype(i,mnum).eq.ntyp1_molec(mnum)&
4586 .or. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) cycle
4588 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
4591 v(j)=incr(j)+0.5d0*d_t(j,i)
4594 incr(j)=incr(j)+d_t(j,i)
4596 call vecpr(pr(1),v(1),vp)
4598 L(j)=L(j)+mp(mnum)*vp(j)
4602 pp(j)=0.5d0*d_t(j,i)
4604 call vecpr(pr(1),pp(1),vp)
4606 L(j)=L(j)+Ip(mnum)*vp(j)
4614 iti=iabs(itype(i,mnum))
4617 pr(j)=c(j,inres)-cm(j)
4619 if (itype(i,1).ne.10 .and.itype(i,mnum).ne.ntyp1_molec(mnum)&
4620 .and.mnum.le.2) then
4622 v(j)=incr(j)+d_t(j,inres)
4629 call vecpr(pr(1),v(1),vp)
4630 !c write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),
4631 !c & " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
4632 ! if (mnum.gt.4) then
4638 L(j)=L(j)+mscab*vp(j)
4640 !c write (iout,*) "L",(l(j),j=1,3)
4641 if (itype(i,1).ne.10 .and.itype(i,mnum).ne.ntyp1_molec(mnum)&
4642 .and.mnum.le.2) then
4644 v(j)=incr(j)+d_t(j,inres)
4646 call vecpr(dc(1,inres),d_t(1,inres),vp)
4648 L(j)=L(j)+Isc(iti,mnum)*vp(j)
4652 incr(j)=incr(j)+d_t(j,i)
4656 end subroutine angmom
4657 !---------------------------------------------------------------
4658 subroutine vcm_vel(vcm)
4659 double precision vcm(3),vv(3),summas,amas
4660 integer i,j,mnum,mnum1
4668 if ((mnum.ge.5).or.(mnum.eq.3))&
4669 mp(mnum)=msc(itype(i,mnum),mnum)
4671 summas=summas+mp(mnum)
4673 vcm(j)=vcm(j)+mp(mnum)*(vv(j)+0.5d0*d_t(j,i))
4677 amas=msc(iabs(itype(i,mnum)),mnum)
4681 ! amas=msc(iabs(itype(i)))
4683 ! if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
4684 if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4685 .and.(mnum.lt.4)) then
4687 vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres))
4691 vcm(j)=vcm(j)+amas*vv(j)
4695 vv(j)=vv(j)+d_t(j,i)
4698 !c write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas
4700 vcm(j)=vcm(j)/summas
4703 end subroutine vcm_vel
4705 subroutine inertia_tensor
4707 ! Calculating the intertia tensor for the entire protein in order to
4708 ! remove the perpendicular components of velocity matrix which cause
4709 ! the molecule to rotate.
4712 ! implicit real*8 (a-h,o-z)
4713 ! include 'DIMENSIONS'
4714 ! include 'COMMON.CONTROL'
4715 ! include 'COMMON.VAR'
4716 ! include 'COMMON.MD'
4717 ! include 'COMMON.CHAIN'
4718 ! include 'COMMON.DERIV'
4719 ! include 'COMMON.GEO'
4720 ! include 'COMMON.LOCAL'
4721 ! include 'COMMON.INTERACT'
4722 ! include 'COMMON.IOUNITS'
4723 ! include 'COMMON.NAMES'
4725 real(kind=8),dimension(3,3) :: Im,Imcp,eigvec,Id
4726 real(kind=8),dimension(3) :: pr,eigval,L,vp,vrot
4727 real(kind=8) :: M_SC,mag,mag2,M_PEP
4728 real(kind=8),dimension(3,0:nres) :: vpp !(3,0:MAXRES)
4729 real(kind=8),dimension(3) :: vs_p,pp,incr,v
4730 real(kind=8),dimension(3,3) :: pr1,pr2
4732 !el common /gucio/ cm
4733 integer :: iti,inres,i,j,k,mnum
4744 ! calculating the center of the mass of the protein
4748 ! if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
4749 write(iout,*) "WTF",itype(i,mnum),i,mnum,mp(mnum)
4750 ! if (itype(i,mnum).eq.ntyp1_molec(mnum)) cycle
4751 M_PEP=M_PEP+mp(mnum)
4754 cm(j)=cm(j)+(c(j,i)+0.5d0*dc(j,i))*mp(mnum)
4763 ! if (mnum.ge.5) cycle
4764 iti=iabs(itype(i,mnum))
4765 M_SC=M_SC+msc(iabs(iti),mnum)
4767 if (mnum.ge.4) inres=i
4769 cm(j)=cm(j)+msc(iabs(iti),mnum)*c(j,inres)
4773 cm(j)=cm(j)/(M_SC+M_PEP)
4775 ! write(iout,*) "Center of mass:",cm
4778 ! if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
4780 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
4782 Im(1,1)=Im(1,1)+mp(mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
4783 Im(1,2)=Im(1,2)-mp(mnum)*pr(1)*pr(2)
4784 Im(1,3)=Im(1,3)-mp(mnum)*pr(1)*pr(3)
4785 Im(2,3)=Im(2,3)-mp(mnum)*pr(2)*pr(3)
4786 Im(2,2)=Im(2,2)+mp(mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
4787 Im(3,3)=Im(3,3)+mp(mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
4790 ! write(iout,*) "The angular momentum before msc add"
4792 ! write (iout,*) (Im(i,j),j=1,3)
4797 iti=iabs(itype(i,mnum))
4798 ! if (mnum.ge.5) cycle
4800 if (mnum.ge.4) inres=i
4802 pr(j)=c(j,inres)-cm(j)
4804 Im(1,1)=Im(1,1)+msc(iabs(iti),mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
4805 Im(1,2)=Im(1,2)-msc(iabs(iti),mnum)*pr(1)*pr(2)
4806 Im(1,3)=Im(1,3)-msc(iabs(iti),mnum)*pr(1)*pr(3)
4807 Im(2,3)=Im(2,3)-msc(iabs(iti),mnum)*pr(2)*pr(3)
4808 Im(2,2)=Im(2,2)+msc(iabs(iti),mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
4809 Im(3,3)=Im(3,3)+msc(iabs(iti),mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
4811 ! write(iout,*) "The angular momentum before Ip add"
4813 ! write (iout,*) (Im(i,j),j=1,3)
4818 Im(1,1)=Im(1,1)+Ip(mnum)*(1-dc_norm(1,i)*dc_norm(1,i))* &
4819 vbld(i+1)*vbld(i+1)*0.25d0
4820 Im(1,2)=Im(1,2)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(2,i))* &
4821 vbld(i+1)*vbld(i+1)*0.25d0
4822 Im(1,3)=Im(1,3)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(3,i))* &
4823 vbld(i+1)*vbld(i+1)*0.25d0
4824 Im(2,3)=Im(2,3)+Ip(mnum)*(-dc_norm(2,i)*dc_norm(3,i))* &
4825 vbld(i+1)*vbld(i+1)*0.25d0
4826 Im(2,2)=Im(2,2)+Ip(mnum)*(1-dc_norm(2,i)*dc_norm(2,i))* &
4827 vbld(i+1)*vbld(i+1)*0.25d0
4828 Im(3,3)=Im(3,3)+Ip(mnum)*(1-dc_norm(3,i)*dc_norm(3,i))* &
4829 vbld(i+1)*vbld(i+1)*0.25d0
4831 ! write(iout,*) "The angular momentum before Isc add"
4833 ! write (iout,*) (Im(i,j),j=1,3)
4839 ! if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
4840 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4841 .and.(mnum.lt.4)) then
4842 iti=iabs(itype(i,mnum))
4844 Im(1,1)=Im(1,1)+Isc(iti,mnum)*(1-dc_norm(1,inres)* &
4845 dc_norm(1,inres))*vbld(inres)*vbld(inres)
4846 Im(1,2)=Im(1,2)-Isc(iti,mnum)*(dc_norm(1,inres)* &
4847 dc_norm(2,inres))*vbld(inres)*vbld(inres)
4848 Im(1,3)=Im(1,3)-Isc(iti,mnum)*(dc_norm(1,inres)* &
4849 dc_norm(3,inres))*vbld(inres)*vbld(inres)
4850 Im(2,3)=Im(2,3)-Isc(iti,mnum)*(dc_norm(2,inres)* &
4851 dc_norm(3,inres))*vbld(inres)*vbld(inres)
4852 Im(2,2)=Im(2,2)+Isc(iti,mnum)*(1-dc_norm(2,inres)* &
4853 dc_norm(2,inres))*vbld(inres)*vbld(inres)
4854 Im(3,3)=Im(3,3)+Isc(iti,mnum)*(1-dc_norm(3,inres)* &
4855 dc_norm(3,inres))*vbld(inres)*vbld(inres)
4859 ! write(iout,*) "The angular momentum before agnom:"
4861 ! write (iout,*) (Im(i,j),j=1,3)
4865 ! write(iout,*) "The angular momentum before adjustment:"
4866 ! write(iout,*) (L(j),j=1,3)
4868 ! write (iout,*) (Im(i,j),j=1,3)
4874 ! Copying the Im matrix for the djacob subroutine
4882 ! Finding the eigenvectors and eignvalues of the inertia tensor
4883 call djacob(3,3,10000,1.0d-10,Imcp,eigvec,eigval)
4884 ! write (iout,*) "Eigenvalues & Eigenvectors"
4885 ! write (iout,'(5x,3f10.5)') (eigval(i),i=1,3)
4888 ! write (iout,'(i5,3f10.5)') i,(eigvec(i,j),j=1,3)
4890 ! Constructing the diagonalized matrix
4892 if (dabs(eigval(i)).gt.1.0d-15) then
4893 Id(i,i)=1.0d0/eigval(i)
4900 Imcp(i,j)=eigvec(j,i)
4906 pr1(i,j)=pr1(i,j)+Id(i,k)*Imcp(k,j)
4913 pr2(i,j)=pr2(i,j)+eigvec(i,k)*pr1(k,j)
4917 ! Calculating the total rotational velocity of the molecule
4920 vrot(i)=vrot(i)+pr2(i,j)*L(j)
4923 ! Resetting the velocities
4925 call vecpr(vrot(1),dc(1,i),vp)
4927 ! print *,"HERE2",d_t(j,i),vp(j)
4928 d_t(j,i)=d_t(j,i)-vp(j)
4929 ! print *,"HERE2",d_t(j,i),vp(j)
4934 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4935 .and.(mnum.lt.4)) then
4936 ! if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
4937 ! if(itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
4939 call vecpr(vrot(1),dc(1,inres),vp)
4941 d_t(j,inres)=d_t(j,inres)-vp(j)
4946 ! write(iout,*) "The angular momentum after adjustment:"
4947 ! write(iout,*) (L(j),j=1,3)
4950 end subroutine inertia_tensor
4951 !-----------------------------------------------------------------------------
4952 subroutine angmom(cm,L)
4955 ! implicit real*8 (a-h,o-z)
4956 ! include 'DIMENSIONS'
4957 ! include 'COMMON.CONTROL'
4958 ! include 'COMMON.VAR'
4959 ! include 'COMMON.MD'
4960 ! include 'COMMON.CHAIN'
4961 ! include 'COMMON.DERIV'
4962 ! include 'COMMON.GEO'
4963 ! include 'COMMON.LOCAL'
4964 ! include 'COMMON.INTERACT'
4965 ! include 'COMMON.IOUNITS'
4966 ! include 'COMMON.NAMES'
4967 real(kind=8) :: mscab
4968 real(kind=8),dimension(3) :: L,cm,pr,vp,vrot,incr,v,pp
4969 integer :: iti,inres,i,j,mnum
4970 ! Calculate the angular momentum
4979 if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
4981 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
4984 v(j)=incr(j)+0.5d0*d_t(j,i)
4987 incr(j)=incr(j)+d_t(j,i)
4989 call vecpr(pr(1),v(1),vp)
4991 L(j)=L(j)+mp(mnum)*vp(j)
4992 ! print *,"HERE3",J,i,L(j),mp(mnum),Ip(mnum),mnum
4996 pp(j)=0.5d0*d_t(j,i)
4998 call vecpr(pr(1),pp(1),vp)
5001 L(j)=L(j)+Ip(mnum)*vp(j)
5009 iti=iabs(itype(i,mnum))
5017 pr(j)=c(j,inres)-cm(j)
5020 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5021 .and.(mnum.lt.4)) then
5023 v(j)=incr(j)+d_t(j,inres)
5030 ! print *,i,pr,"pr",v
5031 call vecpr(pr(1),v(1),vp)
5032 ! write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),&
5033 ! " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
5035 L(j)=L(j)+mscab*vp(j)
5037 ! write (iout,*) "L",(l(j),j=1,3)
5038 ! print *,"L",(l(j),j=1,3),i,vp(1)
5040 if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5041 .and.(mnum.lt.4)) then
5043 v(j)=incr(j)+d_t(j,inres)
5045 call vecpr(dc(1,inres),d_t(1,inres),vp)
5047 L(j)=L(j)+Isc(iti,mnum)*vp(j)
5051 incr(j)=incr(j)+d_t(j,i)
5055 end subroutine angmom
5056 !-----------------------------------------------------------------------------
5057 subroutine vcm_vel(vcm)
5060 ! implicit real*8 (a-h,o-z)
5061 ! include 'DIMENSIONS'
5062 ! include 'COMMON.VAR'
5063 ! include 'COMMON.MD'
5064 ! include 'COMMON.CHAIN'
5065 ! include 'COMMON.DERIV'
5066 ! include 'COMMON.GEO'
5067 ! include 'COMMON.LOCAL'
5068 ! include 'COMMON.INTERACT'
5069 ! include 'COMMON.IOUNITS'
5070 real(kind=8),dimension(3) :: vcm,vv
5071 real(kind=8) :: summas,amas
5081 ! if (mnum.ge.4) mp(mnum)=msc(itype(i,mnum),mnum)
5083 summas=summas+mp(mnum)
5085 vcm(j)=vcm(j)+mp(mnum)*(vv(j)+0.5d0*d_t(j,i))
5086 ! print *,i,j,vv(j),d_t(j,i)
5089 ! if (mnum.ne.4) then
5090 amas=msc(iabs(itype(i,mnum)),mnum)
5095 if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5096 .and.(mnum.lt.4)) then
5098 vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres))
5102 vcm(j)=vcm(j)+amas*vv(j)
5106 vv(j)=vv(j)+d_t(j,i)
5109 ! write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas
5111 vcm(j)=vcm(j)/summas
5114 end subroutine vcm_vel
5116 !-----------------------------------------------------------------------------
5118 !-----------------------------------------------------------------------------
5120 ! RATTLE algorithm for velocity Verlet - step 1, UNRES
5124 ! implicit real*8 (a-h,o-z)
5125 ! include 'DIMENSIONS'
5127 ! include 'COMMON.CONTROL'
5128 ! include 'COMMON.VAR'
5129 ! include 'COMMON.MD'
5131 ! include 'COMMON.LANGEVIN'
5133 ! include 'COMMON.LANGEVIN.lang0'
5135 ! include 'COMMON.CHAIN'
5136 ! include 'COMMON.DERIV'
5137 ! include 'COMMON.GEO'
5138 ! include 'COMMON.LOCAL'
5139 ! include 'COMMON.INTERACT'
5140 ! include 'COMMON.IOUNITS'
5141 ! include 'COMMON.NAMES'
5142 ! include 'COMMON.TIME1'
5143 !el real(kind=8) :: gginv(2*nres,2*nres),&
5144 !el gdc(3,2*nres,2*nres)
5145 real(kind=8) :: dC_uncor(3,2*nres) !,&
5146 !el real(kind=8) :: Cmat(2*nres,2*nres)
5147 real(kind=8) :: x(2*nres),xcorr(3,2*nres) !maxres2=2*maxres
5148 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
5149 !el common /przechowalnia/ nbond
5150 integer :: max_rattle = 5
5151 logical :: lprn = .false., lprn1 = .false., not_done
5152 real(kind=8) :: tol_rattle = 1.0d-5
5154 integer :: ii,i,j,jj,l,ind,ind1,nres2
5157 !el /common/ przechowalnia
5159 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
5160 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
5161 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
5163 if (lprn) write (iout,*) "RATTLE1"
5167 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5168 .and.(mnum.lt.4)) nbond=nbond+1
5170 ! Make a folded form of the Ginv-matrix
5183 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
5188 if (itype(k,1).ne.10 .and. itype(k,mnum).ne.ntyp1_molec(mnum)&
5189 .and.(mnum.lt.4)) then
5193 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
5201 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5212 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
5216 if (itype(k,1).ne.10) then
5220 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
5228 write (iout,*) "Matrix GGinv"
5229 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
5235 if (iter.gt.max_rattle) then
5236 write (iout,*) "Error - too many iterations in RATTLE."
5239 ! Calculate the matrix C = GG**(-1) dC_old o dC
5244 dC_uncor(j,ind1)=dC(j,i)
5248 if (itype(i,1).ne.10) then
5251 dC_uncor(j,ind1)=dC(j,i+nres)
5260 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
5264 if (itype(k,1).ne.10) then
5267 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
5272 ! Calculate deviations from standard virtual-bond lengths
5276 x(ind)=vbld(i+1)**2-vbl**2
5279 if (itype(i,1).ne.10) then
5281 x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
5285 write (iout,*) "Coordinates and violations"
5287 write(iout,'(i5,3f10.5,5x,e15.5)') &
5288 i,(dC_uncor(j,i),j=1,3),x(i)
5290 write (iout,*) "Velocities and violations"
5294 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5295 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
5299 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5300 .and.(mnum.lt.4)) then
5303 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5304 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
5305 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
5308 ! write (iout,*) "gdc"
5310 ! write (iout,*) "i",i
5312 ! write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
5318 if (dabs(x(i)).gt.xmax) then
5322 if (xmax.lt.tol_rattle) then
5326 ! Calculate the matrix of the system of equations
5331 Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
5336 write (iout,*) "Matrix Cmat"
5337 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
5339 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
5340 ! Add constraint term to positions
5347 xx = xx+x(ii)*gdc(j,ind,ii)
5351 d_t_new(j,i)=d_t_new(j,i)-xx/d_time
5355 if (itype(i,1).ne.10) then
5360 xx = xx+x(ii)*gdc(j,ind,ii)
5363 dC(j,i+nres)=dC(j,i+nres)-xx
5364 d_t_new(j,i+nres)=d_t_new(j,i+nres)-xx/d_time
5368 ! Rebuild the chain using the new coordinates
5369 call chainbuild_cart
5371 write (iout,*) "New coordinates, Lagrange multipliers,",&
5372 " and differences between actual and standard bond lengths"
5376 xx=vbld(i+1)**2-vbl**2
5377 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
5378 i,(dC(j,i),j=1,3),x(ind),xx
5382 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5385 xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
5386 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
5387 i,(dC(j,i+nres),j=1,3),x(ind),xx
5390 write (iout,*) "Velocities and violations"
5394 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5395 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
5398 if (itype(i,1).ne.10) then
5400 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5401 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
5402 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
5409 10 write (iout,*) "Error - singularity in solving the system",&
5410 " of equations for Lagrange multipliers."
5414 "RATTLE inactive; use -DRATTLE switch at compile time."
5417 end subroutine rattle1
5418 !-----------------------------------------------------------------------------
5420 ! RATTLE algorithm for velocity Verlet - step 2, UNRES
5424 ! implicit real*8 (a-h,o-z)
5425 ! include 'DIMENSIONS'
5427 ! include 'COMMON.CONTROL'
5428 ! include 'COMMON.VAR'
5429 ! include 'COMMON.MD'
5431 ! include 'COMMON.LANGEVIN'
5433 ! include 'COMMON.LANGEVIN.lang0'
5435 ! include 'COMMON.CHAIN'
5436 ! include 'COMMON.DERIV'
5437 ! include 'COMMON.GEO'
5438 ! include 'COMMON.LOCAL'
5439 ! include 'COMMON.INTERACT'
5440 ! include 'COMMON.IOUNITS'
5441 ! include 'COMMON.NAMES'
5442 ! include 'COMMON.TIME1'
5443 !el real(kind=8) :: gginv(2*nres,2*nres),&
5444 !el gdc(3,2*nres,2*nres)
5445 real(kind=8) :: dC_uncor(3,2*nres) !,&
5446 !el Cmat(2*nres,2*nres)
5447 real(kind=8) :: x(2*nres) !maxres2=2*maxres
5448 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
5449 !el common /przechowalnia/ nbond
5450 integer :: max_rattle = 5
5451 logical :: lprn = .false., lprn1 = .false., not_done
5452 real(kind=8) :: tol_rattle = 1.0d-5
5456 !el /common/ przechowalnia
5457 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
5458 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
5459 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
5461 if (lprn) write (iout,*) "RATTLE2"
5462 if (lprn) write (iout,*) "Velocity correction"
5463 ! Calculate the matrix G dC
5469 gdc(j,i,ind)=GGinv(i,ind)*dC(j,k)
5474 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5475 .and.(mnum.lt.4)) then
5478 gdc(j,i,ind)=GGinv(i,ind)*dC(j,k+nres)
5484 ! write (iout,*) "gdc"
5486 ! write (iout,*) "i",i
5488 ! write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
5492 ! Calculate the matrix of the system of equations
5499 Cmat(ind,j)=Cmat(ind,j)+dC(k,i)*gdc(k,ind,j)
5505 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5506 .and.(mnum.lt.4)) then
5511 Cmat(ind,j)=Cmat(ind,j)+dC(k,i+nres)*gdc(k,ind,j)
5516 ! Calculate the scalar product dC o d_t_new
5520 x(ind)=scalar(d_t(1,i),dC(1,i))
5524 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5525 .and.(mnum.lt.4)) then
5527 x(ind)=scalar(d_t(1,i+nres),dC(1,i+nres))
5531 write (iout,*) "Velocities and violations"
5535 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5536 i,ind,(d_t(j,i),j=1,3),x(ind)
5540 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5541 .and.(mnum.lt.4)) then
5543 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5544 i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind)
5550 if (dabs(x(i)).gt.xmax) then
5554 if (xmax.lt.tol_rattle) then
5559 write (iout,*) "Matrix Cmat"
5560 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
5562 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
5563 ! Add constraint term to velocities
5570 xx = xx+x(ii)*gdc(j,ind,ii)
5572 d_t(j,i)=d_t(j,i)-xx
5577 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5578 .and.(mnum.lt.4)) then
5583 xx = xx+x(ii)*gdc(j,ind,ii)
5585 d_t(j,i+nres)=d_t(j,i+nres)-xx
5591 "New velocities, Lagrange multipliers violations"
5595 if (lprn) write (iout,'(2i5,3f10.5,5x,2e15.5)') &
5596 i,ind,(d_t(j,i),j=1,3),x(ind),scalar(d_t(1,i),dC(1,i))
5600 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5603 write (iout,'(2i5,3f10.5,5x,2e15.5)') &
5604 i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind),&
5605 scalar(d_t(1,i+nres),dC(1,i+nres))
5611 10 write (iout,*) "Error - singularity in solving the system",&
5612 " of equations for Lagrange multipliers."
5616 "RATTLE inactive; use -DRATTLE option at compile time."
5619 end subroutine rattle2
5620 !-----------------------------------------------------------------------------
5621 subroutine rattle_brown
5622 ! RATTLE/LINCS algorithm for Brownian dynamics, UNRES
5626 ! implicit real*8 (a-h,o-z)
5627 ! include 'DIMENSIONS'
5629 ! include 'COMMON.CONTROL'
5630 ! include 'COMMON.VAR'
5631 ! include 'COMMON.MD'
5633 ! include 'COMMON.LANGEVIN'
5635 ! include 'COMMON.LANGEVIN.lang0'
5637 ! include 'COMMON.CHAIN'
5638 ! include 'COMMON.DERIV'
5639 ! include 'COMMON.GEO'
5640 ! include 'COMMON.LOCAL'
5641 ! include 'COMMON.INTERACT'
5642 ! include 'COMMON.IOUNITS'
5643 ! include 'COMMON.NAMES'
5644 ! include 'COMMON.TIME1'
5645 !el real(kind=8) :: gginv(2*nres,2*nres),&
5646 !el gdc(3,2*nres,2*nres)
5647 real(kind=8) :: dC_uncor(3,2*nres) !,&
5648 !el real(kind=8) :: Cmat(2*nres,2*nres)
5649 real(kind=8) :: x(2*nres) !maxres2=2*maxres
5650 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
5651 !el common /przechowalnia/ nbond
5652 integer :: max_rattle = 5
5653 logical :: lprn = .false., lprn1 = .false., not_done
5654 real(kind=8) :: tol_rattle = 1.0d-5
5658 !el /common/ przechowalnia
5659 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
5660 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
5661 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
5664 if (lprn) write (iout,*) "RATTLE_BROWN"
5667 if (itype(i,1).ne.10) nbond=nbond+1
5669 ! Make a folded form of the Ginv-matrix
5682 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
5686 if (itype(k,1).ne.10) then
5690 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
5697 if (itype(i,1).ne.10) then
5707 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
5711 if (itype(k,1).ne.10) then
5715 if (j.eq.1 .and. l.eq.1)GGinv(ii,jj)=fricmat(ind,ind1)
5723 write (iout,*) "Matrix GGinv"
5724 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
5730 if (iter.gt.max_rattle) then
5731 write (iout,*) "Error - too many iterations in RATTLE."
5734 ! Calculate the matrix C = GG**(-1) dC_old o dC
5739 dC_uncor(j,ind1)=dC(j,i)
5743 if (itype(i,1).ne.10) then
5746 dC_uncor(j,ind1)=dC(j,i+nres)
5755 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
5759 if (itype(k,1).ne.10) then
5762 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
5767 ! Calculate deviations from standard virtual-bond lengths
5771 x(ind)=vbld(i+1)**2-vbl**2
5774 if (itype(i,1).ne.10) then
5776 x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
5780 write (iout,*) "Coordinates and violations"
5782 write(iout,'(i5,3f10.5,5x,e15.5)') &
5783 i,(dC_uncor(j,i),j=1,3),x(i)
5785 write (iout,*) "Velocities and violations"
5789 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5790 i,ind,(d_t(j,i),j=1,3),scalar(d_t(1,i),dC_old(1,i))
5793 if (itype(i,1).ne.10) then
5795 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5796 i+nres,ind,(d_t(j,i+nres),j=1,3),&
5797 scalar(d_t(1,i+nres),dC_old(1,i+nres))
5800 write (iout,*) "gdc"
5802 write (iout,*) "i",i
5804 write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
5810 if (dabs(x(i)).gt.xmax) then
5814 if (xmax.lt.tol_rattle) then
5818 ! Calculate the matrix of the system of equations
5823 Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
5828 write (iout,*) "Matrix Cmat"
5829 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
5831 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
5832 ! Add constraint term to positions
5839 xx = xx+x(ii)*gdc(j,ind,ii)
5842 d_t(j,i)=d_t(j,i)+xx/d_time
5847 if (itype(i,1).ne.10) then
5852 xx = xx+x(ii)*gdc(j,ind,ii)
5855 d_t(j,i+nres)=d_t(j,i+nres)+xx/d_time
5856 dC(j,i+nres)=dC(j,i+nres)+xx
5860 ! Rebuild the chain using the new coordinates
5861 call chainbuild_cart
5863 write (iout,*) "New coordinates, Lagrange multipliers,",&
5864 " and differences between actual and standard bond lengths"
5868 xx=vbld(i+1)**2-vbl**2
5869 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
5870 i,(dC(j,i),j=1,3),x(ind),xx
5873 if (itype(i,1).ne.10) then
5875 xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
5876 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
5877 i,(dC(j,i+nres),j=1,3),x(ind),xx
5880 write (iout,*) "Velocities and violations"
5884 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5885 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
5888 if (itype(i,1).ne.10) then
5890 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5891 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
5892 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
5899 10 write (iout,*) "Error - singularity in solving the system",&
5900 " of equations for Lagrange multipliers."
5904 "RATTLE inactive; use -DRATTLE option at compile time"
5907 end subroutine rattle_brown
5908 !-----------------------------------------------------------------------------
5910 !-----------------------------------------------------------------------------
5911 subroutine friction_force
5916 ! implicit real*8 (a-h,o-z)
5917 ! include 'DIMENSIONS'
5918 ! include 'COMMON.VAR'
5919 ! include 'COMMON.CHAIN'
5920 ! include 'COMMON.DERIV'
5921 ! include 'COMMON.GEO'
5922 ! include 'COMMON.LOCAL'
5923 ! include 'COMMON.INTERACT'
5924 ! include 'COMMON.MD'
5926 ! include 'COMMON.LANGEVIN'
5928 ! include 'COMMON.LANGEVIN.lang0'
5930 ! include 'COMMON.IOUNITS'
5931 !el real(kind=8),dimension(6*nres) :: gamvec !(MAXRES6) maxres6=6*maxres
5932 !el common /syfek/ gamvec
5934 integer iposc,ichain,n,innt,inct
5935 double precision rs(nres*2)
5936 real(kind=8) ::v_work(3,6*nres),vvec(2*nres)
5938 real(kind=8) :: v_work(6*nres)
5941 real(kind=8) :: vv(3),vvtot(3,nres)!,v_work(6*nres) !,&
5942 !el ginvfric(2*nres,2*nres) !maxres2=2*maxres
5943 !el common /przechowalnia/ ginvfric
5945 logical :: lprn, checkmode
5946 integer :: i,j,ind,k,nres2,nres6,mnum
5951 ! if (large) lprn=.true.
5952 ! if (large) checkmode=.true.
5954 !c Here accelerations due to friction forces are computed right after forces.
5955 d_t_work(:6*nres)=0.0d0
5957 v_work(j,1)=d_t(j,0)
5958 v_work(j,nnt)=d_t(j,0)
5962 v_work(j,i)=v_work(j,i-1)+d_t(j,i-1)
5967 if (iabs(itype(i,1)).ne.10 .and. iabs(itype(i,mnum)).ne.ntyp1_molec(mnum).and.mnum.lt.3) then
5969 v_work(j,i+nres)=v_work(j,i)+d_t(j,i+nres)
5974 write (iout,*) "v_work"
5976 write (iout,'(i5,3f10.5)') i,(v_work(j,i),j=1,3)
5982 n=dimen_chain(ichain)
5983 iposc=iposd_chain(ichain)
5984 !c write (iout,*) "friction_force j",j," ichain",ichain,
5985 !c & " n",n," iposc",iposc,iposc+n-1
5986 innt=chain_border(1,ichain)
5987 inct=chain_border(2,ichain)
5989 !c innt=chain_border(1,1)
5990 !c inct=chain_border(2,1)
5993 vvec(ind+1)=v_work(j,i)
5995 ! if (iabs(itype(i)).ne.10) then
5996 if (iabs(itype(i,1)).ne.10 .and. iabs(itype(i,mnum)).ne.ntyp1_molec(mnum).and.mnum.lt.3) then
5997 vvec(ind+1)=v_work(j,i+nres)
6002 write (iout,*) "vvec ind",ind," n",n
6003 write (iout,'(f10.5)') (vvec(i),i=iposc,ind)
6005 !c write (iout,*) "chain",i," ind",ind," n",n
6013 ! if (large) print *,"before fivediagmult"
6014 call fivediagmult(n,DMfric(iposc),DU1fric(iposc),&
6015 DU2fric(iposc),vvec(iposc),rs)
6016 ! if (large) print *,"after fivediagmult"
6020 time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
6022 time_fricmatmult=time_fricmatmult+tcpu()-time01
6027 write (iout,'(f10.5)') (rs(i),i=1,n)
6029 do i=iposc,iposc+n-1
6030 ! write (iout,*) "ichain",ichain," i",i," j",j,&
6031 ! "index",3*(i-1)+j,"rs",rs(i-iposc+1)
6032 fric_work(3*(i-1)+j)=-rs(i-iposc+1)
6037 write (iout,*) "Vector fric_work dimen3",dimen3
6038 write (iout,'(3f10.5)') (fric_work(j),j=1,dimen3)
6041 if(.not.allocated(gamvec)) allocate(gamvec(nres6)) !(MAXRES6)
6042 if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
6050 d_t_work(j)=d_t(j,0)
6055 d_t_work(ind+j)=d_t(j,i)
6061 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
6062 .and.(mnum.lt.4)) then
6064 d_t_work(ind+j)=d_t(j,i+nres)
6070 call fricmat_mult(d_t_work,fric_work)
6072 if (.not.checkmode) return
6075 write (iout,*) "d_t_work and fric_work"
6077 write (iout,'(i3,2e15.5)') i,d_t_work(i),fric_work(i)
6081 friction(j,0)=fric_work(j)
6086 friction(j,i)=fric_work(ind+j)
6092 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
6093 .and.(mnum.lt.4)) then
6094 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
6096 friction(j,i+nres)=fric_work(ind+j)
6102 write(iout,*) "Friction backbone"
6104 write(iout,'(i5,3e15.5,5x,3e15.5)') &
6105 i,(friction(j,i),j=1,3),(d_t(j,i),j=1,3)
6107 write(iout,*) "Friction side chain"
6109 write(iout,'(i5,3e15.5,5x,3e15.5)') &
6110 i,(friction(j,i+nres),j=1,3),(d_t(j,i+nres),j=1,3)
6119 vvtot(j,i)=vv(j)+0.5d0*d_t(j,i)
6120 vvtot(j,i+nres)=vv(j)+d_t(j,i+nres)
6121 vv(j)=vv(j)+d_t(j,i)
6124 write (iout,*) "vvtot backbone and sidechain"
6126 write (iout,'(i5,3e15.5,5x,3e15.5)') i,(vvtot(j,i),j=1,3),&
6127 (vvtot(j,i+nres),j=1,3)
6132 v_work(ind+j)=vvtot(j,i)
6138 v_work(ind+j)=vvtot(j,i+nres)
6142 write (iout,*) "v_work gamvec and site-based friction forces"
6144 write (iout,'(i5,3e15.5)') i,v_work(i),gamvec(i),&
6148 ! fric_work1(i)=0.0d0
6150 ! fric_work1(i)=fric_work1(i)-A(j,i)*gamvec(j)*v_work(j)
6153 ! write (iout,*) "fric_work and fric_work1"
6155 ! write (iout,'(i5,2e15.5)') i,fric_work(i),fric_work1(i)
6161 ginvfric(i,j)=ginvfric(i,j)+ginv(i,k)*fricmat(k,j)
6165 write (iout,*) "ginvfric"
6167 write (iout,'(i5,100f8.3)') i,(ginvfric(i,j),j=1,dimen)
6169 write (iout,*) "symmetry check"
6172 write (iout,*) i,j,ginvfric(i,j)-ginvfric(j,i)
6178 end subroutine friction_force
6179 !-----------------------------------------------------------------------------
6180 subroutine setup_fricmat
6184 use control_data, only:time_Bcast
6185 use control, only:tcpu
6187 ! implicit real*8 (a-h,o-z)
6191 real(kind=8) :: time00
6193 ! include 'DIMENSIONS'
6194 ! include 'COMMON.VAR'
6195 ! include 'COMMON.CHAIN'
6196 ! include 'COMMON.DERIV'
6197 ! include 'COMMON.GEO'
6198 ! include 'COMMON.LOCAL'
6199 ! include 'COMMON.INTERACT'
6200 ! include 'COMMON.MD'
6201 ! include 'COMMON.SETUP'
6202 ! include 'COMMON.TIME1'
6203 ! integer licznik /0/
6206 ! include 'COMMON.LANGEVIN'
6208 ! include 'COMMON.LANGEVIN.lang0'
6210 ! include 'COMMON.IOUNITS'
6212 integer :: i,j,ind,ind1,m,ichain,innt,inct
6213 logical :: lprn = .true.
6214 real(kind=8) :: dtdi !el ,gamvec(2*nres)
6215 !el real(kind=8),dimension(2*nres,2*nres) :: ginvfric,fcopy
6216 ! real(kind=8),allocatable,dimension(:,:) :: fcopy
6217 !el real(kind=8),dimension(2*nres*(2*nres+1)/2) :: Ghalf !(mmaxres2) (mmaxres2=(maxres2*(maxres2+1)/2))
6218 !el common /syfek/ gamvec
6219 real(kind=8) :: work(8*2*nres)
6220 integer :: iwork(2*nres)
6221 !el common /przechowalnia/ ginvfric,Ghalf,fcopy
6222 integer :: ii,iti,k,l,nzero,nres2,nres6,ierr,mnum
6227 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
6228 if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
6229 if (fg_rank.ne.king) goto 10
6235 if(.not.allocated(gamvec)) allocate(gamvec(nres2)) !(MAXRES2)
6237 if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
6238 if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
6239 !el allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
6240 if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !maxres2=2*maxres
6242 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
6245 if (.not.allocated(DMfric)) allocate(DMfric(nres2))
6246 if (.not.allocated(DU1fric)) allocate(DU1fric(nres2))
6247 if (.not.allocated(DU2fric)) allocate(DU2fric(nres2))
6248 ! Load the friction coefficients corresponding to peptide groups
6253 gamvec(ind1)=gamp(mnum)
6256 if (molnum(nct).eq.5) then
6259 gamvec(ind1)=gamp(mnum)
6261 ! Load the friction coefficients corresponding to side chains
6265 gamsc(ntyp1_molec(j),j)=1.0d0
6272 gamvec(ii)=gamsc(iabs(iti),mnum)
6274 if (surfarea) call sdarea(gamvec)
6280 innt=chain_border(1,ichain)
6281 inct=chain_border(2,ichain)
6282 !c write (iout,*) "ichain",ichain," innt",innt," inct",inct
6285 DMfric(ind)=gamvec(innt-nnt+1)/4
6286 if (iabs(itype(innt,1)).eq.10.or.mnum.gt.2) then
6287 DMfric(ind)=DMfric(ind)+gamvec(m+innt-nnt+1)
6290 DMfric(ind+1)=gamvec(m+innt-nnt+1)
6293 !c write (iout,*) "DMfric init ind",ind
6297 DMfric(ind)=gamvec(i-nnt+1)/2
6298 if (iabs(itype(i,1)).eq.10.or.mnum.gt.2) then
6299 DMfric(ind)=DMfric(ind)+gamvec(m+i-nnt+1)
6302 DMfric(ind+1)=gamvec(m+i-nnt+1)
6306 !c write (iout,*) "DMfric endloop ind",ind
6307 if (inct.gt.innt) then
6308 DMfric(ind)=gamvec(inct-1-nnt+1)/4
6310 if (iabs(itype(inct,1)).eq.10.or.mnum.gt.2) then
6311 DMfric(ind)=DMfric(ind)+gamvec(inct+m-nnt+1)
6314 DMfric(ind+1)=gamvec(inct+m-nnt+1)
6318 !c write (iout,*) "DMfric end ind",ind
6322 ind=iposd_chain(ichain)
6323 innt=chain_border(1,ichain)
6324 inct=chain_border(2,ichain)
6327 if (iabs(itype(i,1)).ne.10.and.mnum.le.2) then
6330 DU1fric(ind)=gamvec(i-nnt+1)/4
6337 ind=iposd_chain(ichain)
6338 innt=chain_border(1,ichain)
6339 inct=chain_border(2,ichain)
6342 if (iabs(itype(i,1)).ne.10.and.mnum.le.2) then
6343 DU2fric(ind)=gamvec(i-nnt+1)/4
6344 DU2fric(ind+1)=0.0d0
6353 write(iout,*)"The upper part of the five-diagonal friction matrix"
6355 write (iout,'(a,i5)') 'Chain',ichain
6356 innt=iposd_chain(ichain)
6357 inct=iposd_chain(ichain)+dimen_chain(ichain)-1
6359 if (i.lt.inct-1) then
6360 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMfric(i),DU1fric(i),&
6362 else if (i.eq.inct-1) then
6363 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMfric(i),DU1fric(i)
6365 write (iout,'(2i3,3f10.5)') i,i-innt+1,DMfric(i)
6374 ! Zeroing out fricmat
6380 ! Load the friction coefficients corresponding to peptide groups
6385 gamvec(ind1)=gamp(mnum)
6388 if (molnum(nct).eq.5) then
6391 gamvec(ind1)=gamp(mnum)
6393 ! Load the friction coefficients corresponding to side chains
6397 gamsc(ntyp1_molec(j),j)=1.0d0
6404 gamvec(ii)=gamsc(iabs(iti),mnum)
6406 if (surfarea) call sdarea(gamvec)
6408 ! write (iout,*) "Matrix A and vector gamma"
6410 ! write (iout,'(i2,$)') i
6412 ! write (iout,'(f4.1,$)') A(i,j)
6414 ! write (iout,'(f8.3)') gamvec(i)
6418 write (iout,*) "Vector gamvec"
6420 write (iout,'(i5,f10.5)') i, gamvec(i)
6424 ! The friction matrix
6429 dtdi=dtdi+A(j,k)*A(j,i)*gamvec(j)
6436 write (iout,'(//a)') "Matrix fricmat"
6437 call matout2(dimen,dimen,nres2,nres2,fricmat)
6439 if (lang.eq.2 .or. lang.eq.3) then
6440 ! Mass-scale the friction matrix if non-direct integration will be performed
6446 Ginvfric(i,j)=Ginvfric(i,j)+ &
6447 Gsqrm(i,k)*Gsqrm(l,j)*fricmat(k,l)
6452 ! Diagonalize the friction matrix
6457 Ghalf(ind)=Ginvfric(i,j)
6460 call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
6463 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
6464 " mass-scaled friction matrix"
6465 call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
6467 ! Precompute matrices for tinker stochastic integrator
6474 mt1(i,j)=mt1(i,j)+fricvec(k,i)*gsqrm(k,j)
6475 mt2(i,j)=mt2(i,j)+fricvec(k,i)*gsqrp(k,j)
6481 else if (lang.eq.4) then
6482 ! Diagonalize the friction matrix
6487 Ghalf(ind)=fricmat(i,j)
6490 call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
6493 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
6495 call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
6497 ! Determine the number of zero eigenvalues of the friction matrix
6498 nzero=max0(dimen-dimen1,0)
6499 ! do while (fricgam(nzero+1).le.1.0d-5 .and. nzero.lt.dimen)
6502 write (iout,*) "Number of zero eigenvalues:",nzero
6507 fricmat(i,j)=fricmat(i,j) &
6508 +fricvec(i,k)*fricvec(j,k)/fricgam(k)
6513 write (iout,'(//a)') "Generalized inverse of fricmat"
6514 call matout(dimen,dimen,nres6,nres6,fricmat)
6519 if (nfgtasks.gt.1) then
6520 if (fg_rank.eq.0) then
6521 ! The matching BROADCAST for fg processors is called in ERGASTULUM
6527 call MPI_Bcast(10,1,MPI_INTEGER,king,FG_COMM,IERROR)
6529 time_Bcast=time_Bcast+MPI_Wtime()-time00
6531 time_Bcast=time_Bcast+tcpu()-time00
6533 ! print *,"Processor",myrank,
6534 ! & " BROADCAST iorder in SETUP_FRICMAT"
6537 write (iout,*) "setup_fricmat licznik"!,licznik !sp
6543 ! Scatter the friction matrix
6544 call MPI_Scatterv(fricmat(1,1),nginv_counts(0),&
6545 nginv_start(0),MPI_DOUBLE_PRECISION,fcopy(1,1),&
6546 myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
6549 time_scatter=time_scatter+MPI_Wtime()-time00
6550 time_scatter_fmat=time_scatter_fmat+MPI_Wtime()-time00
6552 time_scatter=time_scatter+tcpu()-time00
6553 time_scatter_fmat=time_scatter_fmat+tcpu()-time00
6557 do j=1,2*my_ng_count
6558 fricmat(j,i)=fcopy(i,j)
6561 ! write (iout,*) "My chunk of fricmat"
6562 ! call MATOUT2(my_ng_count,dimen,maxres2,maxres2,fcopy)
6567 end subroutine setup_fricmat
6568 !-----------------------------------------------------------------------------
6569 subroutine stochastic_force(stochforcvec)
6572 use random, only:anorm_distr
6573 ! implicit real*8 (a-h,o-z)
6574 ! include 'DIMENSIONS'
6575 use control, only: tcpu
6580 ! include 'COMMON.VAR'
6581 ! include 'COMMON.CHAIN'
6582 ! include 'COMMON.DERIV'
6583 ! include 'COMMON.GEO'
6584 ! include 'COMMON.LOCAL'
6585 ! include 'COMMON.INTERACT'
6586 ! include 'COMMON.MD'
6587 ! include 'COMMON.TIME1'
6589 ! include 'COMMON.LANGEVIN'
6591 ! include 'COMMON.LANGEVIN.lang0'
6593 ! include 'COMMON.IOUNITS'
6595 real(kind=8) :: x,sig,lowb,highb
6596 real(kind=8) :: ff(3),force(3,0:2*nres),zeta2,lowb2
6597 real(kind=8) :: highb2,sig2,forcvec(6*nres),stochforcvec(6*nres)
6598 real(kind=8) :: time00
6599 logical :: lprn = .false.
6600 integer :: i,j,ind,mnum
6602 integer ichain,innt,inct,iposc
6607 stochforc(j,i)=0.0d0
6617 ! Compute the stochastic forces acting on bodies. Store in force.
6623 force(j,i)=anorm_distr(x,sig,lowb,highb)
6631 force(j,i+nres)=anorm_distr(x,sig2,lowb2,highb2)
6635 time_fsample=time_fsample+MPI_Wtime()-time00
6637 time_fsample=time_fsample+tcpu()-time00
6642 innt=chain_border(1,ichain)
6643 inct=chain_border(2,ichain)
6644 iposc=iposd_chain(ichain)
6645 !c for debugging only
6646 !c innt=chain_border(1,1)
6647 !c inct=chain_border(2,1)
6648 !c iposc=iposd_chain(1)
6649 !c write (iout,*)"stochastic_force ichain=",ichain," innt",innt,
6650 !c & " inct",inct," iposc",iposc
6652 stochforcvec(ind+j)=0.5d0*force(j,innt)
6654 if (iabs(itype(innt,molnum(innt))).eq.10.or.molnum(innt).ge.3) then
6656 stochforcvec(ind+j)=stochforcvec(ind+j)+force(j,innt+nres)
6662 stochforcvec(ind+j)=force(j,innt+nres)
6668 stochforcvec(ind+j)=0.5d0*(force(j,i)+force(j,i-1))
6670 if (iabs(itype(i,1).eq.10).or.molnum(i).ge.3) then
6672 stochforcvec(ind+j)=stochforcvec(ind+j)+force(j,i+nres)
6678 stochforcvec(ind+j)=force(j,i+nres)
6684 stochforcvec(ind+j)=0.5d0*force(j,inct-1)
6686 if (iabs(itype(inct,molnum(inct))).eq.10.or.molnum(inct).ge.3) then
6688 stochforcvec(ind+j)=stochforcvec(ind+j)+force(j,inct+nres)
6694 stochforcvec(ind+j)=force(j,inct+nres)
6698 !c write (iout,*) "chain",ichain," ind",ind
6701 write (iout,*) "stochforcvec"
6702 write (iout,'(3f10.5)') (stochforcvec(j),j=1,ind)
6705 ! Compute the stochastic forces acting on virtual-bond vectors.
6711 stochforc(j,i)=ff(j)+0.5d0*force(j,i)
6714 ff(j)=ff(j)+force(j,i)
6716 ! if (itype(i+1,1).ne.ntyp1) then
6718 if (itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
6720 stochforc(j,i)=stochforc(j,i)+force(j,i+nres+1)
6721 ff(j)=ff(j)+force(j,i+nres+1)
6726 stochforc(j,0)=ff(j)+force(j,nnt+nres)
6730 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
6731 .and.(mnum.lt.4)) then
6732 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
6734 stochforc(j,i+nres)=force(j,i+nres)
6740 stochforcvec(j)=stochforc(j,0)
6745 stochforcvec(ind+j)=stochforc(j,i)
6751 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
6752 .and.(mnum.lt.4)) then
6753 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
6755 stochforcvec(ind+j)=stochforc(j,i+nres)
6761 write (iout,*) "stochforcvec"
6763 write(iout,'(i5,e15.5)') i,stochforcvec(i)
6765 write(iout,*) "Stochastic forces backbone"
6767 write(iout,'(i5,3e15.5)') i,(stochforc(j,i),j=1,3)
6769 write(iout,*) "Stochastic forces side chain"
6771 write(iout,'(i5,3e15.5)') &
6772 i,(stochforc(j,i+nres),j=1,3)
6780 write (iout,*) i,ind
6782 forcvec(ind+j)=force(j,i)
6787 write (iout,*) i,ind
6789 forcvec(j+ind)=force(j,i+nres)
6794 write (iout,*) "forcvec"
6798 write (iout,'(2i3,2f10.5)') i,j,force(j,i),&
6805 write (iout,'(2i3,2f10.5)') i,j,force(j,i+nres),&
6814 end subroutine stochastic_force
6815 !-----------------------------------------------------------------------------
6816 subroutine sdarea(gamvec)
6818 ! Scale the friction coefficients according to solvent accessible surface areas
6819 ! Code adapted from TINKER
6823 ! implicit real*8 (a-h,o-z)
6824 ! include 'DIMENSIONS'
6825 ! include 'COMMON.CONTROL'
6826 ! include 'COMMON.VAR'
6827 ! include 'COMMON.MD'
6829 ! include 'COMMON.LANGEVIN'
6831 ! include 'COMMON.LANGEVIN.lang0'
6833 ! include 'COMMON.CHAIN'
6834 ! include 'COMMON.DERIV'
6835 ! include 'COMMON.GEO'
6836 ! include 'COMMON.LOCAL'
6837 ! include 'COMMON.INTERACT'
6838 ! include 'COMMON.IOUNITS'
6839 ! include 'COMMON.NAMES'
6840 real(kind=8),dimension(2*nres) :: radius,gamvec !(maxres2)
6841 real(kind=8),parameter :: twosix = 1.122462048309372981d0
6842 logical :: lprn = .false.
6843 real(kind=8) :: probe,area,ratio
6844 integer :: i,j,ind,iti,mnum
6846 ! determine new friction coefficients every few SD steps
6848 ! set the atomic radii to estimates of sigma values
6850 ! print *,"Entered sdarea"
6856 ! Load peptide group radii
6859 radius(i)=pstok(mnum)
6861 ! Load side chain radii
6865 radius(i+nres)=restok(iti,mnum)
6868 ! write (iout,*) "i",i," radius",radius(i)
6871 radius(i) = radius(i) / twosix
6872 if (radius(i) .ne. 0.0d0) radius(i) = radius(i) + probe
6875 ! scale atomic friction coefficients by accessible area
6877 if (lprn) write (iout,*) &
6878 "Original gammas, surface areas, scaling factors, new gammas, ",&
6879 "std's of stochastic forces"
6882 if (radius(i).gt.0.0d0) then
6883 call surfatom (i,area,radius)
6884 ratio = dmax1(area/(4.0d0*pi*radius(i)**2),1.0d-1)
6885 if (lprn) write (iout,'(i5,3f10.5,$)') &
6886 i,gamvec(ind+1),area,ratio
6889 gamvec(ind) = ratio * gamvec(ind)
6892 stdforcp(i)=stdfp(mnum)*dsqrt(gamvec(ind))
6893 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcp(i)
6897 if (radius(i+nres).gt.0.0d0) then
6898 call surfatom (i+nres,area,radius)
6899 ratio = dmax1(area/(4.0d0*pi*radius(i+nres)**2),1.0d-1)
6900 if (lprn) write (iout,'(i5,3f10.5,$)') &
6901 i,gamvec(ind+1),area,ratio
6904 gamvec(ind) = ratio * gamvec(ind)
6907 stdforcsc(i)=stdfsc(itype(i,mnum),mnum)*dsqrt(gamvec(ind))
6908 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcsc(i)
6913 end subroutine sdarea
6914 !-----------------------------------------------------------------------------
6916 !-----------------------------------------------------------------------------
6919 ! ###################################################
6920 ! ## COPYRIGHT (C) 1996 by Jay William Ponder ##
6921 ! ## All Rights Reserved ##
6922 ! ###################################################
6924 ! ################################################################
6926 ! ## subroutine surfatom -- exposed surface area of an atom ##
6928 ! ################################################################
6931 ! "surfatom" performs an analytical computation of the surface
6932 ! area of a specified atom; a simplified version of "surface"
6934 ! literature references:
6936 ! T. J. Richmond, "Solvent Accessible Surface Area and
6937 ! Excluded Volume in Proteins", Journal of Molecular Biology,
6940 ! L. Wesson and D. Eisenberg, "Atomic Solvation Parameters
6941 ! Applied to Molecular Dynamics of Proteins in Solution",
6942 ! Protein Science, 1, 227-235 (1992)
6944 ! variables and parameters:
6946 ! ir number of atom for which area is desired
6947 ! area accessible surface area of the atom
6948 ! radius radii of each of the individual atoms
6951 subroutine surfatom(ir,area,radius)
6953 ! implicit real*8 (a-h,o-z)
6954 ! include 'DIMENSIONS'
6956 ! include 'COMMON.GEO'
6957 ! include 'COMMON.IOUNITS'
6959 integer :: nsup,nstart_sup
6960 ! double precision c,dc,dc_old,d_c_work,xloc,xrot,dc_norm
6961 ! common /chain/ c(3,maxres2+2),dc(3,0:maxres2),dc_old(3,0:maxres2),
6962 ! & xloc(3,maxres),xrot(3,maxres),dc_norm(3,0:maxres2),
6963 ! & dc_work(MAXRES6),nres,nres0
6964 integer,parameter :: maxarc=300
6968 integer :: mi,ni,narc
6969 integer :: key(maxarc)
6970 integer :: intag(maxarc)
6971 integer :: intag1(maxarc)
6972 real(kind=8) :: area,arcsum
6973 real(kind=8) :: arclen,exang
6974 real(kind=8) :: delta,delta2
6975 real(kind=8) :: eps,rmove
6976 real(kind=8) :: xr,yr,zr
6977 real(kind=8) :: rr,rrsq
6978 real(kind=8) :: rplus,rminus
6979 real(kind=8) :: axx,axy,axz
6980 real(kind=8) :: ayx,ayy
6981 real(kind=8) :: azx,azy,azz
6982 real(kind=8) :: uxj,uyj,uzj
6983 real(kind=8) :: tx,ty,tz
6984 real(kind=8) :: txb,tyb,td
6985 real(kind=8) :: tr2,tr,txr,tyr
6986 real(kind=8) :: tk1,tk2
6987 real(kind=8) :: thec,the,t,tb
6988 real(kind=8) :: txk,tyk,tzk
6989 real(kind=8) :: t1,ti,tf,tt
6990 real(kind=8) :: txj,tyj,tzj
6991 real(kind=8) :: ccsq,cc,xysq
6992 real(kind=8) :: bsqk,bk,cosine
6993 real(kind=8) :: dsqj,gi,pix2
6994 real(kind=8) :: therk,dk,gk
6995 real(kind=8) :: risqk,rik
6996 real(kind=8) :: radius(2*nres) !(maxatm) (maxatm=maxres2)
6997 real(kind=8) :: ri(maxarc),risq(maxarc)
6998 real(kind=8) :: ux(maxarc),uy(maxarc),uz(maxarc)
6999 real(kind=8) :: xc(maxarc),yc(maxarc),zc(maxarc)
7000 real(kind=8) :: xc1(maxarc),yc1(maxarc),zc1(maxarc)
7001 real(kind=8) :: dsq(maxarc),bsq(maxarc)
7002 real(kind=8) :: dsq1(maxarc),bsq1(maxarc)
7003 real(kind=8) :: arci(maxarc),arcf(maxarc)
7004 real(kind=8) :: ex(maxarc),lt(maxarc),gr(maxarc)
7005 real(kind=8) :: b(maxarc),b1(maxarc),bg(maxarc)
7006 real(kind=8) :: kent(maxarc),kout(maxarc)
7007 real(kind=8) :: ther(maxarc)
7008 logical :: moved,top
7009 logical :: omit(maxarc)
7012 ! maxatm = 2*nres !maxres2 maxres2=2*maxres
7013 ! maxlight = 8*maxatm
7016 ! maxtors = 4*maxatm
7018 ! zero out the surface area for the sphere of interest
7021 ! write (2,*) "ir",ir," radius",radius(ir)
7022 if (radius(ir) .eq. 0.0d0) return
7024 ! set the overlap significance and connectivity shift
7028 delta2 = delta * delta
7033 ! store coordinates and radius of the sphere of interest
7041 ! initialize values of some counters and summations
7050 ! test each sphere to see if it overlaps the sphere of interest
7053 if (i.eq.ir .or. radius(i).eq.0.0d0) goto 30
7054 rplus = rr + radius(i)
7056 if (abs(tx) .ge. rplus) goto 30
7058 if (abs(ty) .ge. rplus) goto 30
7060 if (abs(tz) .ge. rplus) goto 30
7062 ! check for sphere overlap by testing distance against radii
7064 xysq = tx*tx + ty*ty
7065 if (xysq .lt. delta2) then
7072 if (rplus-cc .le. delta) goto 30
7073 rminus = rr - radius(i)
7075 ! check to see if sphere of interest is completely buried
7077 if (cc-abs(rminus) .le. delta) then
7078 if (rminus .le. 0.0d0) goto 170
7082 ! check for too many overlaps with sphere of interest
7084 if (io .ge. maxarc) then
7086 20 format (/,' SURFATOM -- Increase the Value of MAXARC')
7090 ! get overlap between current sphere and sphere of interest
7099 gr(io) = (ccsq+rplus*rminus) / (2.0d0*rr*b1(io))
7105 ! case where no other spheres overlap the sphere of interest
7108 area = 4.0d0 * pi * rrsq
7112 ! case where only one sphere overlaps the sphere of interest
7115 area = pix2 * (1.0d0 + gr(1))
7116 area = mod(area,4.0d0*pi) * rrsq
7120 ! case where many spheres intersect the sphere of interest;
7121 ! sort the intersecting spheres by their degree of overlap
7123 call sort2 (io,gr,key)
7126 intag(i) = intag1(k)
7135 ! get radius of each overlap circle on surface of the sphere
7140 risq(i) = rrsq - gi*gi
7141 ri(i) = sqrt(risq(i))
7142 ther(i) = 0.5d0*pi - asin(min(1.0d0,max(-1.0d0,gr(i))))
7145 ! find boundary of inaccessible area on sphere of interest
7148 if (.not. omit(k)) then
7155 ! check to see if J circle is intersecting K circle;
7156 ! get distance between circle centers and sum of radii
7159 if (omit(j)) goto 60
7160 cc = (txk*xc(j)+tyk*yc(j)+tzk*zc(j))/(bk*b(j))
7161 cc = acos(min(1.0d0,max(-1.0d0,cc)))
7162 td = therk + ther(j)
7164 ! check to see if circles enclose separate regions
7166 if (cc .ge. td) goto 60
7168 ! check for circle J completely inside circle K
7170 if (cc+ther(j) .lt. therk) goto 40
7172 ! check for circles that are essentially parallel
7174 if (cc .gt. delta) goto 50
7179 ! check to see if sphere of interest is completely buried
7182 if (pix2-cc .le. td) goto 170
7188 ! find T value of circle intersections
7191 if (omit(k)) goto 110
7206 ! rotation matrix elements
7218 if (.not. omit(j)) then
7223 ! rotate spheres so K vector colinear with z-axis
7225 uxj = txj*axx + tyj*axy - tzj*axz
7226 uyj = tyj*ayy - txj*ayx
7227 uzj = txj*azx + tyj*azy + tzj*azz
7228 cosine = min(1.0d0,max(-1.0d0,uzj/b(j)))
7229 if (acos(cosine) .lt. therk+ther(j)) then
7230 dsqj = uxj*uxj + uyj*uyj
7235 tr2 = risqk*dsqj - tb*tb
7241 ! get T values of intersection for K circle
7244 tb = min(1.0d0,max(-1.0d0,tb))
7246 if (tyb-txr .lt. 0.0d0) tk1 = pix2 - tk1
7248 tb = min(1.0d0,max(-1.0d0,tb))
7250 if (tyb+txr .lt. 0.0d0) tk2 = pix2 - tk2
7251 thec = (rrsq*uzj-gk*bg(j)) / (rik*ri(j)*b(j))
7252 if (abs(thec) .lt. 1.0d0) then
7254 else if (thec .ge. 1.0d0) then
7256 else if (thec .le. -1.0d0) then
7260 ! see if "tk1" is entry or exit point; check t=0 point;
7261 ! "ti" is exit point, "tf" is entry point
7263 cosine = min(1.0d0,max(-1.0d0, &
7264 (uzj*gk-uxj*rik)/(b(j)*rr)))
7265 if ((acos(cosine)-ther(j))*(tk2-tk1) .le. 0.0d0) then
7273 if (narc .ge. maxarc) then
7275 70 format (/,' SURFATOM -- Increase the Value',&
7279 if (tf .le. ti) then
7300 ! special case; K circle without intersections
7302 if (narc .le. 0) goto 90
7304 ! general case; sum up arclength and set connectivity code
7306 call sort2 (narc,arci,key)
7311 if (narc .gt. 1) then
7314 if (t .lt. arci(j)) then
7315 arcsum = arcsum + arci(j) - t
7316 exang = exang + ex(ni)
7318 if (jb .ge. maxarc) then
7320 80 format (/,' SURFATOM -- Increase the Value',&
7325 kent(jb) = maxarc*i + k
7327 kout(jb) = maxarc*k + i
7336 arcsum = arcsum + pix2 - t
7338 exang = exang + ex(ni)
7341 kent(jb) = maxarc*i + k
7343 kout(jb) = maxarc*k + i
7350 arclen = arclen + gr(k)*arcsum
7353 if (arclen .eq. 0.0d0) goto 170
7354 if (jb .eq. 0) goto 150
7356 ! find number of independent boundaries and check connectivity
7360 if (kout(k) .ne. 0) then
7367 if (m .eq. kent(ii)) then
7370 if (j .eq. jb) goto 150
7382 ! attempt to fix connectivity error by moving atom slightly
7386 140 format (/,' SURFATOM -- Connectivity Error at Atom',i6)
7395 ! compute the exposed surface area for the sphere of interest
7398 area = ib*pix2 + exang + arclen
7399 area = mod(area,4.0d0*pi) * rrsq
7401 ! attempt to fix negative area by moving atom slightly
7403 if (area .lt. 0.0d0) then
7406 160 format (/,' SURFATOM -- Negative Area at Atom',i6)
7417 end subroutine surfatom
7418 !----------------------------------------------------------------
7419 !----------------------------------------------------------------
7420 subroutine alloc_MD_arrays
7421 !EL Allocation of arrays used by MD module
7423 integer :: nres2,nres6
7426 !----------------------
7430 allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
7431 allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
7432 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
7433 allocate(fricvec(nres2,nres2))
7434 allocate(pfric_mat(nres2,nres2),vfric_mat(nres2,nres2))
7435 allocate(afric_mat(nres2,nres2),prand_mat(nres2,nres2))
7436 allocate(vrand_mat1(nres2,nres2),vrand_mat2(nres2,nres2)) !(MAXRES2,MAXRES2)
7437 allocate(pfric0_mat(nres2,nres2,0:maxflag_stoch))
7438 allocate(afric0_mat(nres2,nres2,0:maxflag_stoch))
7439 allocate(vfric0_mat(nres2,nres2,0:maxflag_stoch))
7440 allocate(prand0_mat(nres2,nres2,0:maxflag_stoch))
7441 allocate(vrand0_mat1(nres2,nres2,0:maxflag_stoch))
7442 allocate(vrand0_mat2(nres2,nres2,0:maxflag_stoch)) !(MAXRES2,MAXRES2,0:maxflag_stoch)
7443 allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
7445 allocate(mt1(nres2,nres2),mt2(nres2,nres2),mt3(nres2,nres2)) !(maxres2,maxres2)
7446 !----------------------
7448 ! commom.langevin.lang0
7450 allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
7452 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
7453 allocate(fricvec(nres2,nres2)) !(MAXRES2,MAXRES2)
7455 allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
7456 allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
7459 if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2))
7461 !----------------------
7462 ! commom.hairpin in CSA module
7463 !----------------------
7464 ! common.mce in MCM_MD module
7465 !----------------------
7467 ! common /mdgrad/ in module.energy
7468 ! common /back_constr/ in module.energy
7469 ! common /qmeas/ in module.energy
7472 allocate(potEcomp(0:n_ene+4)) !(0:n_ene+4)
7474 allocate(d_t(3,0:nres2),d_a(3,0:nres2),d_t_old(3,0:nres2)) !(3,0:MAXRES2)
7475 allocate(d_a_work(nres6)) !(6*MAXRES)
7477 allocate(DM(nres2),DU1(nres2),DU2(nres2))
7478 allocate(DMorig(nres2),DU1orig(nres2),DU2orig(nres2))
7480 allocate(Gvec(1300,1300))!maximum number of protein in one chain
7483 write (iout,*) "Before A Allocation",nfgtasks-1
7485 allocate(Gmat(nres2,nres2),A(nres2,nres2))
7486 if(.not.allocated(Ginv)) allocate(Ginv(nres2,nres2)) !in control: ergastulum
7487 allocate(Gsqrp(nres2,nres2),Gsqrm(nres2,nres2),Gvec(nres2,nres2)) !(maxres2,maxres2)
7489 allocate(Geigen(nres2)) !(maxres2)
7490 if(.not.allocated(vtot)) allocate(vtot(nres2)) !(maxres2)
7491 ! common /inertia/ in io_conf: parmread
7492 ! real(kind=8),dimension(:),allocatable :: ISC,msc !(ntyp+1)
7493 ! common /langevin/in io read_MDpar
7494 ! real(kind=8),dimension(:),allocatable :: gamsc !(ntyp1)
7495 ! real(kind=8),dimension(:),allocatable :: stdfsc !(ntyp)
7496 ! in io_conf: parmread
7497 ! real(kind=8),dimension(:),allocatable :: restok !(ntyp+1)
7498 ! common /mdpmpi/ in control: ergastulum
7499 if(.not.allocated(ng_start)) allocate(ng_start(0:nfgtasks-1))
7500 if(.not.allocated(ng_counts)) allocate(ng_counts(0:nfgtasks-1))
7501 if(.not.allocated(nginv_counts)) allocate(nginv_counts(0:nfgtasks-1)) !(0:MaxProcs-1)
7502 if(.not.allocated(nginv_start)) allocate(nginv_start(0:nfgtasks)) !(0:MaxProcs)
7503 !----------------------
7504 ! common.muca in read_muca
7505 ! common /double_muca/
7506 ! real(kind=8) :: elow,ehigh,factor,hbin,factor_min
7507 ! real(kind=8),dimension(:),allocatable :: emuca,nemuca,&
7508 ! nemuca2,hist !(4*maxres)
7509 ! common /integer_muca/
7510 ! integer :: nmuca,imtime,muca_smooth
7512 ! real(kind=8),dimension(:),allocatable :: elowi,ehighi !(maxprocs)
7513 !----------------------
7515 ! common /mdgrad/ in module.energy
7516 ! common /back_constr/ in module.energy
7517 ! common /qmeas/ in module.energy
7521 allocate(d_t_work(nres6),d_t_work_new(nres6),d_af_work(nres6))
7522 allocate(d_as_work(nres6),kinetic_force(nres6)) !(MAXRES6)
7523 allocate(d_t_new(3,0:nres2),d_a_old(3,0:nres2),d_a_short(3,0:nres2)) !,d_a !(3,0:MAXRES2)
7524 allocate(stdforcp(nres),stdforcsc(nres)) !(MAXRES)
7525 !----------------------
7527 allocate(D_ban(nres6)) !(MAXRES6) maxres6=6*maxres
7528 ! common /stochcalc/ stochforcvec
7529 allocate(stochforcvec(nres6)) !(MAXRES6) maxres6=6*maxres
7530 !----------------------
7532 end subroutine alloc_MD_arrays
7533 !-----------------------------------------------------------------------------
7534 !-----------------------------------------------------------------------------