2 !-----------------------------------------------------------------------------
15 !-----------------------------------------------------------------------------
17 ! common /mdgrad/ in module.energy
18 ! common /back_constr/ in module.energy
19 ! common /qmeas/ in module.energy
23 real(kind=8),dimension(:),allocatable :: d_t_work,&
24 d_t_work_new,d_af_work,d_as_work,kinetic_force !(MAXRES6)
25 real(kind=8),dimension(:,:),allocatable :: d_t_new,&
26 d_a_old,d_a_short!,d_a !(3,0:MAXRES2)
27 ! real(kind=8),dimension(:),allocatable :: d_a_work !(6*MAXRES)
28 ! real(kind=8),dimension(:,:),allocatable :: Gmat,Ginv,A,&
29 ! Gsqrp,Gsqrm,Gvec !(maxres2,maxres2)
30 ! real(kind=8),dimension(:),allocatable :: Geigen !(maxres2)
31 ! integer :: dimen,dimen1,dimen3
32 ! integer :: lang,count_reset_moment,count_reset_vel
33 ! logical :: reset_moment,reset_vel,rattle,RESPA
36 ! real(kind=8) :: rwat,etawat,stdfp,cPoise
37 ! real(kind=8),dimension(:),allocatable :: gamsc !(ntyp1)
38 ! real(kind=8),dimension(:),allocatable :: stdfsc !(ntyp)
39 real(kind=8),dimension(:),allocatable :: stdforcp,stdforcsc !(MAXRES)
40 !-----------------------------------------------------------------------------
44 ! ###################################################
45 ! ## COPYRIGHT (C) 1992 by Jay William Ponder ##
46 ! ## All Rights Reserved ##
47 ! ###################################################
49 ! #############################################################
51 ! ## sizes.i -- parameter values to set array dimensions ##
53 ! #############################################################
56 ! "sizes.i" sets values for critical array dimensions used
57 ! throughout the software; these parameters will fix the size
58 ! of the largest systems that can be handled; values too large
59 ! for the computer's memory and/or swap space to accomodate
60 ! will result in poor performance or outright failure
62 ! parameter: maximum allowed number of:
64 ! maxatm atoms in the molecular system
65 ! maxval atoms directly bonded to an atom
66 ! maxgrp ! user-defined groups of atoms
67 ! maxtyp force field atom type definitions
68 ! maxclass force field atom class definitions
69 ! maxkey lines in the keyword file
70 ! maxrot bonds for torsional rotation
71 ! maxvar optimization variables (vector storage)
72 ! maxopt optimization variables (matrix storage)
73 ! maxhess off-diagonal Hessian elements
74 ! maxlight sites for method of lights neighbors
75 ! maxvib vibrational frequencies
76 ! maxgeo distance geometry points
77 ! maxcell unit cells in replicated crystal
78 ! maxring 3-, 4-, or 5-membered rings
79 ! maxfix geometric restraints
80 ! maxbio biopolymer atom definitions
81 ! maxres residues in the macromolecule
82 ! maxamino amino acid residue types
83 ! maxnuc nucleic acid residue types
84 ! maxbnd covalent bonds in molecular system
85 ! maxang bond angles in molecular system
86 ! maxtors torsional angles in molecular system
87 ! maxpi atoms in conjugated pisystem
88 ! maxpib covalent bonds involving pisystem
89 ! maxpit torsional angles involving pisystem
92 !el integer maxatm,maxval,maxgrp
93 !el integer maxtyp,maxclass,maxkey
94 !el integer maxrot,maxopt
95 !el integer maxhess,maxlight,maxvib
96 !el integer maxgeo,maxcell,maxring
97 !el integer maxfix,maxbio
98 !el integer maxamino,maxnuc,maxbnd
99 !el integer maxang,maxtors,maxpi
100 !el integer maxpib,maxpit
101 ! integer :: maxatm !=2*nres !maxres2 maxres2=2*maxres
102 ! integer,parameter :: maxval=8
103 ! integer,parameter :: maxgrp=1000
104 ! integer,parameter :: maxtyp=3000
105 ! integer,parameter :: maxclass=500
106 ! integer,parameter :: maxkey=10000
107 ! integer,parameter :: maxrot=1000
108 ! integer,parameter :: maxopt=1000
109 ! integer,parameter :: maxhess=1000000
110 ! integer :: maxlight !=8*maxatm
111 ! integer,parameter :: maxvib=1000
112 ! integer,parameter :: maxgeo=1000
113 ! integer,parameter :: maxcell=10000
114 ! integer,parameter :: maxring=10000
115 ! integer,parameter :: maxfix=10000
116 ! integer,parameter :: maxbio=10000
117 ! integer,parameter :: maxamino=31
118 ! integer,parameter :: maxnuc=12
119 ! integer :: maxbnd !=2*maxatm
120 ! integer :: maxang !=3*maxatm
121 ! integer :: maxtors !=4*maxatm
122 ! integer,parameter :: maxpi=100
123 ! integer,parameter :: maxpib=2*maxpi
124 ! integer,parameter :: maxpit=4*maxpi
125 !-----------------------------------------------------------------------------
126 ! Maximum number of seed
127 ! integer,parameter :: max_seed=1
128 !-----------------------------------------------------------------------------
129 real(kind=8),dimension(:),allocatable :: stochforcvec !(MAXRES6) maxres6=6*maxres
130 ! common /stochcalc/ stochforcvec
131 !-----------------------------------------------------------------------------
132 ! common /przechowalnia/ subroutines: rattle1,rattle2,rattle_brown
133 real(kind=8),dimension(:,:),allocatable :: GGinv !(2*nres,2*nres) maxres2=2*maxres
134 real(kind=8),dimension(:,:,:),allocatable :: gdc !(3,2*nres,2*nres) maxres2=2*maxres
135 real(kind=8),dimension(:,:),allocatable :: Cmat !(2*nres,2*nres) maxres2=2*maxres
136 !-----------------------------------------------------------------------------
137 ! common /syfek/ subroutines: friction_force,setup_fricmat
138 !el real(kind=8),dimension(:),allocatable :: gamvec !(MAXRES6) or (MAXRES2)
139 !-----------------------------------------------------------------------------
140 ! common /przechowalnia/ subroutines: friction_force,setup_fricmat
141 real(kind=8),dimension(:,:),allocatable :: ginvfric !(2*nres,2*nres) !maxres2=2*maxres
142 !-----------------------------------------------------------------------------
143 ! common /przechowalnia/ subroutine: setup_fricmat
144 real(kind=8),dimension(:,:),allocatable :: fcopy !(2*nres,2*nres)
145 !-----------------------------------------------------------------------------
148 !-----------------------------------------------------------------------------
150 !-----------------------------------------------------------------------------
152 !-----------------------------------------------------------------------------
153 subroutine brown_step(itime)
154 !------------------------------------------------
155 ! Perform a single Euler integration step of Brownian dynamics
156 !------------------------------------------------
157 ! implicit real*8 (a-h,o-z)
159 use control, only: tcpu
162 ! use io_conf, only:cartprint
163 ! include 'DIMENSIONS'
167 ! include 'COMMON.CONTROL'
168 ! include 'COMMON.VAR'
169 ! include 'COMMON.MD'
171 ! include 'COMMON.LANGEVIN'
173 ! include 'COMMON.LANGEVIN.lang0'
175 ! include 'COMMON.CHAIN'
176 ! include 'COMMON.DERIV'
177 ! include 'COMMON.GEO'
178 ! include 'COMMON.LOCAL'
179 ! include 'COMMON.INTERACT'
180 ! include 'COMMON.IOUNITS'
181 ! include 'COMMON.NAMES'
182 ! include 'COMMON.TIME1'
183 real(kind=8),dimension(6*nres) :: zapas !(MAXRES6) maxres6=6*maxres
184 integer :: rstcount !ilen,
186 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
187 real(kind=8),dimension(6*nres,2*nres) :: Bmat,GBmat,Tmat !(MAXRES6,MAXRES2) (maxres2=2*maxres,maxres6=6*maxres)
188 real(kind=8),dimension(2*nres,2*nres) :: Cmat_,Cinv !(maxres2,maxres2) maxres2=2*maxres
189 real(kind=8),dimension(6*nres,6*nres) :: Pmat !(maxres6,maxres6) maxres6=6*maxres
190 ! real(kind=8),dimension(:,:),allocatable :: Bmat,GBmat,Tmat !(MAXRES6,MAXRES2) (maxres2=2*maxres,maxres6=6*maxres)
191 ! real(kind=8),dimension(:,:),allocatable :: Cmat_,Cinv !(maxres2,maxres2) maxres2=2*maxres
192 ! real(kind=8),dimension(:,:),allocatable :: Pmat !(maxres6,maxres6) maxres6=6*maxres
193 real(kind=8),dimension(6*nres) :: Td !(maxres6) maxres6=6*maxres
194 real(kind=8),dimension(2*nres) :: ppvec !(maxres2) maxres2=2*maxres
195 !el common /stochcalc/ stochforcvec
196 !el real(kind=8),dimension(3) :: cm !el
197 !el common /gucio/ cm
199 logical :: lprn = .false.,lprn1 = .false.
200 integer :: maxiter = 5
201 real(kind=8) :: difftol = 1.0d-5
202 real(kind=8) :: xx,diffmax,blen2,diffbond,tt0
203 integer :: i,j,nbond,k,ind,ind1,iter
204 integer :: nres2,nres6
208 ! if (.not.allocated(Bmat)) allocate(Bmat(nres6,nres2))
209 ! if (.not.allocated(GBmat)) allocate (GBmat(nres6,nres2))
210 ! if (.not.allocated(Tmat)) allocate (Tmat(nres6,nres2))
211 ! if (.not.allocated(Cmat_)) allocate(Cmat_(nres2,nres2))
212 ! if (.not.allocated(Cinv)) allocate (Cinv(nres2,nres2))
213 ! if (.not.allocated(Pmat)) allocate(Pmat(6*nres,6*nres))
215 if (.not.allocated(stochforcvec)) allocate(stochforcvec(nres6)) !(MAXRES6) maxres6=6*maxres
219 if (itype(i,1).ne.10) nbond=nbond+1
223 write (iout,*) "Generalized inverse of fricmat"
224 call matout(dimen,dimen,nres6,nres6,fricmat)
236 Bmat(ind+j,ind1)=dC_norm(j,i)
241 if (itype(i,1).ne.10) then
244 Bmat(ind+j,ind1)=dC_norm(j,i+nres)
250 write (iout,*) "Matrix Bmat"
251 call MATOUT(nbond,dimen,nres6,nres6,Bmat)
257 GBmat(i,j)=GBmat(i,j)+fricmat(i,k)*Bmat(k,j)
262 write (iout,*) "Matrix GBmat"
263 call MATOUT(nbond,dimen,nres6,nres2,Gbmat)
269 Cmat_(i,j)=Cmat_(i,j)+Bmat(k,i)*GBmat(k,j)
274 write (iout,*) "Matrix Cmat"
275 call MATOUT(nbond,nbond,nres2,nres2,Cmat_)
277 call matinvert(nbond,nres2,Cmat_,Cinv,osob)
279 write (iout,*) "Matrix Cinv"
280 call MATOUT(nbond,nbond,nres2,nres2,Cinv)
286 Tmat(i,j)=Tmat(i,j)+GBmat(i,k)*Cinv(k,j)
291 write (iout,*) "Matrix Tmat"
292 call MATOUT(nbond,dimen,nres6,nres2,Tmat)
302 Pmat(i,j)=Pmat(i,j)-Tmat(i,k)*Bmat(j,k)
307 write (iout,*) "Matrix Pmat"
308 call MATOUT(dimen,dimen,nres6,nres6,Pmat)
315 Td(i)=Td(i)+vbl*Tmat(i,ind)
318 if (itype(k,1).ne.10) then
320 Td(i)=Td(i)+vbldsc0(1,itype(k,1))*Tmat(i,ind)
325 write (iout,*) "Vector Td"
327 write (iout,'(i5,f10.5)') i,Td(i)
330 call stochastic_force(stochforcvec)
332 write (iout,*) "stochforcvec"
334 write (iout,*) i,stochforcvec(i)
338 zapas(j)=-gcart(j,0)+stochforcvec(j)
340 dC_work(j)=dC_old(j,0)
346 zapas(ind)=-gcart(j,i)+stochforcvec(ind)
347 dC_work(ind)=dC_old(j,i)
351 if (itype(i,1).ne.10) then
354 zapas(ind)=-gxcart(j,i)+stochforcvec(ind)
355 dC_work(ind)=dC_old(j,i+nres)
361 write (iout,*) "Initial d_t_work"
363 write (iout,*) i,d_t_work(i)
370 d_t_work(i)=d_t_work(i)+fricmat(i,j)*zapas(j)
377 zapas(i)=zapas(i)+Pmat(i,j)*(dC_work(j)+d_t_work(j)*d_time)
381 write (iout,*) "Final d_t_work and zapas"
383 write (iout,*) i,d_t_work(i),zapas(i)
397 dc_work(ind+j)=dc(j,i)
403 d_t(j,i+nres)=d_t_work(ind+j)
404 dc(j,i+nres)=zapas(ind+j)
405 dc_work(ind+j)=dc(j,i+nres)
411 write (iout,*) "Before correction for rotational lengthening"
412 write (iout,*) "New coordinates",&
413 " and differences between actual and standard bond lengths"
418 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
422 if (itype(i,1).ne.10) then
424 xx=vbld(i+nres)-vbldsc0(1,itype(i,1))
425 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
426 i,(dC(j,i+nres),j=1,3),xx
430 ! Second correction (rotational lengthening)
436 blen2 = scalar(dc(1,i),dc(1,i))
437 ppvec(ind)=2*vbl**2-blen2
438 diffbond=dabs(vbl-dsqrt(blen2))
439 if (diffbond.gt.diffmax) diffmax=diffbond
440 if (ppvec(ind).gt.0.0d0) then
441 ppvec(ind)=dsqrt(ppvec(ind))
446 write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind)
450 if (itype(i,1).ne.10) then
452 blen2 = scalar(dc(1,i+nres),dc(1,i+nres))
453 ppvec(ind)=2*vbldsc0(1,itype(i,1))**2-blen2
454 diffbond=dabs(vbldsc0(1,itype(i,1))-dsqrt(blen2))
455 if (diffbond.gt.diffmax) diffmax=diffbond
456 if (ppvec(ind).gt.0.0d0) then
457 ppvec(ind)=dsqrt(ppvec(ind))
462 write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind)
466 if (lprn) write (iout,*) "iter",iter," diffmax",diffmax
467 if (diffmax.lt.difftol) goto 10
471 Td(i)=Td(i)+ppvec(j)*Tmat(i,j)
477 zapas(i)=zapas(i)+Pmat(i,j)*dc_work(j)
488 dc_work(ind+j)=zapas(ind+j)
493 if (itype(i,1).ne.10) then
495 dc(j,i+nres)=zapas(ind+j)
496 dc_work(ind+j)=zapas(ind+j)
501 ! Building the chain from the newly calculated coordinates
504 if (large.and. mod(itime,ntwe).eq.0) then
505 write (iout,*) "Cartesian and internal coordinates: step 1"
508 write (iout,'(a)') "Potential forces"
510 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(-gcart(j,i),j=1,3),&
513 write (iout,'(a)') "Stochastic forces"
515 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(stochforc(j,i),j=1,3),&
516 (stochforc(j,i+nres),j=1,3)
518 write (iout,'(a)') "Velocities"
520 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
521 (d_t(j,i+nres),j=1,3)
526 write (iout,*) "After correction for rotational lengthening"
527 write (iout,*) "New coordinates",&
528 " and differences between actual and standard bond lengths"
533 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
537 if (itype(i,1).ne.10) then
539 xx=vbld(i+nres)-vbldsc0(1,itype(i,1))
540 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
541 i,(dC(j,i+nres),j=1,3),xx
546 ! write (iout,*) "Too many attempts at correcting the bonds"
554 ! Calculate energy and forces
556 call etotal(potEcomp)
557 potE=potEcomp(0)-potEcomp(51)
561 ! Calculate the kinetic and total energy and the kinetic temperature
564 t_enegrad=t_enegrad+MPI_Wtime()-tt0
566 t_enegrad=t_enegrad+tcpu()-tt0
569 kinetic_T=2.0d0/(dimen*Rb)*EK
571 end subroutine brown_step
572 !-----------------------------------------------------------------------------
574 !-----------------------------------------------------------------------------
575 subroutine gauss(RO,AP,MT,M,N,*)
577 ! CALCULATES (RO**(-1))*AP BY GAUSS ELIMINATION
578 ! RO IS A SQUARE MATRIX
579 ! THE CALCULATED PRODUCT IS STORED IN AP
580 ! ABNORMAL EXIT IF RO IS SINGULAR
582 integer :: MT, M, N, M1,I,J,IM,&
584 real(kind=8) :: RO(MT,M),AP(MT,N),X,RM,PR,Y
590 if(dabs(X).le.1.0D-13) return 1
602 if(DABS(RO(J,I)).LE.RM) goto 2
616 if(dabs(X).le.1.0E-13) return 1
625 8 AP(J,K)=AP(J,K)-Y*AP(I,K)
627 9 RO(J,K)=RO(J,K)-Y*RO(I,K)
631 if(dabs(X).le.1.0E-13) return 1
641 15 X=X-AP(K,J)*RO(MI,K)
646 !-----------------------------------------------------------------------------
649 subroutine kinetic(KE_total)
650 !c----------------------------------------------------------------
651 !c This subroutine calculates the total kinetic energy of the chain
652 !c-----------------------------------------------------------------
653 !c 3/5/2020 AL Corrected for multichain systems, no fake peptide groups
654 !c inside, implemented with five-diagonal inertia matrix
658 include 'COMMON.CHAIN'
659 include 'COMMON.DERIV'
661 include 'COMMON.LOCAL'
662 include 'COMMON.INTERACT'
664 include 'COMMON.LAGRANGE.5diag'
665 include 'COMMON.IOUNITS'
666 double precision KE_total
668 double precision KEt_p,KEt_sc,KEr_p,KEr_sc,incr(3),&
675 !c write (iout,*) "ISC",(isc(itype(i)),i=1,nres)
676 !c The translational part for peptide virtual bonds
681 !c write (iout,*) "Kinetic trp:",i,(incr(j),j=1,3
682 !c Skip dummy peptide groups
683 if (itype(i).ne.ntyp1 .and. itype(i+1).ne.ntyp1) then
685 v(j)=incr(j)+0.5d0*d_t(j,i)
687 !c write (iout,*) "Kinetic trp:",i,(v(j),j=1,3)
688 vtot(i)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
689 KEt_p=KEt_p+(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
692 incr(j)=incr(j)+d_t(j,i)
695 !c write(iout,*) 'KEt_p', KEt_p
696 !c The translational part for the side chain virtual bond
697 !c Only now we can initialize incr with zeros. It must be equal
698 !c to the velocities of the first Calpha.
704 if (itype(i).eq.10 .and. itype(i).ne.ntyp1) then
710 v(j)=incr(j)+d_t(j,nres+i)
713 !c write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3)
714 !c write (iout,*) "i",i," msc",msc(iti)," v",(v(j),j=1,3)
715 KEt_sc=KEt_sc+msc(iti)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
716 vtot(i+nres)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
718 incr(j)=incr(j)+d_t(j,i)
722 ! write(iout,*) 'KEt_sc', KEt_sc
723 ! The part due to stretching and rotation of the peptide groups
725 if (itype(i).ne.ntyp1.and.itype(i+1).ne.ntyp1) then
726 ! write (iout,*) "i",i
727 ! write (iout,*) "i",i," mag1",mag1," mag2",mag2
731 !c write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3)
732 KEr_p=KEr_p+(incr(1)*incr(1)+incr(2)*incr(2) &
737 !c write(iout,*) 'KEr_p', KEr_p
738 !c The rotational part of the side chain virtual bond
741 if (itype(i).ne.10.and.itype(i).ne.ntyp1) then
743 incr(j)=d_t(j,nres+i)
745 !c write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3)
746 KEr_sc=KEr_sc+Isc(iti)*(incr(1)*incr(1)+incr(2)*incr(2)+&
750 !c The total kinetic energy
752 !c write(iout,*) ' KEt_p',KEt_p,' KEt_sc',KEt_sc,' KEr_p',KEr_p,
753 !c & ' KEr_sc', KEr_sc
754 KE_total=0.5d0*(mp*KEt_p+KEt_sc+0.25d0*Ip*KEr_p+KEr_sc)
755 !c write (iout,*) "KE_total",KE_total
757 end subroutine kinetic
760 !-----------------------------------------------------------------------------
761 subroutine kinetic(KE_total)
762 !----------------------------------------------------------------
763 ! This subroutine calculates the total kinetic energy of the chain
764 !-----------------------------------------------------------------
766 ! implicit real*8 (a-h,o-z)
767 ! include 'DIMENSIONS'
768 ! include 'COMMON.VAR'
769 ! include 'COMMON.CHAIN'
770 ! include 'COMMON.DERIV'
771 ! include 'COMMON.GEO'
772 ! include 'COMMON.LOCAL'
773 ! include 'COMMON.INTERACT'
774 ! include 'COMMON.MD'
775 ! include 'COMMON.IOUNITS'
776 real(kind=8) :: KE_total,mscab
778 integer :: i,j,k,iti,mnum,term
779 real(kind=8) :: KEt_p,KEt_sc,KEr_p,KEr_sc,incr(3),&
782 write (iout,*) "Velocities, kietic"
784 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
785 (d_t(j,i+nres),j=1,3)
790 ! write (iout,*) "ISC",(isc(itype(i,1)),i=1,nres)
791 ! The translational part for peptide virtual bonds
796 ! if (molnum(nct).gt.3) term=nct
799 if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
800 ! write (iout,*) "Kinetic trp:",i,(incr(j),j=1,3),mp(mnum)
803 v(j)=incr(j)+0.5d0*d_t(j,i)
807 v(j)=incr(j)+0.5d0*d_t(j,i)
810 vtot(i)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
811 KEt_p=KEt_p+mp(mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
813 incr(j)=incr(j)+d_t(j,i)
816 ! write(iout,*) 'KEt_p', KEt_p
817 ! The translational part for the side chain virtual bond
818 ! Only now we can initialize incr with zeros. It must be equal
819 ! to the velocities of the first Calpha.
825 iti=iabs(itype(i,mnum))
826 ! if (mnum.ge.4) then
831 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
832 if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
833 .or.(mnum.ge.4)) then
839 v(j)=incr(j)+d_t(j,nres+i)
842 ! write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3)
843 ! write (iout,*) "i",i," msc",msc(iti,mnum)," v",(v(j),j=1,3)
844 KEt_sc=KEt_sc+mscab*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
845 vtot(i+nres)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
847 incr(j)=incr(j)+d_t(j,i)
851 ! write(iout,*) 'KEt_sc', KEt_sc
852 ! The part due to stretching and rotation of the peptide groups
856 ! write (iout,*) "i",i
857 ! write (iout,*) "i",i," mag1",mag1," mag2",mag2
861 ! write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3)
862 KEr_p=KEr_p+Ip(mnum)*(incr(1)*incr(1)+incr(2)*incr(2) &
866 ! write(iout,*) 'KEr_p', KEr_p
867 ! The rotational part of the side chain virtual bond
871 iti=iabs(itype(i,mnum))
872 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
873 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
874 .and.(mnum.lt.4)) then
876 incr(j)=d_t(j,nres+i)
878 ! write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3)
879 KEr_sc=KEr_sc+Isc(iti,mnum)*(incr(1)*incr(1)+incr(2)*incr(2)+ &
883 ! The total kinetic energy
885 ! write(iout,*) 'KEr_sc', KEr_sc
886 KE_total=0.5d0*(KEt_p+KEt_sc+0.25d0*KEr_p+KEr_sc)
887 ! write (iout,*) "KE_total",KE_total
889 end subroutine kinetic
890 !-----------------------------------------------------------------------------
892 subroutine kinetic_CASC(KE_total)
893 !c----------------------------------------------------------------
894 !c Compute the kinetic energy of the system using the Calpha-SC
896 !c-----------------------------------------------------------------
898 double precision KE_total
900 integer i,j,k,iti,ichain,innt,inct
901 double precision KEt_p,KEt_sc,KEr_p,KEr_sc,incr(3),&
908 !c write (iout,*) "ISC",(isc(itype(i)),i=1,nres)
909 !c The translational part for peptide virtual bonds
912 innt=chain_border(1,ichain)
913 inct=chain_border(2,ichain)
914 !c write (iout,*) "Kinetic_CASC chain",ichain," innt",innt,
918 !c write (iout,*) i,(d_t(j,i),j=1,3),(d_t(j,i+1),j=1,3)
920 v(j)=0.5d0*(d_t(j,i)+d_t(j,i+1))
922 !c write (iout,*) "Kinetic trp i",i," v",(v(j),j=1,3)
923 KEt_p=KEt_p+(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
925 !c write(iout,*) 'KEt_p', KEt_p
926 !c The translational part for the side chain virtual bond
927 !c Only now we can initialize incr with zeros. It must be equal
928 !c to the velocities of the first Calpha.
932 !c write (iout,*) i,iti,(d_t(j,i),j=1,3)
937 !c write (iout,*) i,iti,(d_t(j,nres+i),j=1,3)
942 !c write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3)
943 !c write (iout,*) "i",i," msc",msc(iti)," v",(v(j),j=1,3)
944 KEt_sc=KEt_sc+msc(iti)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
947 !c write(iout,*) 'KEt_sc', KEt_sc
948 !c The part due to stretching and rotation of the peptide groups
951 incr(j)=d_t(j,i+1)-d_t(j,i)
953 !c write (iout,*) i,(incr(j),j=1,3)
954 !c write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3)
955 KEr_p=KEr_p+(incr(1)*incr(1)+incr(2)*incr(2)
959 !c write(iout,*) 'KEr_p', KEr_p
960 !c The rotational part of the side chain virtual bond
965 incr(j)=d_t(j,nres+i)-d_t(j,i)
967 !c write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3)
968 KEr_sc=KEr_sc+Isc(iti)*(incr(1)*incr(1)+incr(2)*incr(2)+
974 !c The total kinetic energy
976 !c write(iout,*) ' KEt_p',KEt_p,' KEt_sc',KEt_sc,' KEr_p',KEr_p,
977 !c & ' KEr_sc', KEr_sc
978 KE_total=0.5d0*(mp*KEt_p+KEt_sc+0.25d0*Ip*KEr_p+KEr_sc)
979 !c write (iout,*) "KE_total",KE_tota
981 write (iout,*) "Need to compile with -DFIVEDIAG to use this sub!"
985 end subroutine kinetic_CASC
988 !-----------------------------------------------------------------------------
990 !------------------------------------------------
991 ! The driver for molecular dynamics subroutines
992 !------------------------------------------------
995 use control, only:tcpu,ovrtim
996 ! use io_comm, only:ilen
998 use compare, only:secondary2,hairpin
999 use io, only:cartout,statout
1000 ! implicit real*8 (a-h,o-z)
1001 ! include 'DIMENSIONS'
1004 integer :: IERROR,ERRCODE
1006 ! include 'COMMON.SETUP'
1007 ! include 'COMMON.CONTROL'
1008 ! include 'COMMON.VAR'
1009 ! include 'COMMON.MD'
1011 ! include 'COMMON.LANGEVIN'
1013 ! include 'COMMON.LANGEVIN.lang0'
1015 ! include 'COMMON.CHAIN'
1016 ! include 'COMMON.DERIV'
1017 ! include 'COMMON.GEO'
1018 ! include 'COMMON.LOCAL'
1019 ! include 'COMMON.INTERACT'
1020 ! include 'COMMON.IOUNITS'
1021 ! include 'COMMON.NAMES'
1022 ! include 'COMMON.TIME1'
1023 ! include 'COMMON.HAIRPIN'
1024 real(kind=8),dimension(3) :: L,vcm
1026 real(kind=8),dimension(6*nres) :: v_work,v_transf !(maxres6) maxres6=6*maxres
1028 integer :: rstcount !ilen,
1030 character(len=50) :: tytul
1031 !el common /gucio/ cm
1032 integer :: i,j,nharp
1033 integer,dimension(4,nres) :: iharp !(4,nres/3)(4,maxres/3)
1035 real(kind=8) :: tt0,scalfac
1036 integer :: nres2,itime
1041 print *,"MY tmpdir",tmpdir,ilen(tmpdir)
1042 if (ilen(tmpdir).gt.0) &
1043 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_" &
1044 //liczba(:ilen(liczba))//'.rst')
1046 if (ilen(tmpdir).gt.0) &
1047 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"//'.rst')
1054 write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
1060 print *,"just befor setup matix",nres
1061 ! Determine the inverse of the inertia matrix.
1062 call setup_MD_matrices
1064 print *,"AFTER SETUP MATRICES"
1066 print *,"AFTER INIT MD"
1069 t_MDsetup = MPI_Wtime()-tt0
1071 t_MDsetup = tcpu()-tt0
1074 ! Entering the MD loop
1080 if (lang.eq.2 .or. lang.eq.3) then
1084 call sd_verlet_p_setup
1086 call sd_verlet_ciccotti_setup
1090 pfric0_mat(i,j,0)=pfric_mat(i,j)
1091 afric0_mat(i,j,0)=afric_mat(i,j)
1092 vfric0_mat(i,j,0)=vfric_mat(i,j)
1093 prand0_mat(i,j,0)=prand_mat(i,j)
1094 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
1095 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
1098 flag_stoch(0)=.true.
1099 do i=1,maxflag_stoch
1100 flag_stoch(i)=.false.
1104 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
1106 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1110 else if (lang.eq.1 .or. lang.eq.4) then
1111 print *,"before setup_fricmat"
1113 print *,"after setup_fricmat"
1116 t_langsetup=MPI_Wtime()-tt0
1119 t_langsetup=tcpu()-tt0
1122 do itime=1,n_timestep
1124 if (large.and. mod(itime,ntwe).eq.0) &
1125 write (iout,*) "itime",itime
1127 if (lang.gt.0 .and. surfarea .and. &
1128 mod(itime,reset_fricmat).eq.0) then
1129 if (lang.eq.2 .or. lang.eq.3) then
1133 call sd_verlet_p_setup
1135 call sd_verlet_ciccotti_setup
1139 pfric0_mat(i,j,0)=pfric_mat(i,j)
1140 afric0_mat(i,j,0)=afric_mat(i,j)
1141 vfric0_mat(i,j,0)=vfric_mat(i,j)
1142 prand0_mat(i,j,0)=prand_mat(i,j)
1143 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
1144 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
1147 flag_stoch(0)=.true.
1148 do i=1,maxflag_stoch
1149 flag_stoch(i)=.false.
1152 else if (lang.eq.1 .or. lang.eq.4) then
1155 write (iout,'(a,i10)') &
1156 "Friction matrix reset based on surface area, itime",itime
1158 if (reset_vel .and. tbf .and. lang.eq.0 &
1159 .and. mod(itime,count_reset_vel).eq.0) then
1161 write(iout,'(a,f20.2)') &
1162 "Velocities reset to random values, time",totT
1165 d_t_old(j,i)=d_t(j,i)
1169 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
1173 d_t(j,0)=d_t(j,0)-vcm(j)
1176 kinetic_T=2.0d0/(dimen3*Rb)*EK
1177 scalfac=dsqrt(T_bath/kinetic_T)
1178 write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
1181 d_t_old(j,i)=scalfac*d_t(j,i)
1187 ! Time-reversible RESPA algorithm
1188 ! (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
1189 call RESPA_step(itime)
1191 ! Variable time step algorithm.
1192 call velverlet_step(itime)
1196 call brown_step(itime)
1198 print *,"Brown dynamics not here!"
1200 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1207 if (mod(itime,ntwe).eq.0) then
1211 ! call check_ecartint
1221 v_work(ind)=d_t(j,i)
1226 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1.and.mnum.lt.4) then
1229 v_work(ind)=d_t(j,i+nres)
1234 write (66,'(80f10.5)') &
1235 ((d_t(j,i),j=1,3),i=0,nres-1),((d_t(j,i+nres),j=1,3),i=1,nres)
1239 v_transf(i)=v_transf(i)+gvec(j,i)*v_work(j)
1241 v_transf(i)= v_transf(i)*dsqrt(geigen(i))
1243 write (67,'(80f10.5)') (v_transf(i),i=1,ind)
1246 if (mod(itime,ntwx).eq.0) then
1248 write (tytul,'("time",f8.2)') totT
1250 call hairpin(.true.,nharp,iharp)
1251 call secondary2(.true.)
1252 call pdbout(potE,tytul,ipdb)
1253 call enerprint(potEcomp)
1258 write(iout,*) "starting fodstep"
1259 call fodstep(nfodstep)
1260 write(iout,*) "after fodstep"
1263 call hairpin(.true.,nharp,iharp)
1264 call secondary2(.true.)
1265 call pdbout(potE,tytul,ipdb)
1272 if (rstcount.eq.1000.or.itime.eq.n_timestep) then
1273 open(irest2,file=rest2name,status='unknown')
1274 write(irest2,*) totT,EK,potE,totE,t_bath
1276 ! AL 4/17/17: Now writing d_t(0,:) too
1278 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
1280 ! AL 4/17/17: Now writing d_c(0,:) too
1282 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
1290 t_MD=MPI_Wtime()-tt0
1294 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))') &
1296 'MD calculations setup:',t_MDsetup,&
1297 'Energy & gradient evaluation:',t_enegrad,&
1298 'Stochastic MD setup:',t_langsetup,&
1299 'Stochastic MD step setup:',t_sdsetup,&
1301 write (iout,'(/28(1h=),a25,27(1h=))') &
1302 ' End of MD calculation '
1304 write (iout,*) "time for etotal",t_etotal," elong",t_elong,&
1306 write (iout,*) "time_fric",time_fric," time_stoch",time_stoch,&
1307 " time_fricmatmult",time_fricmatmult," time_fsample ",&
1312 !-----------------------------------------------------------------------------
1313 subroutine velverlet_step(itime)
1314 !-------------------------------------------------------------------------------
1315 ! Perform a single velocity Verlet step; the time step can be rescaled if
1316 ! increments in accelerations exceed the threshold
1317 !-------------------------------------------------------------------------------
1318 ! implicit real*8 (a-h,o-z)
1319 ! include 'DIMENSIONS'
1321 use control, only:tcpu
1323 use minimm, only:minim_dc
1326 integer :: ierror,ierrcode
1327 real(kind=8) :: errcode
1329 ! include 'COMMON.SETUP'
1330 ! include 'COMMON.VAR'
1331 ! include 'COMMON.MD'
1333 ! include 'COMMON.LANGEVIN'
1335 ! include 'COMMON.LANGEVIN.lang0'
1337 ! include 'COMMON.CHAIN'
1338 ! include 'COMMON.DERIV'
1339 ! include 'COMMON.GEO'
1340 ! include 'COMMON.LOCAL'
1341 ! include 'COMMON.INTERACT'
1342 ! include 'COMMON.IOUNITS'
1343 ! include 'COMMON.NAMES'
1344 ! include 'COMMON.TIME1'
1345 ! include 'COMMON.MUCA'
1346 real(kind=8),dimension(3) :: vcm,incr
1347 real(kind=8),dimension(3) :: L
1348 integer :: count,rstcount !ilen,
1350 character(len=50) :: tytul
1351 integer :: maxcount_scale = 30
1352 !el common /gucio/ cm
1353 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1354 !el common /stochcalc/ stochforcvec
1355 integer :: icount_scale,itime_scal,i,j,ifac_time,iretcode,itime
1357 real(kind=8) :: epdrift,tt0,fac_time
1359 if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres)) !(MAXRES6) maxres6=6*maxres
1365 else if (lang.eq.2 .or. lang.eq.3) then
1367 call stochastic_force(stochforcvec)
1370 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
1372 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1379 icount_scale=icount_scale+1
1380 ! write(iout,*) "icount_scale",icount_scale,scalel
1381 if (icount_scale.gt.maxcount_scale) then
1383 "ERROR: too many attempts at scaling down the time step. ",&
1384 "amax=",amax,"epdrift=",epdrift,&
1385 "damax=",damax,"edriftmax=",edriftmax,&
1389 call MPI_Abort(MPI_COMM_WORLD,IERROR,IERRCODE)
1393 ! First step of the velocity Verlet algorithm
1398 else if (lang.eq.3) then
1400 call sd_verlet1_ciccotti
1402 else if (lang.eq.1) then
1407 ! Build the chain from the newly calculated coordinates
1408 call chainbuild_cart
1409 if (rattle) call rattle1
1411 if (large.and. mod(itime,ntwe).eq.0) then
1412 write (iout,*) "Cartesian and internal coordinates: step 1"
1417 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),&
1418 (dc(j,i+nres),j=1,3)
1420 write (iout,*) "Accelerations"
1422 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1423 (d_a(j,i+nres),j=1,3)
1425 write (iout,*) "Velocities, step 1"
1427 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1428 (d_t(j,i+nres),j=1,3)
1437 ! Calculate energy and forces
1439 call etotal(potEcomp)
1440 ! AL 4/17/17: Reduce the steps if NaNs occurred.
1441 if (potEcomp(0).gt.0.99e18 .or. isnan(potEcomp(0)).gt.0) then
1442 call enerprint(potEcomp)
1444 if (icount_scale.gt.15) then
1445 write (iout,*) "Tu jest problem",potEcomp(0),d_time
1446 ! call gen_rand_conf(1,*335)
1447 ! call minim_dc(potEcomp(0),iretcode,100)
1450 ! call etotal(potEcomp)
1451 ! write(iout,*) "needed to repara,",potEcomp
1454 ! 335 write(iout,*) "Failed genrand"
1458 if (large.and. mod(itime,ntwe).eq.0) &
1459 call enerprint(potEcomp)
1462 t_etotal=t_etotal+MPI_Wtime()-tt0
1464 t_etotal=t_etotal+tcpu()-tt0
1467 potE=potEcomp(0)-potEcomp(51)
1469 ! Get the new accelerations
1472 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1474 t_enegrad=t_enegrad+tcpu()-tt0
1476 ! Determine maximum acceleration and scale down the timestep if needed
1478 amax=amax/(itime_scal+1)**2
1479 call predict_edrift(epdrift)
1480 ! write(iout,*) "amax=",amax,damax,epdrift,edriftmax,amax/(itime_scal+1)
1482 ! write (iout,*) "before enter if",scalel,icount_scale
1483 if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then
1484 ! write(iout,*) "I enter if"
1485 ! Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
1487 ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax)) &
1489 itime_scal=itime_scal+ifac_time
1490 ! fac_time=dmin1(damax/amax,0.5d0)
1491 fac_time=0.5d0**ifac_time
1492 d_time=d_time*fac_time
1493 if (lang.eq.2 .or. lang.eq.3) then
1495 ! write (iout,*) "Calling sd_verlet_setup: 1"
1496 ! Rescale the stochastic forces and recalculate or restore
1497 ! the matrices of tinker integrator
1498 if (itime_scal.gt.maxflag_stoch) then
1499 if (large) write (iout,'(a,i5,a)') &
1500 "Calculate matrices for stochastic step;",&
1501 " itime_scal ",itime_scal
1503 call sd_verlet_p_setup
1505 call sd_verlet_ciccotti_setup
1507 write (iout,'(2a,i3,a,i3,1h.)') &
1508 "Warning: cannot store matrices for stochastic",&
1509 " integration because the index",itime_scal,&
1510 " is greater than",maxflag_stoch
1511 write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",&
1512 " integration Langevin algorithm for better efficiency."
1513 else if (flag_stoch(itime_scal)) then
1514 if (large) write (iout,'(a,i5,a,l1)') &
1515 "Restore matrices for stochastic step; itime_scal ",&
1516 itime_scal," flag ",flag_stoch(itime_scal)
1519 pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
1520 afric_mat(i,j)=afric0_mat(i,j,itime_scal)
1521 vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
1522 prand_mat(i,j)=prand0_mat(i,j,itime_scal)
1523 vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
1524 vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
1528 if (large) write (iout,'(2a,i5,a,l1)') &
1529 "Calculate & store matrices for stochastic step;",&
1530 " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
1532 call sd_verlet_p_setup
1534 call sd_verlet_ciccotti_setup
1536 flag_stoch(ifac_time)=.true.
1539 pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
1540 afric0_mat(i,j,itime_scal)=afric_mat(i,j)
1541 vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
1542 prand0_mat(i,j,itime_scal)=prand_mat(i,j)
1543 vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
1544 vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
1548 fac_time=1.0d0/dsqrt(fac_time)
1550 stochforcvec(i)=fac_time*stochforcvec(i)
1553 else if (lang.eq.1) then
1554 ! Rescale the accelerations due to stochastic forces
1555 fac_time=1.0d0/dsqrt(fac_time)
1557 d_as_work(i)=d_as_work(i)*fac_time
1560 if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)') &
1561 "itime",itime," Timestep scaled down to ",&
1562 d_time," ifac_time",ifac_time," itime_scal",itime_scal
1564 ! Second step of the velocity Verlet algorithm
1569 else if (lang.eq.3) then
1571 call sd_verlet2_ciccotti
1573 else if (lang.eq.1) then
1578 if (rattle) call rattle2
1581 if (d_time.ne.d_time0) then
1584 if (lang.eq.2 .or. lang.eq.3) then
1585 if (large) write (iout,'(a)') &
1586 "Restore original matrices for stochastic step"
1587 ! write (iout,*) "Calling sd_verlet_setup: 2"
1588 ! Restore the matrices of tinker integrator if the time step has been restored
1591 pfric_mat(i,j)=pfric0_mat(i,j,0)
1592 afric_mat(i,j)=afric0_mat(i,j,0)
1593 vfric_mat(i,j)=vfric0_mat(i,j,0)
1594 prand_mat(i,j)=prand0_mat(i,j,0)
1595 vrand_mat1(i,j)=vrand0_mat1(i,j,0)
1596 vrand_mat2(i,j)=vrand0_mat2(i,j,0)
1604 ! Calculate the kinetic and the total energy and the kinetic temperature
1608 ! call kinetic1(EK1)
1609 ! write (iout,*) "step",itime," EK",EK," EK1",EK1
1611 ! Couple the system to Berendsen bath if needed
1612 if (tbf .and. lang.eq.0) then
1615 kinetic_T=2.0d0/(dimen3*Rb)*EK
1616 ! Backup the coordinates, velocities, and accelerations
1620 d_t_old(j,i)=d_t(j,i)
1621 d_a_old(j,i)=d_a(j,i)
1625 if (mod(itime,ntwe).eq.0 .and. large) then
1626 write (iout,*) "Velocities, step 2"
1628 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1629 (d_t(j,i+nres),j=1,3)
1634 end subroutine velverlet_step
1635 !-----------------------------------------------------------------------------
1636 subroutine RESPA_step(itime)
1637 !-------------------------------------------------------------------------------
1638 ! Perform a single RESPA step.
1639 !-------------------------------------------------------------------------------
1640 ! implicit real*8 (a-h,o-z)
1641 ! include 'DIMENSIONS'
1645 use control, only:tcpu
1647 ! use io_conf, only:cartprint
1650 integer :: IERROR,ERRCODE
1652 ! include 'COMMON.SETUP'
1653 ! include 'COMMON.CONTROL'
1654 ! include 'COMMON.VAR'
1655 ! include 'COMMON.MD'
1657 ! include 'COMMON.LANGEVIN'
1659 ! include 'COMMON.LANGEVIN.lang0'
1661 ! include 'COMMON.CHAIN'
1662 ! include 'COMMON.DERIV'
1663 ! include 'COMMON.GEO'
1664 ! include 'COMMON.LOCAL'
1665 ! include 'COMMON.INTERACT'
1666 ! include 'COMMON.IOUNITS'
1667 ! include 'COMMON.NAMES'
1668 ! include 'COMMON.TIME1'
1669 real(kind=8),dimension(0:n_ene) :: energia_short,energia_long
1670 real(kind=8),dimension(3) :: L,vcm,incr
1671 real(kind=8),dimension(3,0:2*nres) :: dc_old0,d_t_old0,d_a_old0 !(3,0:maxres2) maxres2=2*maxres
1672 logical :: PRINT_AMTS_MSG = .false.
1673 integer :: count,rstcount !ilen,
1675 character(len=50) :: tytul
1676 integer :: maxcount_scale = 10
1677 !el common /gucio/ cm
1678 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1679 !el common /stochcalc/ stochforcvec
1680 integer :: itt,i,j,itsplit,itime
1682 !el common /cipiszcze/ itt
1684 real(kind=8) :: epdrift,tt0,epdriftmax
1687 if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres)) !(MAXRES6) maxres6=6*maxres
1691 if (large.and. mod(itime,ntwe).eq.0) then
1692 write (iout,*) "***************** RESPA itime",itime
1693 write (iout,*) "Cartesian and internal coordinates: step 0"
1695 call pdbout(0.0d0,"cipiszcze",iout)
1697 write (iout,*) "Accelerations from long-range forces"
1699 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1700 (d_a(j,i+nres),j=1,3)
1702 write (iout,*) "Velocities, step 0"
1704 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1705 (d_t(j,i+nres),j=1,3)
1710 ! Perform the initial RESPA step (increment velocities)
1711 ! write (iout,*) "*********************** RESPA ini"
1714 if (mod(itime,ntwe).eq.0 .and. large) then
1715 write (iout,*) "Velocities, end"
1717 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1718 (d_t(j,i+nres),j=1,3)
1722 ! Compute the short-range forces
1728 ! 7/2/2009 commented out
1730 ! call etotal_short(energia_short)
1733 ! 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
1736 d_a(j,i)=d_a_short(j,i)
1740 if (large.and. mod(itime,ntwe).eq.0) then
1741 write (iout,*) "energia_short",energia_short(0)
1742 write (iout,*) "Accelerations from short-range forces"
1744 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1745 (d_a(j,i+nres),j=1,3)
1750 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1752 t_enegrad=t_enegrad+tcpu()-tt0
1757 d_t_old(j,i)=d_t(j,i)
1758 d_a_old(j,i)=d_a(j,i)
1761 ! 6/30/08 A-MTS: attempt at increasing the split number
1764 dc_old0(j,i)=dc_old(j,i)
1765 d_t_old0(j,i)=d_t_old(j,i)
1766 d_a_old0(j,i)=d_a_old(j,i)
1769 if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
1770 if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
1777 ! write (iout,*) "itime",itime," ntime_split",ntime_split
1778 ! Split the time step
1779 d_time=d_time0/ntime_split
1780 ! Perform the short-range RESPA steps (velocity Verlet increments of
1781 ! positions and velocities using short-range forces)
1782 ! write (iout,*) "*********************** RESPA split"
1783 do itsplit=1,ntime_split
1786 else if (lang.eq.2 .or. lang.eq.3) then
1788 call stochastic_force(stochforcvec)
1791 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
1793 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1798 ! First step of the velocity Verlet algorithm
1803 else if (lang.eq.3) then
1805 call sd_verlet1_ciccotti
1807 else if (lang.eq.1) then
1812 ! Build the chain from the newly calculated coordinates
1813 call chainbuild_cart
1814 if (rattle) call rattle1
1816 if (large.and. mod(itime,ntwe).eq.0) then
1817 write (iout,*) "***** ITSPLIT",itsplit
1818 write (iout,*) "Cartesian and internal coordinates: step 1"
1819 call pdbout(0.0d0,"cipiszcze",iout)
1822 write (iout,*) "Velocities, step 1"
1824 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1825 (d_t(j,i+nres),j=1,3)
1834 ! Calculate energy and forces
1836 call etotal_short(energia_short)
1837 ! AL 4/17/17: Exit itime_split loop when energy goes infinite
1838 if (energia_short(0).gt.0.99e20 .or. isnan(energia_short(0)) ) then
1839 if (PRINT_AMTS_MSG) &
1840 write (iout,*) "Infinities/NaNs in energia_short",energia_short(0),"; increasing ntime_split to",ntime_split
1841 ntime_split=ntime_split*2
1842 if (ntime_split.gt.maxtime_split) then
1845 "Cannot rescue the run; aborting job. Retry with a smaller time step"
1847 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1850 "Cannot rescue the run; terminating. Retry with a smaller time step"
1856 if (large.and. mod(itime,ntwe).eq.0) &
1857 call enerprint(energia_short)
1860 t_eshort=t_eshort+MPI_Wtime()-tt0
1862 t_eshort=t_eshort+tcpu()-tt0
1866 ! Get the new accelerations
1868 ! 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
1871 d_a_short(j,i)=d_a(j,i)
1875 if (large.and. mod(itime,ntwe).eq.0) then
1876 write (iout,*)"energia_short",energia_short(0)
1877 write (iout,*) "Accelerations from short-range forces"
1879 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1880 (d_a(j,i+nres),j=1,3)
1885 ! Determine maximum acceleration and scale down the timestep if needed
1887 amax=amax/ntime_split**2
1888 call predict_edrift(epdrift)
1889 if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0) &
1890 write (iout,*) "amax",amax," damax",damax,&
1891 " epdrift",epdrift," epdriftmax",epdriftmax
1892 ! Exit loop and try with increased split number if the change of
1893 ! acceleration is too big
1894 if (amax.gt.damax .or. epdrift.gt.edriftmax) then
1895 if (ntime_split.lt.maxtime_split) then
1897 ntime_split=ntime_split*2
1898 ! AL 4/17/17: We should exit the itime_split loop when acceleration change is too big
1902 dc_old(j,i)=dc_old0(j,i)
1903 d_t_old(j,i)=d_t_old0(j,i)
1904 d_a_old(j,i)=d_a_old0(j,i)
1907 if (PRINT_AMTS_MSG) then
1908 write (iout,*) "acceleration/energy drift too large",amax,&
1909 epdrift," split increased to ",ntime_split," itime",itime,&
1915 "Uh-hu. Bumpy landscape. Maximum splitting number",&
1917 " already reached!!! Trying to carry on!"
1921 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1923 t_enegrad=t_enegrad+tcpu()-tt0
1925 ! Second step of the velocity Verlet algorithm
1930 else if (lang.eq.3) then
1932 call sd_verlet2_ciccotti
1934 else if (lang.eq.1) then
1939 if (rattle) call rattle2
1940 ! Backup the coordinates, velocities, and accelerations
1944 d_t_old(j,i)=d_t(j,i)
1945 d_a_old(j,i)=d_a(j,i)
1952 ! Restore the time step
1954 ! Compute long-range forces
1961 call etotal_long(energia_long)
1962 if (energia_long(0).gt.0.99e20 .or. isnan(energia_long(0))) then
1965 "Infinitied/NaNs in energia_long, Aborting MPI job."
1967 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1969 write (iout,*) "Infinitied/NaNs in energia_long, terminating."
1973 if (large.and. mod(itime,ntwe).eq.0) &
1974 call enerprint(energia_long)
1977 t_elong=t_elong+MPI_Wtime()-tt0
1979 t_elong=t_elong+tcpu()-tt0
1982 potE=potEcomp(0)-potEcomp(51)
1986 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1988 t_enegrad=t_enegrad+tcpu()-tt0
1990 ! Compute accelerations from long-range forces
1992 if (large.and. mod(itime,ntwe).eq.0) then
1993 write (iout,*) "energia_long",energia_long(0)
1994 write (iout,*) "Cartesian and internal coordinates: step 2"
1996 call pdbout(0.0d0,"cipiszcze",iout)
1998 write (iout,*) "Accelerations from long-range forces"
2000 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
2001 (d_a(j,i+nres),j=1,3)
2003 write (iout,*) "Velocities, step 2"
2005 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2006 (d_t(j,i+nres),j=1,3)
2010 ! Compute the final RESPA step (increment velocities)
2011 ! write (iout,*) "*********************** RESPA fin"
2013 ! Compute the complete potential energy
2015 potEcomp(i)=energia_short(i)+energia_long(i)
2017 potE=potEcomp(0)-potEcomp(51)
2018 ! potE=energia_short(0)+energia_long(0)
2021 ! Calculate the kinetic and the total energy and the kinetic temperature
2024 ! Couple the system to Berendsen bath if needed
2025 if (tbf .and. lang.eq.0) then
2028 kinetic_T=2.0d0/(dimen3*Rb)*EK
2029 ! Backup the coordinates, velocities, and accelerations
2031 if (mod(itime,ntwe).eq.0 .and. large) then
2032 write (iout,*) "Velocities, end"
2034 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2035 (d_t(j,i+nres),j=1,3)
2040 end subroutine RESPA_step
2041 !-----------------------------------------------------------------------------
2042 subroutine RESPA_vel
2043 ! First and last RESPA step (incrementing velocities using long-range
2046 ! implicit real*8 (a-h,o-z)
2047 ! include 'DIMENSIONS'
2048 ! include 'COMMON.CONTROL'
2049 ! include 'COMMON.VAR'
2050 ! include 'COMMON.MD'
2051 ! include 'COMMON.CHAIN'
2052 ! include 'COMMON.DERIV'
2053 ! include 'COMMON.GEO'
2054 ! include 'COMMON.LOCAL'
2055 ! include 'COMMON.INTERACT'
2056 ! include 'COMMON.IOUNITS'
2057 ! include 'COMMON.NAMES'
2058 integer :: i,j,inres,mnum
2061 d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
2065 d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
2070 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2071 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2072 .and.(mnum.lt.4)) then
2075 d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
2080 end subroutine RESPA_vel
2081 !-----------------------------------------------------------------------------
2083 ! Applying velocity Verlet algorithm - step 1 to coordinates
2085 ! implicit real*8 (a-h,o-z)
2086 ! include 'DIMENSIONS'
2087 ! include 'COMMON.CONTROL'
2088 ! include 'COMMON.VAR'
2089 ! include 'COMMON.MD'
2090 ! include 'COMMON.CHAIN'
2091 ! include 'COMMON.DERIV'
2092 ! include 'COMMON.GEO'
2093 ! include 'COMMON.LOCAL'
2094 ! include 'COMMON.INTERACT'
2095 ! include 'COMMON.IOUNITS'
2096 ! include 'COMMON.NAMES'
2097 real(kind=8) :: adt,adt2
2098 integer :: i,j,inres,mnum
2101 write (iout,*) "VELVERLET1 START: DC"
2103 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),&
2104 (dc(j,i+nres),j=1,3)
2108 adt=d_a_old(j,0)*d_time
2110 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
2111 d_t_new(j,0)=d_t_old(j,0)+adt2
2112 d_t(j,0)=d_t_old(j,0)+adt
2116 adt=d_a_old(j,i)*d_time
2118 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
2119 d_t_new(j,i)=d_t_old(j,i)+adt2
2120 d_t(j,i)=d_t_old(j,i)+adt
2125 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2126 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2127 .and.(mnum.lt.4)) then
2130 adt=d_a_old(j,inres)*d_time
2132 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
2133 d_t_new(j,inres)=d_t_old(j,inres)+adt2
2134 d_t(j,inres)=d_t_old(j,inres)+adt
2139 write (iout,*) "VELVERLET1 END: DC"
2141 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),&
2142 (dc(j,i+nres),j=1,3)
2146 end subroutine verlet1
2147 !-----------------------------------------------------------------------------
2149 ! Step 2 of the velocity Verlet algorithm: update velocities
2151 ! implicit real*8 (a-h,o-z)
2152 ! include 'DIMENSIONS'
2153 ! include 'COMMON.CONTROL'
2154 ! include 'COMMON.VAR'
2155 ! include 'COMMON.MD'
2156 ! include 'COMMON.CHAIN'
2157 ! include 'COMMON.DERIV'
2158 ! include 'COMMON.GEO'
2159 ! include 'COMMON.LOCAL'
2160 ! include 'COMMON.INTERACT'
2161 ! include 'COMMON.IOUNITS'
2162 ! include 'COMMON.NAMES'
2163 integer :: i,j,inres,mnum
2166 d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
2170 d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
2175 ! iti=iabs(itype(i,mnum))
2176 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2177 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2178 .and.(mnum.lt.4)) then
2181 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
2186 end subroutine verlet2
2187 !-----------------------------------------------------------------------------
2188 subroutine sddir_precalc
2189 ! Applying velocity Verlet algorithm - step 1 to coordinates
2190 ! implicit real*8 (a-h,o-z)
2191 ! include 'DIMENSIONS'
2197 ! include 'COMMON.CONTROL'
2198 ! include 'COMMON.VAR'
2199 ! include 'COMMON.MD'
2201 ! include 'COMMON.LANGEVIN'
2203 ! include 'COMMON.LANGEVIN.lang0'
2205 ! include 'COMMON.CHAIN'
2206 ! include 'COMMON.DERIV'
2207 ! include 'COMMON.GEO'
2208 ! include 'COMMON.LOCAL'
2209 ! include 'COMMON.INTERACT'
2210 ! include 'COMMON.IOUNITS'
2211 ! include 'COMMON.NAMES'
2212 ! include 'COMMON.TIME1'
2213 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
2214 !el common /stochcalc/ stochforcvec
2215 real(kind=8) :: time00
2217 ! Compute friction and stochastic forces
2222 time_fric=time_fric+MPI_Wtime()-time00
2224 call stochastic_force(stochforcvec)
2225 time_stoch=time_stoch+MPI_Wtime()-time00
2228 ! Compute the acceleration due to friction forces (d_af_work) and stochastic
2229 ! forces (d_as_work)
2231 ! call ginv_mult(fric_work, d_af_work)
2232 ! call ginv_mult(stochforcvec, d_as_work)
2234 call fivediaginv_mult(dimen,fric_work, d_af_work)
2235 call fivediaginv_mult(dimen,stochforcvec, d_as_work)
2237 call ginv_mult(fric_work, d_af_work)
2238 call ginv_mult(stochforcvec, d_as_work)
2242 end subroutine sddir_precalc
2243 !-----------------------------------------------------------------------------
2244 subroutine sddir_verlet1
2245 ! Applying velocity Verlet algorithm - step 1 to velocities
2248 ! implicit real*8 (a-h,o-z)
2249 ! include 'DIMENSIONS'
2250 ! include 'COMMON.CONTROL'
2251 ! include 'COMMON.VAR'
2252 ! include 'COMMON.MD'
2254 ! include 'COMMON.LANGEVIN'
2256 ! include 'COMMON.LANGEVIN.lang0'
2258 ! include 'COMMON.CHAIN'
2259 ! include 'COMMON.DERIV'
2260 ! include 'COMMON.GEO'
2261 ! include 'COMMON.LOCAL'
2262 ! include 'COMMON.INTERACT'
2263 ! include 'COMMON.IOUNITS'
2264 ! include 'COMMON.NAMES'
2265 ! Revised 3/31/05 AL: correlation between random contributions to
2266 ! position and velocity increments included.
2267 real(kind=8) :: sqrt13 = 0.57735026918962576451d0 ! 1/sqrt(3)
2268 real(kind=8) :: adt,adt2
2269 integer :: i,j,ind,inres,mnum
2271 ! Add the contribution from BOTH friction and stochastic force to the
2272 ! coordinates, but ONLY the contribution from the friction forces to velocities
2275 adt=(d_a_old(j,0)+d_af_work(j))*d_time
2276 adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
2277 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
2278 d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
2279 d_t(j,0)=d_t_old(j,0)+adt
2284 adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
2285 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
2286 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
2287 d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
2288 d_t(j,i)=d_t_old(j,i)+adt
2294 ! iti=iabs(itype(i,mnum))
2295 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2296 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2297 .and.(mnum.lt.4)) then
2300 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
2301 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
2302 ! write(iout,*) "adt",adt,"ads",adt2
2303 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
2304 d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
2305 d_t(j,inres)=d_t_old(j,inres)+adt
2311 end subroutine sddir_verlet1
2312 !-----------------------------------------------------------------------------
2313 subroutine sddir_verlet2
2314 ! Calculating the adjusted velocities for accelerations
2317 ! implicit real*8 (a-h,o-z)
2318 ! include 'DIMENSIONS'
2319 ! include 'COMMON.CONTROL'
2320 ! include 'COMMON.VAR'
2321 ! include 'COMMON.MD'
2323 ! include 'COMMON.LANGEVIN'
2325 ! include 'COMMON.LANGEVIN.lang0'
2327 ! include 'COMMON.CHAIN'
2328 ! include 'COMMON.DERIV'
2329 ! include 'COMMON.GEO'
2330 ! include 'COMMON.LOCAL'
2331 ! include 'COMMON.INTERACT'
2332 ! include 'COMMON.IOUNITS'
2333 ! include 'COMMON.NAMES'
2334 real(kind=8),dimension(6*nres) :: stochforcvec,d_as_work1 !(MAXRES6) maxres6=6*maxres
2335 real(kind=8) :: cos60 = 0.5d0, sin60 = 0.86602540378443864676d0
2336 integer :: i,j,ind,inres,mnum
2337 ! Revised 3/31/05 AL: correlation between random contributions to
2338 ! position and velocity increments included.
2339 ! The correlation coefficients are calculated at low-friction limit.
2340 ! Also, friction forces are now not calculated with new velocities.
2342 ! call friction_force
2343 call stochastic_force(stochforcvec)
2345 ! Compute the acceleration due to friction forces (d_af_work) and stochastic
2346 ! forces (d_as_work)
2349 call fivediaginv_mult(maxres6,stochforcvec, d_as_work1)
2351 call ginv_mult(stochforcvec, d_as_work1)
2358 d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j)) &
2359 +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
2364 d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j)) &
2365 +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
2371 ! iti=iabs(itype(i,mnum))
2372 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2373 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2374 .and.(mnum.lt.4)) then
2377 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres) &
2378 +d_af_work(ind+j))+sin60*d_as_work(ind+j) &
2379 +cos60*d_as_work1(ind+j))*d_time
2385 end subroutine sddir_verlet2
2386 !-----------------------------------------------------------------------------
2387 subroutine max_accel
2389 ! Find the maximum difference in the accelerations of the the sites
2390 ! at the beginning and the end of the time step.
2393 ! implicit real*8 (a-h,o-z)
2394 ! include 'DIMENSIONS'
2395 ! include 'COMMON.CONTROL'
2396 ! include 'COMMON.VAR'
2397 ! include 'COMMON.MD'
2398 ! include 'COMMON.CHAIN'
2399 ! include 'COMMON.DERIV'
2400 ! include 'COMMON.GEO'
2401 ! include 'COMMON.LOCAL'
2402 ! include 'COMMON.INTERACT'
2403 ! include 'COMMON.IOUNITS'
2404 real(kind=8),dimension(3) :: aux,accel,accel_old
2405 real(kind=8) :: dacc
2409 ! aux(j)=d_a(j,0)-d_a_old(j,0)
2410 accel_old(j)=d_a_old(j,0)
2417 ! 7/3/08 changed to asymmetric difference
2419 ! accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
2420 accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
2421 accel(j)=accel(j)+0.5d0*d_a(j,i)
2422 ! if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
2423 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
2424 dacc=dabs(accel(j)-accel_old(j))
2425 ! write (iout,*) i,dacc
2426 if (dacc.gt.amax) amax=dacc
2434 accel_old(j)=d_a_old(j,0)
2439 accel_old(j)=accel_old(j)+d_a_old(j,1)
2440 accel(j)=accel(j)+d_a(j,1)
2445 ! iti=iabs(itype(i,mnum))
2446 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2447 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2448 .and.(mnum.lt.4)) then
2450 ! accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
2451 accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
2452 accel(j)=accel(j)+d_a(j,i+nres)
2456 ! if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
2457 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
2458 dacc=dabs(accel(j)-accel_old(j))
2459 ! write (iout,*) "side-chain",i,dacc
2460 if (dacc.gt.amax) amax=dacc
2464 accel_old(j)=accel_old(j)+d_a_old(j,i)
2465 accel(j)=accel(j)+d_a(j,i)
2466 ! aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
2470 end subroutine max_accel
2471 !-----------------------------------------------------------------------------
2472 subroutine predict_edrift(epdrift)
2474 ! Predict the drift of the potential energy
2477 use control_data, only: lmuca
2478 ! implicit real*8 (a-h,o-z)
2479 ! include 'DIMENSIONS'
2480 ! include 'COMMON.CONTROL'
2481 ! include 'COMMON.VAR'
2482 ! include 'COMMON.MD'
2483 ! include 'COMMON.CHAIN'
2484 ! include 'COMMON.DERIV'
2485 ! include 'COMMON.GEO'
2486 ! include 'COMMON.LOCAL'
2487 ! include 'COMMON.INTERACT'
2488 ! include 'COMMON.IOUNITS'
2489 ! include 'COMMON.MUCA'
2490 real(kind=8) :: epdrift,epdriftij
2492 ! Drift of the potential energy
2498 epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
2499 if (lmuca) epdriftij=epdriftij*factor
2500 ! write (iout,*) "back",i,j,epdriftij
2501 if (epdriftij.gt.epdrift) epdrift=epdriftij
2505 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1.and.&
2506 molnum(i).lt.4) then
2509 dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
2510 if (lmuca) epdriftij=epdriftij*factor
2511 ! write (iout,*) "side",i,j,epdriftij
2512 if (epdriftij.gt.epdrift) epdrift=epdriftij
2516 epdrift=0.5d0*epdrift*d_time*d_time
2517 ! write (iout,*) "epdrift",epdrift
2519 end subroutine predict_edrift
2520 !-----------------------------------------------------------------------------
2521 subroutine verlet_bath
2523 ! Coupling to the thermostat by using the Berendsen algorithm
2526 ! implicit real*8 (a-h,o-z)
2527 ! include 'DIMENSIONS'
2528 ! include 'COMMON.CONTROL'
2529 ! include 'COMMON.VAR'
2530 ! include 'COMMON.MD'
2531 ! include 'COMMON.CHAIN'
2532 ! include 'COMMON.DERIV'
2533 ! include 'COMMON.GEO'
2534 ! include 'COMMON.LOCAL'
2535 ! include 'COMMON.INTERACT'
2536 ! include 'COMMON.IOUNITS'
2537 ! include 'COMMON.NAMES'
2538 real(kind=8) :: T_half,fact
2539 integer :: i,j,inres,mnum
2541 T_half=2.0d0/(dimen3*Rb)*EK
2542 fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
2543 ! write(iout,*) "T_half", T_half
2544 ! write(iout,*) "EK", EK
2545 ! write(iout,*) "fact", fact
2547 d_t(j,0)=fact*d_t(j,0)
2551 d_t(j,i)=fact*d_t(j,i)
2556 ! iti=iabs(itype(i,mnum))
2557 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2558 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2559 .and.(mnum.lt.4)) then
2562 d_t(j,inres)=fact*d_t(j,inres)
2567 end subroutine verlet_bath
2568 !-----------------------------------------------------------------------------
2570 ! Set up the initial conditions of a MD simulation
2573 use control, only:tcpu
2574 !el use io_basic, only:ilen
2577 use minimm, only:minim_dc,minimize,sc_move
2578 use io_config, only:readrst
2579 use io, only:statout
2580 use random, only: iran_num
2581 ! implicit real*8 (a-h,o-z)
2582 ! include 'DIMENSIONS'
2585 character(len=16) :: form
2586 integer :: IERROR,ERRCODE
2588 integer :: iranmin,itrial,itmp,n_model_try,k, &
2590 integer, dimension(:),allocatable :: list_model_try
2591 integer, dimension(0:nodes-1) :: i_start_models
2592 ! include 'COMMON.SETUP'
2593 ! include 'COMMON.CONTROL'
2594 ! include 'COMMON.VAR'
2595 ! include 'COMMON.MD'
2597 ! include 'COMMON.LANGEVIN'
2599 ! include 'COMMON.LANGEVIN.lang0'
2601 ! include 'COMMON.CHAIN'
2602 ! include 'COMMON.DERIV'
2603 ! include 'COMMON.GEO'
2604 ! include 'COMMON.LOCAL'
2605 ! include 'COMMON.INTERACT'
2606 ! include 'COMMON.IOUNITS'
2607 ! include 'COMMON.NAMES'
2608 ! include 'COMMON.REMD'
2609 real(kind=8),dimension(0:n_ene) :: energia_long,energia_short,energia
2610 real(kind=8),dimension(3) :: vcm,incr,L
2611 real(kind=8) :: xv,sigv,lowb,highb
2612 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
2613 character(len=256) :: qstr
2616 character(len=50) :: tytul
2617 logical :: file_exist
2618 !el common /gucio/ cm
2619 integer :: i,j,ipos,iq,iw,nft_sc,iretcode,nfun,ierr,mnum,itime
2620 real(kind=8) :: etot,tt0
2624 ! write(iout,*) "d_time", d_time
2625 ! Compute the standard deviations of stochastic forces for Langevin dynamics
2626 ! if the friction coefficients do not depend on surface area
2627 if (lang.gt.0 .and. .not.surfarea) then
2630 stdforcp(i)=stdfp(mnum)*dsqrt(gamp(mnum))
2634 stdforcsc(i)=stdfsc(iabs(itype(i,mnum)),mnum) &
2635 *dsqrt(gamsc(iabs(itype(i,mnum)),mnum))
2639 ! Open the pdb file for snapshotshots
2642 if (ilen(tmpdir).gt.0) &
2643 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2644 liczba(:ilen(liczba))//".pdb")
2646 file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2650 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) &
2651 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2652 liczba(:ilen(liczba))//".x")
2653 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2656 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) &
2657 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2658 liczba(:ilen(liczba))//".cx")
2659 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2665 if (ilen(tmpdir).gt.0) &
2666 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
2667 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
2669 if (ilen(tmpdir).gt.0) &
2670 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
2671 cartname=prefix(:ilen(prefix))//"_MD.cx"
2675 write (qstr,'(256(1h ))')
2678 iq = qinfrag(i,iset)*10
2679 iw = wfrag(i,iset)/100
2681 if(me.eq.king.or..not.out1file) &
2682 write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
2683 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
2688 iq = qinpair(i,iset)*10
2689 iw = wpair(i,iset)/100
2691 if(me.eq.king.or..not.out1file) &
2692 write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
2693 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
2697 ! pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
2699 ! cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
2701 ! cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
2703 ! statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
2707 if (restart1file) then
2709 inquire(file=mremd_rst_name,exist=file_exist)
2710 write (*,*) me," Before broadcast: file_exist",file_exist
2712 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,&
2715 write (*,*) me," After broadcast: file_exist",file_exist
2716 ! inquire(file=mremd_rst_name,exist=file_exist)
2717 if(me.eq.king.or..not.out1file) &
2718 write(iout,*) "Initial state read by master and distributed"
2720 if (ilen(tmpdir).gt.0) &
2721 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_' &
2722 //liczba(:ilen(liczba))//'.rst')
2723 inquire(file=rest2name,exist=file_exist)
2726 if(.not.restart1file) then
2727 if(me.eq.king.or..not.out1file) &
2728 write(iout,*) "Initial state will be read from file ",&
2729 rest2name(:ilen(rest2name))
2732 call rescale_weights(t_bath)
2734 if(me.eq.king.or..not.out1file)then
2735 if (restart1file) then
2736 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),&
2739 write(iout,*) "File ",rest2name(:ilen(rest2name)),&
2742 write(iout,*) "Initial velocities randomly generated"
2749 ! Generate initial velocities
2750 if(me.eq.king.or..not.out1file) &
2751 write(iout,*) "Initial velocities randomly generated"
2756 ! rest2name = prefix(:ilen(prefix))//'.rst'
2757 if(me.eq.king.or..not.out1file)then
2758 write (iout,*) "Initial velocities"
2760 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2761 (d_t(j,i+nres),j=1,3)
2763 ! Zeroing the total angular momentum of the system
2764 write(iout,*) "Calling the zero-angular momentum subroutine"
2767 ! Getting the potential energy and forces and velocities and accelerations
2769 ! write (iout,*) "velocity of the center of the mass:"
2770 ! write (iout,*) (vcm(j),j=1,3)
2772 d_t(j,0)=d_t(j,0)-vcm(j)
2774 ! Removing the velocity of the center of mass
2776 if(me.eq.king.or..not.out1file)then
2777 write (iout,*) "vcm right after adjustment:"
2778 write (iout,*) (vcm(j),j=1,3)
2785 if ((.not.rest).or.(forceminim)) then
2786 if (forceminim) call chainbuild_cart
2788 if(iranconf.ne.0 .or.indpdb.gt.0.and..not.unres_pdb .or.preminim) then
2790 print *, 'Calling OVERLAP_SC'
2791 call overlap_sc(fail)
2792 print *,'after OVERLAP'
2795 print *,'call SC_MOVE'
2796 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
2797 print *,'SC_move',nft_sc,etot
2798 if(me.eq.king.or..not.out1file) &
2799 write(iout,*) 'SC_move',nft_sc,etot
2803 print *, 'Calling MINIM_DC'
2804 call minim_dc(etot,iretcode,nfun)
2806 call geom_to_var(nvar,varia)
2807 print *,'Calling MINIMIZE.'
2808 call minimize(etot,varia,iretcode,nfun)
2809 call var_to_geom(nvar,varia)
2811 write(iout,*) "just before minimin"
2813 if(me.eq.king.or..not.out1file) &
2814 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2816 write(iout,*) "just after minimin"
2818 if(iranconf.ne.0) then
2819 !c 8/22/17 AL Loop to produce a low-energy random conformation
2822 if(me.eq.king.or..not.out1file) &
2823 write (iout,*) 'Calling OVERLAP_SC'
2824 call overlap_sc(fail)
2825 endif !endif overlap
2828 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
2829 print *,'SC_move',nft_sc,etot
2830 if(me.eq.king.or..not.out1file) &
2831 write(iout,*) 'SC_move',nft_sc,etot
2835 print *, 'Calling MINIM_DC'
2836 call minim_dc(etot,iretcode,nfun)
2837 call int_from_cart1(.false.)
2839 call geom_to_var(nvar,varia)
2840 print *,'Calling MINIMIZE.'
2841 call minimize(etot,varia,iretcode,nfun)
2842 call var_to_geom(nvar,varia)
2844 if(me.eq.king.or..not.out1file) &
2845 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2846 write(iout,*) "just after minimin"
2848 if (isnan(etot) .or. etot.gt.4.0d6) then
2849 write (iout,*) "Energy too large",etot, &
2850 " trying another random conformation"
2853 call gen_rand_conf(itmp,*30)
2855 30 write (iout,*) 'Failed to generate random conformation', &
2857 write (*,*) 'Processor:',me, &
2858 ' Failed to generate random conformation',&
2867 write (iout,'(a,i3,a)') 'Processor:',me, &
2868 ' error in generating random conformation.'
2869 write (*,'(a,i3,a)') 'Processor:',me, &
2870 ' error in generating random conformation.'
2873 ! call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
2874 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
2884 write (iout,'(a,i3,a)') 'Processor:',me, &
2885 ' failed to generate a low-energy random conformation.'
2886 write (*,'(a,i3,a,f10.3)') 'Processor:',me, &
2887 ' failed to generate a low-energy random conformation.',etot
2891 ! call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
2892 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
2897 else if (preminim) then
2898 if (start_from_model) then
2902 do while (fail .and. n_model_try.lt.nmodel_start)
2903 write (iout,*) "n_model_try",n_model_try
2905 i_model=iran_num(1,nmodel_start)
2907 if (i_model.eq.list_model_try(k)) exit
2909 if (k.gt.n_model_try) exit
2911 n_model_try=n_model_try+1
2912 list_model_try(n_model_try)=i_model
2913 if (me.eq.king .or. .not. out1file) &
2914 write (iout,*) 'Trying to start from model ',&
2915 pdbfiles_chomo(i_model)(:ilen(pdbfiles_chomo(i_model)))
2918 c(j,i)=chomo(j,i,i_model)
2921 call int_from_cart(.true.,.false.)
2922 call sc_loc_geom(.false.)
2926 dc(j,i)=c(j,i+1)-c(j,i)
2927 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
2932 dc(j,i+nres)=c(j,i+nres)-c(j,i)
2933 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
2936 if (me.eq.king.or..not.out1file) then
2937 write (iout,*) "Energies before removing overlaps"
2938 call etotal(energia(0))
2939 call enerprint(energia(0))
2941 ! Remove SC overlaps if requested
2943 write (iout,*) 'Calling OVERLAP_SC'
2944 call overlap_sc(fail)
2947 "Failed to remove overlap from model",i_model
2951 if (me.eq.king.or..not.out1file) then
2952 write (iout,*) "Energies after removing overlaps"
2953 call etotal(energia(0))
2954 call enerprint(energia(0))
2957 ! Search for better SC rotamers if requested
2959 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
2960 print *,'SC_move',nft_sc,etot
2961 if (me.eq.king.or..not.out1file)&
2962 write(iout,*) 'SC_move',nft_sc,etot
2964 call etotal(energia(0))
2967 call MPI_Gather(i_model,1,MPI_INTEGER,i_start_models(0),&
2968 1,MPI_INTEGER,king,CG_COMM,IERROR)
2969 if (n_model_try.gt.nmodel_start .and.&
2970 (me.eq.king .or. out1file)) then
2972 "All models have irreparable overlaps. Trying randoms starts."
2974 i_model=nmodel_start+1
2978 ! Remove SC overlaps if requested
2980 write (iout,*) 'Calling OVERLAP_SC'
2981 call overlap_sc(fail)
2984 "Failed to remove overlap"
2987 if (me.eq.king.or..not.out1file) then
2988 write (iout,*) "Energies after removing overlaps"
2989 call etotal(energia(0))
2990 call enerprint(energia(0))
2993 ! 8/22/17 AL Minimize initial structure
2995 if (me.eq.king.or..not.out1file) write(iout,*)&
2996 'Minimizing initial PDB structure: Calling MINIM_DC'
2997 call minim_dc(etot,iretcode,nfun)
2999 call geom_to_var(nvar,varia)
3000 if(me.eq.king.or..not.out1file) write (iout,*)&
3001 'Minimizing initial PDB structure: Calling MINIMIZE.'
3002 call minimize(etot,varia,iretcode,nfun)
3003 call var_to_geom(nvar,varia)
3005 if (me.eq.king.or..not.out1file)&
3006 write(iout,*) 'LBFGS return code is ',status,' eval ',nfun
3007 if(me.eq.king.or..not.out1file)&
3008 write(iout,*) 'LBFGS return code is ',status,' eval ',nfun
3010 if (me.eq.king.or..not.out1file)&
3011 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
3012 if(me.eq.king.or..not.out1file)&
3013 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
3017 if (nmodel_start.gt.0 .and. me.eq.king) then
3018 write (iout,'(a)') "Task Starting model"
3020 if (i_start_models(i).gt.nmodel_start) then
3021 write (iout,'(i4,2x,a)') i,"RANDOM STRUCTURE"
3023 write(iout,'(i4,2x,a)')i,pdbfiles_chomo(i_start_models(i)) &
3024 (:ilen(pdbfiles_chomo(i_start_models(i))))
3029 call chainbuild_cart
3031 write(iout,*) "just after kinetic"
3036 kinetic_T=2.0d0/(dimen3*Rb)*EK
3037 if(me.eq.king.or..not.out1file)then
3038 write(iout,*) "just after verlet_bath"
3048 write(iout,*) "before ETOTAL"
3049 call etotal(potEcomp)
3050 if (large) call enerprint(potEcomp)
3053 t_etotal=t_etotal+MPI_Wtime()-tt0
3055 t_etotal=t_etotal+tcpu()-tt0
3062 if (amax*d_time .gt. dvmax) then
3063 d_time=d_time*dvmax/amax
3064 if(me.eq.king.or..not.out1file) write (iout,*) &
3065 "Time step reduced to",d_time,&
3066 " because of too large initial acceleration."
3068 if(me.eq.king.or..not.out1file)then
3069 write(iout,*) "Potential energy and its components"
3070 call enerprint(potEcomp)
3071 ! write(iout,*) (potEcomp(i),i=0,n_ene)
3073 potE=potEcomp(0)-potEcomp(51)
3077 if (ntwe.ne.0) call statout(itime)
3078 if(me.eq.king.or..not.out1file) &
3079 write (iout,'(/a/3(a25,1pe14.5/))') "Initial:", &
3080 " Kinetic energy",EK," Potential energy",potE, &
3081 " Total energy",totE," Maximum acceleration ", &
3084 write (iout,*) "Initial coordinates"
3086 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),&
3089 write (iout,*) "Initial dC"
3091 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),&
3092 (dc(j,i+nres),j=1,3)
3094 write (iout,*) "Initial velocities"
3095 write (iout,"(13x,' backbone ',23x,' side chain')")
3097 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
3098 (d_t(j,i+nres),j=1,3)
3100 write (iout,*) "Initial accelerations"
3102 ! write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
3103 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),&
3104 (d_a(j,i+nres),j=1,3)
3110 d_t_old(j,i)=d_t(j,i)
3111 d_a_old(j,i)=d_a(j,i)
3113 ! write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
3122 call etotal_short(energia_short)
3123 if (large) call enerprint(potEcomp)
3126 t_eshort=t_eshort+MPI_Wtime()-tt0
3128 t_eshort=t_eshort+tcpu()-tt0
3133 if(.not.out1file .and. large) then
3134 write (iout,*) "energia_long",energia_long(0),&
3135 " energia_short",energia_short(0),&
3136 " total",energia_long(0)+energia_short(0)
3137 write (iout,*) "Initial fast-force accelerations"
3139 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
3140 (d_a(j,i+nres),j=1,3)
3143 ! 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
3146 d_a_short(j,i)=d_a(j,i)
3155 call etotal_long(energia_long)
3156 if (large) call enerprint(potEcomp)
3159 t_elong=t_elong+MPI_Wtime()-tt0
3161 t_elong=t_elong+tcpu()-tt0
3166 if(.not.out1file .and. large) then
3167 write (iout,*) "energia_long",energia_long(0)
3168 write (iout,*) "Initial slow-force accelerations"
3170 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
3171 (d_a(j,i+nres),j=1,3)
3175 t_enegrad=t_enegrad+MPI_Wtime()-tt0
3177 t_enegrad=t_enegrad+tcpu()-tt0
3181 end subroutine init_MD
3182 !-----------------------------------------------------------------------------
3183 subroutine random_vel
3185 ! implicit real*8 (a-h,o-z)
3187 use random, only:anorm_distr
3189 ! include 'DIMENSIONS'
3190 ! include 'COMMON.CONTROL'
3191 ! include 'COMMON.VAR'
3192 ! include 'COMMON.MD'
3194 ! include 'COMMON.LANGEVIN'
3196 ! include 'COMMON.LANGEVIN.lang0'
3198 ! include 'COMMON.CHAIN'
3199 ! include 'COMMON.DERIV'
3200 ! include 'COMMON.GEO'
3201 ! include 'COMMON.LOCAL'
3202 ! include 'COMMON.INTERACT'
3203 ! include 'COMMON.IOUNITS'
3204 ! include 'COMMON.NAMES'
3205 ! include 'COMMON.TIME1'
3206 real(kind=8) :: xv,sigv,lowb,highb ,Ek1
3208 integer ichain,n,innt,inct,ibeg,ierr
3209 double precision work(48*nres)
3210 integer iwork(maxres6)
3211 double precision Ghalf(mmaxres2_chain),Geigen(maxres2_chain),
3212 & Gvec(maxres2_chain,maxres2_chain)
3213 common /przechowalnia/Ghalf,Geigen,Gvec
3215 double precision inertia(maxres2_chain,maxres2_chain)
3220 real(kind=8) ,allocatable, dimension(:) :: DDU1,DDU2,DL2,DL1,xsolv,DML,rs
3221 real(kind=8) :: sumx
3223 real(kind=8) ,allocatable, dimension(:) :: rsold
3224 real (kind=8),allocatable,dimension(:,:) :: matold
3228 integer :: i,j,ii,k,ind,mark,imark,mnum
3229 ! Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
3230 ! First generate velocities in the eigenspace of the G matrix
3231 ! write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
3238 sigv=dsqrt((Rb*t_bath)/geigen(i))
3241 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
3243 write (iout,*) "i",i," ii",ii," geigen",geigen(i),&
3244 " d_t_work_new",d_t_work_new(ii)
3255 Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
3258 write (iout,*) "Ek from eigenvectors",Ek1
3259 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
3263 ! Transform velocities to UNRES coordinate space
3264 allocate (DL1(2*nres))
3265 allocate (DDU1(2*nres))
3266 allocate (DL2(2*nres))
3267 allocate (DDU2(2*nres))
3268 allocate (xsolv(2*nres))
3269 allocate (DML(2*nres))
3270 allocate (rs(2*nres))
3272 allocate (rsold(2*nres))
3273 allocate (matold(2*nres,2*nres))
3275 matold(1,1)=DMorig(1)
3276 matold(1,2)=DU1orig(1)
3277 matold(1,3)=DU2orig(1)
3278 write (*,*) DMorig(1),DU1orig(1),DU2orig(1)
3283 matold(i,j)=DMorig(i)
3284 matold(i,j-1)=DU1orig(i-1)
3286 matold(i,j-2)=DU2orig(i-2)
3294 matold(i,j+1)=DU1orig(i)
3300 matold(i,j+2)=DU2orig(i)
3304 matold(dimen,dimen)=DMorig(dimen)
3305 matold(dimen,dimen-1)=DU1orig(dimen-1)
3306 matold(dimen,dimen-2)=DU2orig(dimen-2)
3307 write(iout,*) "old gmatrix"
3308 call matout(dimen,dimen,2*nres,2*nres,matold)
3312 ! Find the ith eigenvector of the pentadiagonal inertiq matrix
3316 DML(j)=DMorig(j)-geigen(i)
3319 DML(j-1)=DMorig(j)-geigen(i)
3324 DDU1(imark-1)=DU2orig(imark-1)
3325 do j=imark+1,dimen-1
3326 DDU1(j-1)=DU1orig(j)
3334 DDU2(j)=DU2orig(j+1)
3343 write (iout,*) "DL2,DL1,DML,DDU1,DDU2"
3344 write(iout,'(10f10.5)') (DL2(k),k=3,dimen-1)
3345 write(iout,'(10f10.5)') (DL1(k),k=2,dimen-1)
3346 write(iout,'(10f10.5)') (DML(k),k=1,dimen-1)
3347 write(iout,'(10f10.5)') (DDU1(k),k=1,dimen-2)
3348 write(iout,'(10f10.5)') (DDU2(k),k=1,dimen-3)
3351 if (imark.gt.2) rs(imark-2)=-DU2orig(imark-2)
3352 if (imark.gt.1) rs(imark-1)=-DU1orig(imark-1)
3353 if (imark.lt.dimen) rs(imark)=-DU1orig(imark)
3354 if (imark.lt.dimen-1) rs(imark+1)=-DU2orig(imark)
3358 ! write (iout,*) "Vector rs"
3360 ! write (iout,*) j,rs(j)
3363 call FDIAG(dimen-1,DL2,DL1,DML,DDU1,DDU2,rs,xsolv,mark)
3370 sumx=-geigen(i)*xsolv(j)
3372 sumx=sumx+matold(j,k)*xsolv(k)
3375 sumx=sumx+matold(j,k)*xsolv(k-1)
3377 write(iout,'(i5,3f10.5)') j,sumx,rsold(j),sumx-rsold(j)
3380 sumx=-geigen(i)*xsolv(j-1)
3382 sumx=sumx+matold(j,k)*xsolv(k)
3385 sumx=sumx+matold(j,k)*xsolv(k-1)
3387 write(iout,'(i5,3f10.5)') j-1,sumx,rsold(j-1),sumx-rsold(j-1)
3391 "Solution of equations system",i," for eigenvalue",geigen(i)
3393 write(iout,'(i5,f10.5)') j,xsolv(j)
3396 do j=dimen-1,imark,-1
3401 write (iout,*) "Un-normalized eigenvector",i," for eigenvalue",geigen(i)
3403 write(iout,'(i5,f10.5)') j,xsolv(j)
3406 ! Normalize ith eigenvector
3409 sumx=sumx+xsolv(j)**2
3413 xsolv(j)=xsolv(j)/sumx
3416 write (iout,*) "Eigenvector",i," for eigenvalue",geigen(i)
3418 write(iout,'(i5,3f10.5)') j,xsolv(j)
3421 ! All done at this point for eigenvector i, exit loop
3429 write (iout,*) "Unable to find eigenvector",i
3432 ! write (iout,*) "i=",i
3434 ! write (iout,*) "k=",k
3437 ! write(iout,*) "ind",ind," ind1",3*(i-1)+k
3438 d_t_work(ind)=d_t_work(ind) &
3439 +xsolv(j)*d_t_work_new(3*(i-1)+k)
3442 enddo ! i (loop over the eigenvectors)
3445 write (iout,*) "d_t_work"
3447 write (iout,"(i5,f10.5)") i,d_t_work(i)
3452 ! if (itype(i,1).eq.10) then
3454 if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
3455 iti=iabs(itype(i,mnum))
3456 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3457 if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
3458 .or.(mnum.ge.4)) then
3465 ! write (iout,*) "k",k," ii+k",ii+k," ii+j+k",ii+j+k,"EK1",Ek1
3466 Ek1=Ek1+0.5d0*mp(mnum)*((d_t_work(ii+j+k)+d_t_work_new(ii+k))/2)**2+&
3467 0.5d0*0.25d0*IP(mnum)*(d_t_work(ii+j+k)-d_t_work_new(ii+k))**2
3470 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3471 .and.(mnum.lt.4)) ii=ii+3
3472 write (iout,*) "i",i," itype",itype(i,mnum)," mass",msc(itype(i,mnum),mnum)
3473 write (iout,*) "ii",ii
3476 write (iout,*) "k",k," ii",ii,"EK1",EK1
3477 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3479 Ek1=Ek1+0.5d0*Isc(iabs(itype(i,mnum)),mnum)*(d_t_work(ii)-d_t_work(ii-3))**2
3480 Ek1=Ek1+0.5d0*msc(iabs(itype(i,mnum)),mnum)*d_t_work(ii)**2
3482 write (iout,*) "i",i," ii",ii
3484 write (iout,*) "Ek from d_t_work",Ek1
3485 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
3487 if(allocated(DDU1)) deallocate(DDU1)
3488 if(allocated(DDU2)) deallocate(DDU2)
3489 if(allocated(DL2)) deallocate(DL2)
3490 if(allocated(DL1)) deallocate(DL1)
3491 if(allocated(xsolv)) deallocate(xsolv)
3492 if(allocated(DML)) deallocate(DML)
3493 if(allocated(rs)) deallocate(rs)
3495 if(allocated(matold)) deallocate(matold)
3496 if(allocated(rsold)) deallocate(rsold)
3501 d_t(k,j)=d_t_work(ind)
3505 if (itype(j,1).ne.10 .and. itype(j,mnum).ne.ntyp1_molec(mnum)&
3506 .and.(mnum.lt.4)) then
3508 d_t(k,j+nres)=d_t_work(ind)
3514 write (iout,*) "Random velocities in the Calpha,SC space"
3516 write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3)
3519 write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3)
3526 ! if (itype(i,1).eq.10) then
3528 if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
3529 .or.(mnum.ge.4)) then
3531 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
3535 d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
3536 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
3541 write (iout,*)"Random velocities in the virtual-bond-vector space"
3543 write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3)
3546 write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3)
3549 write (iout,*) "Ek from d_t_work",Ek1
3550 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
3558 d_t_work(ind)=d_t_work(ind) &
3559 +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
3561 ! write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
3565 ! Transfer to the d_t vector
3567 d_t(j,0)=d_t_work(j)
3573 d_t(j,i)=d_t_work(ind)
3578 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3579 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3580 .and.(mnum.lt.4)) then
3583 d_t(j,i+nres)=d_t_work(ind)
3589 ! write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",&
3590 ! 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
3592 ! write(iout,*) "end init MD"
3594 end subroutine random_vel
3595 !-----------------------------------------------------------------------------
3597 subroutine sd_verlet_p_setup
3598 ! Sets up the parameters of stochastic Verlet algorithm
3599 ! implicit real*8 (a-h,o-z)
3600 ! include 'DIMENSIONS'
3601 use control, only: tcpu
3606 ! include 'COMMON.CONTROL'
3607 ! include 'COMMON.VAR'
3608 ! include 'COMMON.MD'
3610 ! include 'COMMON.LANGEVIN'
3612 ! include 'COMMON.LANGEVIN.lang0'
3614 ! include 'COMMON.CHAIN'
3615 ! include 'COMMON.DERIV'
3616 ! include 'COMMON.GEO'
3617 ! include 'COMMON.LOCAL'
3618 ! include 'COMMON.INTERACT'
3619 ! include 'COMMON.IOUNITS'
3620 ! include 'COMMON.NAMES'
3621 ! include 'COMMON.TIME1'
3622 real(kind=8),dimension(6*nres) :: emgdt !(MAXRES6) maxres6=6*maxres
3623 real(kind=8) :: pterm,vterm,rho,rhoc,vsig
3624 real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
3625 prand_vec,vrand_vec1,vrand_vec2 !(MAXRES6) maxres6=6*maxres
3626 logical :: lprn = .false.
3627 real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
3628 real(kind=8) :: ktm,gdt,egdt,gdt2,gdt3,gdt4,gdt5,gdt6,gdt7,gdt8,&
3630 integer :: i,maxres2
3637 ! AL 8/17/04 Code adapted from tinker
3639 ! Get the frictional and random terms for stochastic dynamics in the
3640 ! eigenspace of mass-scaled UNRES friction matrix
3644 gdt = fricgam(i) * d_time
3646 ! Stochastic dynamics reduces to simple MD for zero friction
3648 if (gdt .le. zero) then
3649 pfric_vec(i) = 1.0d0
3650 vfric_vec(i) = d_time
3651 afric_vec(i) = 0.5d0 * d_time * d_time
3652 prand_vec(i) = 0.0d0
3653 vrand_vec1(i) = 0.0d0
3654 vrand_vec2(i) = 0.0d0
3656 ! Analytical expressions when friction coefficient is large
3659 if (gdt .ge. gdt_radius) then
3662 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
3663 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
3664 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
3665 vterm = 1.0d0 - egdt**2
3666 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
3668 ! Use series expansions when friction coefficient is small
3679 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0 &
3680 - gdt5/120.0d0 + gdt6/720.0d0 &
3681 - gdt7/5040.0d0 + gdt8/40320.0d0 &
3682 - gdt9/362880.0d0) / fricgam(i)**2
3683 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
3684 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
3685 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0 &
3686 + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0 &
3687 + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0 &
3688 + 127.0d0*gdt9/90720.0d0
3689 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0 &
3690 - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0 &
3691 - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0 &
3692 - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
3693 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0 &
3694 - 17.0d0*gdt2/1280.0d0 &
3695 + 17.0d0*gdt3/6144.0d0 &
3696 + 40967.0d0*gdt4/34406400.0d0 &
3697 - 57203.0d0*gdt5/275251200.0d0 &
3698 - 1429487.0d0*gdt6/13212057600.0d0)
3701 ! Compute the scaling factors of random terms for the nonzero friction case
3703 ktm = 0.5d0*d_time/fricgam(i)
3704 psig = dsqrt(ktm*pterm) / fricgam(i)
3705 vsig = dsqrt(ktm*vterm)
3706 rhoc = dsqrt(1.0d0 - rho*rho)
3708 vrand_vec1(i) = vsig * rho
3709 vrand_vec2(i) = vsig * rhoc
3714 "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",&
3717 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),&
3718 afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
3722 ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
3725 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
3726 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
3727 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
3728 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
3729 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
3730 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
3733 t_sdsetup=t_sdsetup+MPI_Wtime()
3735 t_sdsetup=t_sdsetup+tcpu()-tt0
3738 end subroutine sd_verlet_p_setup
3739 !-----------------------------------------------------------------------------
3740 subroutine eigtransf1(n,ndim,ab,d,c)
3744 real(kind=8) :: ab(ndim,ndim,n),c(ndim,n),d(ndim)
3750 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
3755 end subroutine eigtransf1
3756 !-----------------------------------------------------------------------------
3757 subroutine eigtransf(n,ndim,a,b,d,c)
3761 real(kind=8) :: a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
3767 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
3772 end subroutine eigtransf
3773 !-----------------------------------------------------------------------------
3774 subroutine sd_verlet1
3776 ! Applying stochastic velocity Verlet algorithm - step 1 to velocities
3778 ! implicit real*8 (a-h,o-z)
3779 ! include 'DIMENSIONS'
3780 ! include 'COMMON.CONTROL'
3781 ! include 'COMMON.VAR'
3782 ! include 'COMMON.MD'
3784 ! include 'COMMON.LANGEVIN'
3786 ! include 'COMMON.LANGEVIN.lang0'
3788 ! include 'COMMON.CHAIN'
3789 ! include 'COMMON.DERIV'
3790 ! include 'COMMON.GEO'
3791 ! include 'COMMON.LOCAL'
3792 ! include 'COMMON.INTERACT'
3793 ! include 'COMMON.IOUNITS'
3794 ! include 'COMMON.NAMES'
3795 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
3796 !el common /stochcalc/ stochforcvec
3797 logical :: lprn = .false.
3798 real(kind=8) :: ddt1,ddt2
3799 integer :: i,j,ind,inres
3801 ! write (iout,*) "dc_old"
3803 ! write (iout,'(i5,3f10.5,5x,3f10.5)')
3804 ! & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3807 dc_work(j)=dc_old(j,0)
3808 d_t_work(j)=d_t_old(j,0)
3809 d_a_work(j)=d_a_old(j,0)
3814 dc_work(ind+j)=dc_old(j,i)
3815 d_t_work(ind+j)=d_t_old(j,i)
3816 d_a_work(ind+j)=d_a_old(j,i)
3822 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3823 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3824 .and.(mnum.lt.4)) then
3826 dc_work(ind+j)=dc_old(j,i+nres)
3827 d_t_work(ind+j)=d_t_old(j,i+nres)
3828 d_a_work(ind+j)=d_a_old(j,i+nres)
3836 "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",&
3840 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
3841 vfric_mat(i,j),afric_mat(i,j),&
3842 prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3850 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) &
3851 +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3852 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3853 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3855 d_t_work_new(i)=ddt1+0.5d0*ddt2
3856 d_t_work(i)=ddt1+ddt2
3861 d_t(j,0)=d_t_work(j)
3866 dc(j,i)=dc_work(ind+j)
3867 d_t(j,i)=d_t_work(ind+j)
3873 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3874 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3875 .and.(mnum.lt.4)) then
3878 dc(j,inres)=dc_work(ind+j)
3879 d_t(j,inres)=d_t_work(ind+j)
3885 end subroutine sd_verlet1
3886 !-----------------------------------------------------------------------------
3887 subroutine sd_verlet2
3889 ! Calculating the adjusted velocities for accelerations
3891 ! implicit real*8 (a-h,o-z)
3892 ! include 'DIMENSIONS'
3893 ! include 'COMMON.CONTROL'
3894 ! include 'COMMON.VAR'
3895 ! include 'COMMON.MD'
3897 ! include 'COMMON.LANGEVIN'
3899 ! include 'COMMON.LANGEVIN.lang0'
3901 ! include 'COMMON.CHAIN'
3902 ! include 'COMMON.DERIV'
3903 ! include 'COMMON.GEO'
3904 ! include 'COMMON.LOCAL'
3905 ! include 'COMMON.INTERACT'
3906 ! include 'COMMON.IOUNITS'
3907 ! include 'COMMON.NAMES'
3908 !el real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV !(MAXRES6) maxres6=6*maxres
3909 real(kind=8),dimension(6*nres) :: stochforcvecV !(MAXRES6) maxres6=6*maxres
3910 !el common /stochcalc/ stochforcvec
3912 real(kind=8) :: ddt1,ddt2
3913 integer :: i,j,ind,inres
3914 ! Compute the stochastic forces which contribute to velocity change
3916 call stochastic_force(stochforcvecV)
3923 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3924 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+ &
3925 vrand_mat2(i,j)*stochforcvecV(j)
3927 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
3931 d_t(j,0)=d_t_work(j)
3936 d_t(j,i)=d_t_work(ind+j)
3942 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3943 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3944 .and.(mnum.lt.4)) then
3947 d_t(j,inres)=d_t_work(ind+j)
3953 end subroutine sd_verlet2
3954 !-----------------------------------------------------------------------------
3955 subroutine sd_verlet_ciccotti_setup
3957 ! Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
3959 ! implicit real*8 (a-h,o-z)
3960 ! include 'DIMENSIONS'
3961 use control, only: tcpu
3966 ! include 'COMMON.CONTROL'
3967 ! include 'COMMON.VAR'
3968 ! include 'COMMON.MD'
3970 ! include 'COMMON.LANGEVIN'
3972 ! include 'COMMON.LANGEVIN.lang0'
3974 ! include 'COMMON.CHAIN'
3975 ! include 'COMMON.DERIV'
3976 ! include 'COMMON.GEO'
3977 ! include 'COMMON.LOCAL'
3978 ! include 'COMMON.INTERACT'
3979 ! include 'COMMON.IOUNITS'
3980 ! include 'COMMON.NAMES'
3981 ! include 'COMMON.TIME1'
3982 real(kind=8),dimension(6*nres) :: emgdt !(MAXRES6) maxres6=6*maxres
3983 real(kind=8) :: pterm,vterm,rho,rhoc,vsig
3984 real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
3985 prand_vec,vrand_vec1,vrand_vec2 !(MAXRES6) maxres6=6*maxres
3986 logical :: lprn = .false.
3987 real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
3988 real(kind=8) :: ktm,gdt,egdt,tt0
3989 integer :: i,maxres2
3996 ! AL 8/17/04 Code adapted from tinker
3998 ! Get the frictional and random terms for stochastic dynamics in the
3999 ! eigenspace of mass-scaled UNRES friction matrix
4003 write (iout,*) "i",i," fricgam",fricgam(i)
4004 gdt = fricgam(i) * d_time
4006 ! Stochastic dynamics reduces to simple MD for zero friction
4008 if (gdt .le. zero) then
4009 pfric_vec(i) = 1.0d0
4010 vfric_vec(i) = d_time
4011 afric_vec(i) = 0.5d0*d_time*d_time
4012 prand_vec(i) = afric_vec(i)
4013 vrand_vec2(i) = vfric_vec(i)
4015 ! Analytical expressions when friction coefficient is large
4020 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
4021 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
4022 prand_vec(i) = afric_vec(i)
4023 vrand_vec2(i) = vfric_vec(i)
4025 ! Compute the scaling factors of random terms for the nonzero friction case
4027 ! ktm = 0.5d0*d_time/fricgam(i)
4028 ! psig = dsqrt(ktm*pterm) / fricgam(i)
4029 ! vsig = dsqrt(ktm*vterm)
4030 ! prand_vec(i) = psig*afric_vec(i)
4031 ! vrand_vec2(i) = vsig*vfric_vec(i)
4036 "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",&
4039 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),&
4040 afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
4044 ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
4046 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
4047 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
4048 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
4049 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
4050 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
4052 t_sdsetup=t_sdsetup+MPI_Wtime()
4054 t_sdsetup=t_sdsetup+tcpu()-tt0
4057 end subroutine sd_verlet_ciccotti_setup
4058 !-----------------------------------------------------------------------------
4059 subroutine sd_verlet1_ciccotti
4061 ! Applying stochastic velocity Verlet algorithm - step 1 to velocities
4062 ! implicit real*8 (a-h,o-z)
4064 ! include 'DIMENSIONS'
4068 ! include 'COMMON.CONTROL'
4069 ! include 'COMMON.VAR'
4070 ! include 'COMMON.MD'
4072 ! include 'COMMON.LANGEVIN'
4074 ! include 'COMMON.LANGEVIN.lang0'
4076 ! include 'COMMON.CHAIN'
4077 ! include 'COMMON.DERIV'
4078 ! include 'COMMON.GEO'
4079 ! include 'COMMON.LOCAL'
4080 ! include 'COMMON.INTERACT'
4081 ! include 'COMMON.IOUNITS'
4082 ! include 'COMMON.NAMES'
4083 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
4084 !el common /stochcalc/ stochforcvec
4085 logical :: lprn = .false.
4086 real(kind=8) :: ddt1,ddt2
4087 integer :: i,j,ind,inres
4088 ! write (iout,*) "dc_old"
4090 ! write (iout,'(i5,3f10.5,5x,3f10.5)')
4091 ! & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
4094 dc_work(j)=dc_old(j,0)
4095 d_t_work(j)=d_t_old(j,0)
4096 d_a_work(j)=d_a_old(j,0)
4101 dc_work(ind+j)=dc_old(j,i)
4102 d_t_work(ind+j)=d_t_old(j,i)
4103 d_a_work(ind+j)=d_a_old(j,i)
4108 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
4110 dc_work(ind+j)=dc_old(j,i+nres)
4111 d_t_work(ind+j)=d_t_old(j,i+nres)
4112 d_a_work(ind+j)=d_a_old(j,i+nres)
4121 "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",&
4125 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
4126 vfric_mat(i,j),afric_mat(i,j),&
4127 prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
4135 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) &
4136 +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
4137 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
4138 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
4140 d_t_work_new(i)=ddt1+0.5d0*ddt2
4141 d_t_work(i)=ddt1+ddt2
4146 d_t(j,0)=d_t_work(j)
4151 dc(j,i)=dc_work(ind+j)
4152 d_t(j,i)=d_t_work(ind+j)
4158 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
4159 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4160 .and.(mnum.lt.4)) then
4163 dc(j,inres)=dc_work(ind+j)
4164 d_t(j,inres)=d_t_work(ind+j)
4170 end subroutine sd_verlet1_ciccotti
4171 !-----------------------------------------------------------------------------
4172 subroutine sd_verlet2_ciccotti
4174 ! Calculating the adjusted velocities for accelerations
4176 ! implicit real*8 (a-h,o-z)
4177 ! include 'DIMENSIONS'
4178 ! include 'COMMON.CONTROL'
4179 ! include 'COMMON.VAR'
4180 ! include 'COMMON.MD'
4182 ! include 'COMMON.LANGEVIN'
4184 ! include 'COMMON.LANGEVIN.lang0'
4186 ! include 'COMMON.CHAIN'
4187 ! include 'COMMON.DERIV'
4188 ! include 'COMMON.GEO'
4189 ! include 'COMMON.LOCAL'
4190 ! include 'COMMON.INTERACT'
4191 ! include 'COMMON.IOUNITS'
4192 ! include 'COMMON.NAMES'
4193 !el real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV !(MAXRES6) maxres6=6*maxres
4194 real(kind=8),dimension(6*nres) :: stochforcvecV !(MAXRES6) maxres6=6*maxres
4195 !el common /stochcalc/ stochforcvec
4196 real(kind=8) :: ddt1,ddt2
4197 integer :: i,j,ind,inres
4199 ! Compute the stochastic forces which contribute to velocity change
4201 call stochastic_force(stochforcvecV)
4208 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
4209 ! ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
4210 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
4212 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
4216 d_t(j,0)=d_t_work(j)
4221 d_t(j,i)=d_t_work(ind+j)
4227 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4229 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
4232 d_t(j,inres)=d_t_work(ind+j)
4238 end subroutine sd_verlet2_ciccotti
4240 !-----------------------------------------------------------------------------
4242 !-----------------------------------------------------------------------------
4243 subroutine inertia_tensor
4245 ! Calculating the intertia tensor for the entire protein in order to
4246 ! remove the perpendicular components of velocity matrix which cause
4247 ! the molecule to rotate.
4250 ! implicit real*8 (a-h,o-z)
4251 ! include 'DIMENSIONS'
4252 ! include 'COMMON.CONTROL'
4253 ! include 'COMMON.VAR'
4254 ! include 'COMMON.MD'
4255 ! include 'COMMON.CHAIN'
4256 ! include 'COMMON.DERIV'
4257 ! include 'COMMON.GEO'
4258 ! include 'COMMON.LOCAL'
4259 ! include 'COMMON.INTERACT'
4260 ! include 'COMMON.IOUNITS'
4261 ! include 'COMMON.NAMES'
4263 real(kind=8),dimension(3,3) :: Im,Imcp,eigvec,Id
4264 real(kind=8),dimension(3) :: pr,eigval,L,vp,vrot
4265 real(kind=8) :: M_SC,mag,mag2,M_PEP
4266 real(kind=8),dimension(3,0:nres) :: vpp !(3,0:MAXRES)
4267 real(kind=8),dimension(3) :: vs_p,pp,incr,v
4268 real(kind=8),dimension(3,3) :: pr1,pr2
4270 !el common /gucio/ cm
4271 integer :: iti,inres,i,j,k,mnum
4282 ! calculating the center of the mass of the protein
4286 if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
4287 write(iout,*) "WTF",itype(i,mnum),i,mnum,mp(mnum)
4288 ! if (itype(i,mnum).eq.ntyp1_molec(mnum)) cycle
4289 M_PEP=M_PEP+mp(mnum)
4292 cm(j)=cm(j)+(c(j,i)+0.5d0*dc(j,i))*mp(mnum)
4301 if (mnum.ge.5) cycle
4302 iti=iabs(itype(i,mnum))
4303 M_SC=M_SC+msc(iabs(iti),mnum)
4306 cm(j)=cm(j)+msc(iabs(iti),mnum)*c(j,inres)
4310 cm(j)=cm(j)/(M_SC+M_PEP)
4312 ! write(iout,*) "Center of mass:",cm
4315 if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
4317 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
4319 Im(1,1)=Im(1,1)+mp(mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
4320 Im(1,2)=Im(1,2)-mp(mnum)*pr(1)*pr(2)
4321 Im(1,3)=Im(1,3)-mp(mnum)*pr(1)*pr(3)
4322 Im(2,3)=Im(2,3)-mp(mnum)*pr(2)*pr(3)
4323 Im(2,2)=Im(2,2)+mp(mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
4324 Im(3,3)=Im(3,3)+mp(mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
4327 ! write(iout,*) "The angular momentum before msc add"
4329 ! write (iout,*) (Im(i,j),j=1,3)
4334 iti=iabs(itype(i,mnum))
4335 if (mnum.ge.5) cycle
4338 pr(j)=c(j,inres)-cm(j)
4340 Im(1,1)=Im(1,1)+msc(iabs(iti),mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
4341 Im(1,2)=Im(1,2)-msc(iabs(iti),mnum)*pr(1)*pr(2)
4342 Im(1,3)=Im(1,3)-msc(iabs(iti),mnum)*pr(1)*pr(3)
4343 Im(2,3)=Im(2,3)-msc(iabs(iti),mnum)*pr(2)*pr(3)
4344 Im(2,2)=Im(2,2)+msc(iabs(iti),mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
4345 Im(3,3)=Im(3,3)+msc(iabs(iti),mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
4347 ! write(iout,*) "The angular momentum before Ip add"
4349 ! write (iout,*) (Im(i,j),j=1,3)
4354 Im(1,1)=Im(1,1)+Ip(mnum)*(1-dc_norm(1,i)*dc_norm(1,i))* &
4355 vbld(i+1)*vbld(i+1)*0.25d0
4356 Im(1,2)=Im(1,2)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(2,i))* &
4357 vbld(i+1)*vbld(i+1)*0.25d0
4358 Im(1,3)=Im(1,3)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(3,i))* &
4359 vbld(i+1)*vbld(i+1)*0.25d0
4360 Im(2,3)=Im(2,3)+Ip(mnum)*(-dc_norm(2,i)*dc_norm(3,i))* &
4361 vbld(i+1)*vbld(i+1)*0.25d0
4362 Im(2,2)=Im(2,2)+Ip(mnum)*(1-dc_norm(2,i)*dc_norm(2,i))* &
4363 vbld(i+1)*vbld(i+1)*0.25d0
4364 Im(3,3)=Im(3,3)+Ip(mnum)*(1-dc_norm(3,i)*dc_norm(3,i))* &
4365 vbld(i+1)*vbld(i+1)*0.25d0
4367 ! write(iout,*) "The angular momentum before Isc add"
4369 ! write (iout,*) (Im(i,j),j=1,3)
4375 ! if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
4376 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4377 .and.(mnum.lt.4)) then
4378 iti=iabs(itype(i,mnum))
4380 Im(1,1)=Im(1,1)+Isc(iti,mnum)*(1-dc_norm(1,inres)* &
4381 dc_norm(1,inres))*vbld(inres)*vbld(inres)
4382 Im(1,2)=Im(1,2)-Isc(iti,mnum)*(dc_norm(1,inres)* &
4383 dc_norm(2,inres))*vbld(inres)*vbld(inres)
4384 Im(1,3)=Im(1,3)-Isc(iti,mnum)*(dc_norm(1,inres)* &
4385 dc_norm(3,inres))*vbld(inres)*vbld(inres)
4386 Im(2,3)=Im(2,3)-Isc(iti,mnum)*(dc_norm(2,inres)* &
4387 dc_norm(3,inres))*vbld(inres)*vbld(inres)
4388 Im(2,2)=Im(2,2)+Isc(iti,mnum)*(1-dc_norm(2,inres)* &
4389 dc_norm(2,inres))*vbld(inres)*vbld(inres)
4390 Im(3,3)=Im(3,3)+Isc(iti,mnum)*(1-dc_norm(3,inres)* &
4391 dc_norm(3,inres))*vbld(inres)*vbld(inres)
4395 ! write(iout,*) "The angular momentum before agnom:"
4397 ! write (iout,*) (Im(i,j),j=1,3)
4401 ! write(iout,*) "The angular momentum before adjustment:"
4402 ! write(iout,*) (L(j),j=1,3)
4404 ! write (iout,*) (Im(i,j),j=1,3)
4410 ! Copying the Im matrix for the djacob subroutine
4418 ! Finding the eigenvectors and eignvalues of the inertia tensor
4419 call djacob(3,3,10000,1.0d-10,Imcp,eigvec,eigval)
4420 ! write (iout,*) "Eigenvalues & Eigenvectors"
4421 ! write (iout,'(5x,3f10.5)') (eigval(i),i=1,3)
4424 ! write (iout,'(i5,3f10.5)') i,(eigvec(i,j),j=1,3)
4426 ! Constructing the diagonalized matrix
4428 if (dabs(eigval(i)).gt.1.0d-15) then
4429 Id(i,i)=1.0d0/eigval(i)
4436 Imcp(i,j)=eigvec(j,i)
4442 pr1(i,j)=pr1(i,j)+Id(i,k)*Imcp(k,j)
4449 pr2(i,j)=pr2(i,j)+eigvec(i,k)*pr1(k,j)
4453 ! Calculating the total rotational velocity of the molecule
4456 vrot(i)=vrot(i)+pr2(i,j)*L(j)
4459 ! Resetting the velocities
4461 call vecpr(vrot(1),dc(1,i),vp)
4463 ! print *,"HERE2",d_t(j,i),vp(j)
4464 d_t(j,i)=d_t(j,i)-vp(j)
4465 ! print *,"HERE2",d_t(j,i),vp(j)
4470 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4471 .and.(mnum.lt.4)) then
4472 ! if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
4473 ! if(itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
4475 call vecpr(vrot(1),dc(1,inres),vp)
4477 d_t(j,inres)=d_t(j,inres)-vp(j)
4482 ! write(iout,*) "The angular momentum after adjustment:"
4483 ! write(iout,*) (L(j),j=1,3)
4486 end subroutine inertia_tensor
4487 !-----------------------------------------------------------------------------
4488 subroutine angmom(cm,L)
4491 ! implicit real*8 (a-h,o-z)
4492 ! include 'DIMENSIONS'
4493 ! include 'COMMON.CONTROL'
4494 ! include 'COMMON.VAR'
4495 ! include 'COMMON.MD'
4496 ! include 'COMMON.CHAIN'
4497 ! include 'COMMON.DERIV'
4498 ! include 'COMMON.GEO'
4499 ! include 'COMMON.LOCAL'
4500 ! include 'COMMON.INTERACT'
4501 ! include 'COMMON.IOUNITS'
4502 ! include 'COMMON.NAMES'
4503 real(kind=8) :: mscab
4504 real(kind=8),dimension(3) :: L,cm,pr,vp,vrot,incr,v,pp
4505 integer :: iti,inres,i,j,mnum
4506 ! Calculate the angular momentum
4515 if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
4517 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
4520 v(j)=incr(j)+0.5d0*d_t(j,i)
4523 incr(j)=incr(j)+d_t(j,i)
4525 call vecpr(pr(1),v(1),vp)
4527 L(j)=L(j)+mp(mnum)*vp(j)
4528 ! print *,"HERE3",J,i,L(j),mp(mnum),Ip(mnum),mnum
4532 pp(j)=0.5d0*d_t(j,i)
4534 call vecpr(pr(1),pp(1),vp)
4537 L(j)=L(j)+Ip(mnum)*vp(j)
4545 iti=iabs(itype(i,mnum))
4553 pr(j)=c(j,inres)-cm(j)
4556 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4557 .and.(mnum.lt.4)) then
4559 v(j)=incr(j)+d_t(j,inres)
4566 ! print *,i,pr,"pr",v
4567 call vecpr(pr(1),v(1),vp)
4568 ! write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),&
4569 ! " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
4571 L(j)=L(j)+mscab*vp(j)
4573 ! write (iout,*) "L",(l(j),j=1,3)
4574 ! print *,"L",(l(j),j=1,3),i,vp(1)
4576 if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4577 .and.(mnum.lt.4)) then
4579 v(j)=incr(j)+d_t(j,inres)
4581 call vecpr(dc(1,inres),d_t(1,inres),vp)
4583 L(j)=L(j)+Isc(iti,mnum)*vp(j)
4587 incr(j)=incr(j)+d_t(j,i)
4591 end subroutine angmom
4592 !-----------------------------------------------------------------------------
4593 subroutine vcm_vel(vcm)
4596 ! implicit real*8 (a-h,o-z)
4597 ! include 'DIMENSIONS'
4598 ! include 'COMMON.VAR'
4599 ! include 'COMMON.MD'
4600 ! include 'COMMON.CHAIN'
4601 ! include 'COMMON.DERIV'
4602 ! include 'COMMON.GEO'
4603 ! include 'COMMON.LOCAL'
4604 ! include 'COMMON.INTERACT'
4605 ! include 'COMMON.IOUNITS'
4606 real(kind=8),dimension(3) :: vcm,vv
4607 real(kind=8) :: summas,amas
4617 if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
4619 summas=summas+mp(mnum)
4621 vcm(j)=vcm(j)+mp(mnum)*(vv(j)+0.5d0*d_t(j,i))
4622 ! print *,i,j,vv(j),d_t(j,i)
4626 amas=msc(iabs(itype(i,mnum)),mnum)
4631 if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4632 .and.(mnum.lt.4)) then
4634 vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres))
4638 vcm(j)=vcm(j)+amas*vv(j)
4642 vv(j)=vv(j)+d_t(j,i)
4645 ! write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas
4647 vcm(j)=vcm(j)/summas
4650 end subroutine vcm_vel
4651 !-----------------------------------------------------------------------------
4653 !-----------------------------------------------------------------------------
4655 ! RATTLE algorithm for velocity Verlet - step 1, UNRES
4659 ! implicit real*8 (a-h,o-z)
4660 ! include 'DIMENSIONS'
4662 ! include 'COMMON.CONTROL'
4663 ! include 'COMMON.VAR'
4664 ! include 'COMMON.MD'
4666 ! include 'COMMON.LANGEVIN'
4668 ! include 'COMMON.LANGEVIN.lang0'
4670 ! include 'COMMON.CHAIN'
4671 ! include 'COMMON.DERIV'
4672 ! include 'COMMON.GEO'
4673 ! include 'COMMON.LOCAL'
4674 ! include 'COMMON.INTERACT'
4675 ! include 'COMMON.IOUNITS'
4676 ! include 'COMMON.NAMES'
4677 ! include 'COMMON.TIME1'
4678 !el real(kind=8) :: gginv(2*nres,2*nres),&
4679 !el gdc(3,2*nres,2*nres)
4680 real(kind=8) :: dC_uncor(3,2*nres) !,&
4681 !el real(kind=8) :: Cmat(2*nres,2*nres)
4682 real(kind=8) :: x(2*nres),xcorr(3,2*nres) !maxres2=2*maxres
4683 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
4684 !el common /przechowalnia/ nbond
4685 integer :: max_rattle = 5
4686 logical :: lprn = .false., lprn1 = .false., not_done
4687 real(kind=8) :: tol_rattle = 1.0d-5
4689 integer :: ii,i,j,jj,l,ind,ind1,nres2
4692 !el /common/ przechowalnia
4694 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4695 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4696 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4698 if (lprn) write (iout,*) "RATTLE1"
4702 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4703 .and.(mnum.lt.4)) nbond=nbond+1
4705 ! Make a folded form of the Ginv-matrix
4718 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4723 if (itype(k,1).ne.10 .and. itype(k,mnum).ne.ntyp1_molec(mnum)&
4724 .and.(mnum.lt.4)) then
4728 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4736 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4747 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4751 if (itype(k,1).ne.10) then
4755 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4763 write (iout,*) "Matrix GGinv"
4764 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
4770 if (iter.gt.max_rattle) then
4771 write (iout,*) "Error - too many iterations in RATTLE."
4774 ! Calculate the matrix C = GG**(-1) dC_old o dC
4779 dC_uncor(j,ind1)=dC(j,i)
4783 if (itype(i,1).ne.10) then
4786 dC_uncor(j,ind1)=dC(j,i+nres)
4795 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
4799 if (itype(k,1).ne.10) then
4802 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
4807 ! Calculate deviations from standard virtual-bond lengths
4811 x(ind)=vbld(i+1)**2-vbl**2
4814 if (itype(i,1).ne.10) then
4816 x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4820 write (iout,*) "Coordinates and violations"
4822 write(iout,'(i5,3f10.5,5x,e15.5)') &
4823 i,(dC_uncor(j,i),j=1,3),x(i)
4825 write (iout,*) "Velocities and violations"
4829 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4830 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4834 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4835 .and.(mnum.lt.4)) then
4838 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4839 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4840 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4843 ! write (iout,*) "gdc"
4845 ! write (iout,*) "i",i
4847 ! write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4853 if (dabs(x(i)).gt.xmax) then
4857 if (xmax.lt.tol_rattle) then
4861 ! Calculate the matrix of the system of equations
4866 Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
4871 write (iout,*) "Matrix Cmat"
4872 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4874 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4875 ! Add constraint term to positions
4882 xx = xx+x(ii)*gdc(j,ind,ii)
4886 d_t_new(j,i)=d_t_new(j,i)-xx/d_time
4890 if (itype(i,1).ne.10) then
4895 xx = xx+x(ii)*gdc(j,ind,ii)
4898 dC(j,i+nres)=dC(j,i+nres)-xx
4899 d_t_new(j,i+nres)=d_t_new(j,i+nres)-xx/d_time
4903 ! Rebuild the chain using the new coordinates
4904 call chainbuild_cart
4906 write (iout,*) "New coordinates, Lagrange multipliers,",&
4907 " and differences between actual and standard bond lengths"
4911 xx=vbld(i+1)**2-vbl**2
4912 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4913 i,(dC(j,i),j=1,3),x(ind),xx
4917 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4920 xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4921 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4922 i,(dC(j,i+nres),j=1,3),x(ind),xx
4925 write (iout,*) "Velocities and violations"
4929 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4930 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4933 if (itype(i,1).ne.10) then
4935 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4936 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4937 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4944 10 write (iout,*) "Error - singularity in solving the system",&
4945 " of equations for Lagrange multipliers."
4949 "RATTLE inactive; use -DRATTLE switch at compile time."
4952 end subroutine rattle1
4953 !-----------------------------------------------------------------------------
4955 ! RATTLE algorithm for velocity Verlet - step 2, UNRES
4959 ! implicit real*8 (a-h,o-z)
4960 ! include 'DIMENSIONS'
4962 ! include 'COMMON.CONTROL'
4963 ! include 'COMMON.VAR'
4964 ! include 'COMMON.MD'
4966 ! include 'COMMON.LANGEVIN'
4968 ! include 'COMMON.LANGEVIN.lang0'
4970 ! include 'COMMON.CHAIN'
4971 ! include 'COMMON.DERIV'
4972 ! include 'COMMON.GEO'
4973 ! include 'COMMON.LOCAL'
4974 ! include 'COMMON.INTERACT'
4975 ! include 'COMMON.IOUNITS'
4976 ! include 'COMMON.NAMES'
4977 ! include 'COMMON.TIME1'
4978 !el real(kind=8) :: gginv(2*nres,2*nres),&
4979 !el gdc(3,2*nres,2*nres)
4980 real(kind=8) :: dC_uncor(3,2*nres) !,&
4981 !el Cmat(2*nres,2*nres)
4982 real(kind=8) :: x(2*nres) !maxres2=2*maxres
4983 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
4984 !el common /przechowalnia/ nbond
4985 integer :: max_rattle = 5
4986 logical :: lprn = .false., lprn1 = .false., not_done
4987 real(kind=8) :: tol_rattle = 1.0d-5
4991 !el /common/ przechowalnia
4992 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4993 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4994 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4996 if (lprn) write (iout,*) "RATTLE2"
4997 if (lprn) write (iout,*) "Velocity correction"
4998 ! Calculate the matrix G dC
5004 gdc(j,i,ind)=GGinv(i,ind)*dC(j,k)
5009 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5010 .and.(mnum.lt.4)) then
5013 gdc(j,i,ind)=GGinv(i,ind)*dC(j,k+nres)
5019 ! write (iout,*) "gdc"
5021 ! write (iout,*) "i",i
5023 ! write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
5027 ! Calculate the matrix of the system of equations
5034 Cmat(ind,j)=Cmat(ind,j)+dC(k,i)*gdc(k,ind,j)
5040 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5041 .and.(mnum.lt.4)) then
5046 Cmat(ind,j)=Cmat(ind,j)+dC(k,i+nres)*gdc(k,ind,j)
5051 ! Calculate the scalar product dC o d_t_new
5055 x(ind)=scalar(d_t(1,i),dC(1,i))
5059 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5060 .and.(mnum.lt.4)) then
5062 x(ind)=scalar(d_t(1,i+nres),dC(1,i+nres))
5066 write (iout,*) "Velocities and violations"
5070 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5071 i,ind,(d_t(j,i),j=1,3),x(ind)
5075 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5076 .and.(mnum.lt.4)) then
5078 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5079 i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind)
5085 if (dabs(x(i)).gt.xmax) then
5089 if (xmax.lt.tol_rattle) then
5094 write (iout,*) "Matrix Cmat"
5095 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
5097 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
5098 ! Add constraint term to velocities
5105 xx = xx+x(ii)*gdc(j,ind,ii)
5107 d_t(j,i)=d_t(j,i)-xx
5112 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5113 .and.(mnum.lt.4)) then
5118 xx = xx+x(ii)*gdc(j,ind,ii)
5120 d_t(j,i+nres)=d_t(j,i+nres)-xx
5126 "New velocities, Lagrange multipliers violations"
5130 if (lprn) write (iout,'(2i5,3f10.5,5x,2e15.5)') &
5131 i,ind,(d_t(j,i),j=1,3),x(ind),scalar(d_t(1,i),dC(1,i))
5135 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
5138 write (iout,'(2i5,3f10.5,5x,2e15.5)') &
5139 i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind),&
5140 scalar(d_t(1,i+nres),dC(1,i+nres))
5146 10 write (iout,*) "Error - singularity in solving the system",&
5147 " of equations for Lagrange multipliers."
5151 "RATTLE inactive; use -DRATTLE option at compile time."
5154 end subroutine rattle2
5155 !-----------------------------------------------------------------------------
5156 subroutine rattle_brown
5157 ! RATTLE/LINCS algorithm for Brownian dynamics, UNRES
5161 ! implicit real*8 (a-h,o-z)
5162 ! include 'DIMENSIONS'
5164 ! include 'COMMON.CONTROL'
5165 ! include 'COMMON.VAR'
5166 ! include 'COMMON.MD'
5168 ! include 'COMMON.LANGEVIN'
5170 ! include 'COMMON.LANGEVIN.lang0'
5172 ! include 'COMMON.CHAIN'
5173 ! include 'COMMON.DERIV'
5174 ! include 'COMMON.GEO'
5175 ! include 'COMMON.LOCAL'
5176 ! include 'COMMON.INTERACT'
5177 ! include 'COMMON.IOUNITS'
5178 ! include 'COMMON.NAMES'
5179 ! include 'COMMON.TIME1'
5180 !el real(kind=8) :: gginv(2*nres,2*nres),&
5181 !el gdc(3,2*nres,2*nres)
5182 real(kind=8) :: dC_uncor(3,2*nres) !,&
5183 !el real(kind=8) :: Cmat(2*nres,2*nres)
5184 real(kind=8) :: x(2*nres) !maxres2=2*maxres
5185 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
5186 !el common /przechowalnia/ nbond
5187 integer :: max_rattle = 5
5188 logical :: lprn = .false., lprn1 = .false., not_done
5189 real(kind=8) :: tol_rattle = 1.0d-5
5193 !el /common/ przechowalnia
5194 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
5195 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
5196 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
5199 if (lprn) write (iout,*) "RATTLE_BROWN"
5202 if (itype(i,1).ne.10) nbond=nbond+1
5204 ! Make a folded form of the Ginv-matrix
5217 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
5221 if (itype(k,1).ne.10) then
5225 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
5232 if (itype(i,1).ne.10) then
5242 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
5246 if (itype(k,1).ne.10) then
5250 if (j.eq.1 .and. l.eq.1)GGinv(ii,jj)=fricmat(ind,ind1)
5258 write (iout,*) "Matrix GGinv"
5259 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
5265 if (iter.gt.max_rattle) then
5266 write (iout,*) "Error - too many iterations in RATTLE."
5269 ! Calculate the matrix C = GG**(-1) dC_old o dC
5274 dC_uncor(j,ind1)=dC(j,i)
5278 if (itype(i,1).ne.10) then
5281 dC_uncor(j,ind1)=dC(j,i+nres)
5290 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
5294 if (itype(k,1).ne.10) then
5297 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
5302 ! Calculate deviations from standard virtual-bond lengths
5306 x(ind)=vbld(i+1)**2-vbl**2
5309 if (itype(i,1).ne.10) then
5311 x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
5315 write (iout,*) "Coordinates and violations"
5317 write(iout,'(i5,3f10.5,5x,e15.5)') &
5318 i,(dC_uncor(j,i),j=1,3),x(i)
5320 write (iout,*) "Velocities and violations"
5324 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5325 i,ind,(d_t(j,i),j=1,3),scalar(d_t(1,i),dC_old(1,i))
5328 if (itype(i,1).ne.10) then
5330 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5331 i+nres,ind,(d_t(j,i+nres),j=1,3),&
5332 scalar(d_t(1,i+nres),dC_old(1,i+nres))
5335 write (iout,*) "gdc"
5337 write (iout,*) "i",i
5339 write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
5345 if (dabs(x(i)).gt.xmax) then
5349 if (xmax.lt.tol_rattle) then
5353 ! Calculate the matrix of the system of equations
5358 Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
5363 write (iout,*) "Matrix Cmat"
5364 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
5366 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
5367 ! Add constraint term to positions
5374 xx = xx+x(ii)*gdc(j,ind,ii)
5377 d_t(j,i)=d_t(j,i)+xx/d_time
5382 if (itype(i,1).ne.10) then
5387 xx = xx+x(ii)*gdc(j,ind,ii)
5390 d_t(j,i+nres)=d_t(j,i+nres)+xx/d_time
5391 dC(j,i+nres)=dC(j,i+nres)+xx
5395 ! Rebuild the chain using the new coordinates
5396 call chainbuild_cart
5398 write (iout,*) "New coordinates, Lagrange multipliers,",&
5399 " and differences between actual and standard bond lengths"
5403 xx=vbld(i+1)**2-vbl**2
5404 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
5405 i,(dC(j,i),j=1,3),x(ind),xx
5408 if (itype(i,1).ne.10) then
5410 xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
5411 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
5412 i,(dC(j,i+nres),j=1,3),x(ind),xx
5415 write (iout,*) "Velocities and violations"
5419 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5420 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
5423 if (itype(i,1).ne.10) then
5425 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5426 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
5427 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
5434 10 write (iout,*) "Error - singularity in solving the system",&
5435 " of equations for Lagrange multipliers."
5439 "RATTLE inactive; use -DRATTLE option at compile time"
5442 end subroutine rattle_brown
5443 !-----------------------------------------------------------------------------
5445 !-----------------------------------------------------------------------------
5446 subroutine friction_force
5451 ! implicit real*8 (a-h,o-z)
5452 ! include 'DIMENSIONS'
5453 ! include 'COMMON.VAR'
5454 ! include 'COMMON.CHAIN'
5455 ! include 'COMMON.DERIV'
5456 ! include 'COMMON.GEO'
5457 ! include 'COMMON.LOCAL'
5458 ! include 'COMMON.INTERACT'
5459 ! include 'COMMON.MD'
5461 ! include 'COMMON.LANGEVIN'
5463 ! include 'COMMON.LANGEVIN.lang0'
5465 ! include 'COMMON.IOUNITS'
5466 !el real(kind=8),dimension(6*nres) :: gamvec !(MAXRES6) maxres6=6*maxres
5467 !el common /syfek/ gamvec
5469 integer iposc,ichain,n,innt,inct
5470 double precision rs(nres+2)
5473 real(kind=8) :: vv(3),vvtot(3,nres),v_work(6*nres) !,&
5474 !el ginvfric(2*nres,2*nres) !maxres2=2*maxres
5475 !el common /przechowalnia/ ginvfric
5477 logical :: lprn, checkmode
5478 integer :: i,j,ind,k,nres2,nres6,mnum
5483 ! if (large) lprn=.true.
5484 ! if (large) checkmode=.true.
5486 c Here accelerations due to friction forces are computed right after forces.
5487 d_t_work(:6*nres)=0.0d0
5489 v_work(j,1)=d_t(j,0)
5490 v_work(j,nnt)=d_t(j,0)
5494 v_work(j,i)=v_work(j,i-1)+d_t(j,i-1)
5499 if (iabs(itype(i,1)).ne.10 .and. iabs(itype(i,mnum)).ne.ntyp1_molec(mnum).and.mnum.lt.3) then
5501 v_work(j,i+nres)=v_work(j,i)+d_t(j,i+nres)
5506 write (iout,*) "v_work"
5508 write (iout,'(i5,3f10.5)') i,(v_work(j,i),j=1,3)
5514 n=dimen_chain(ichain)
5515 iposc=iposd_chain(ichain)
5516 c write (iout,*) "friction_force j",j," ichain",ichain,
5517 c & " n",n," iposc",iposc,iposc+n-1
5518 innt=chain_border(1,ichain)
5519 inct=chain_border(2,ichain)
5521 c innt=chain_border(1,1)
5522 c inct=chain_border(2,1)
5524 vvec(ind+1)=v_work(j,i)
5526 ! if (iabs(itype(i)).ne.10) then
5527 if (iabs(itype(i,1)).ne.10 .and. iabs(itype(i,mnum)).ne.ntyp1_molec(mnum).and.mnum.lt.3) then
5528 vvec(ind+1)=v_work(j,i+nres)
5533 write (iout,*) "vvec ind",ind," n",n
5534 write (iout,'(f10.5)') (vvec(i),i=iposc,ind)
5536 c write (iout,*) "chain",i," ind",ind," n",n
5544 call fivediagmult(n,DMfric(iposc),DU1fric(iposc),
5545 & DU2fric(iposc),vvec(iposc),rs)
5548 time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
5550 time_fricmatmult=time_fricmatmult+tcpu()-time01
5555 write (iout,'(f10.5)') (rs(i),i=1,n)
5557 do i=iposc,iposc+n-1
5558 c write (iout,*) "ichain",ichain," i",i," j",j,
5559 c & "index",3*(i-1)+j,"rs",rs(i-iposc+1)
5560 fric_work(3*(i-1)+j)=-rs(i-iposc+1)
5565 write (iout,*) "Vector fric_work dimen3",dimen3
5566 write (iout,'(3f10.5)') (fric_work(j),j=1,dimen3)
5569 if(.not.allocated(gamvec)) allocate(gamvec(nres6)) !(MAXRES6)
5570 if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
5578 d_t_work(j)=d_t(j,0)
5583 d_t_work(ind+j)=d_t(j,i)
5589 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5590 .and.(mnum.lt.4)) then
5592 d_t_work(ind+j)=d_t(j,i+nres)
5598 call fricmat_mult(d_t_work,fric_work)
5600 if (.not.checkmode) return
5603 write (iout,*) "d_t_work and fric_work"
5605 write (iout,'(i3,2e15.5)') i,d_t_work(i),fric_work(i)
5609 friction(j,0)=fric_work(j)
5614 friction(j,i)=fric_work(ind+j)
5620 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5621 .and.(mnum.lt.4)) then
5622 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
5624 friction(j,i+nres)=fric_work(ind+j)
5630 write(iout,*) "Friction backbone"
5632 write(iout,'(i5,3e15.5,5x,3e15.5)') &
5633 i,(friction(j,i),j=1,3),(d_t(j,i),j=1,3)
5635 write(iout,*) "Friction side chain"
5637 write(iout,'(i5,3e15.5,5x,3e15.5)') &
5638 i,(friction(j,i+nres),j=1,3),(d_t(j,i+nres),j=1,3)
5647 vvtot(j,i)=vv(j)+0.5d0*d_t(j,i)
5648 vvtot(j,i+nres)=vv(j)+d_t(j,i+nres)
5649 vv(j)=vv(j)+d_t(j,i)
5652 write (iout,*) "vvtot backbone and sidechain"
5654 write (iout,'(i5,3e15.5,5x,3e15.5)') i,(vvtot(j,i),j=1,3),&
5655 (vvtot(j,i+nres),j=1,3)
5660 v_work(ind+j)=vvtot(j,i)
5666 v_work(ind+j)=vvtot(j,i+nres)
5670 write (iout,*) "v_work gamvec and site-based friction forces"
5672 write (iout,'(i5,3e15.5)') i,v_work(i),gamvec(i),&
5676 ! fric_work1(i)=0.0d0
5678 ! fric_work1(i)=fric_work1(i)-A(j,i)*gamvec(j)*v_work(j)
5681 ! write (iout,*) "fric_work and fric_work1"
5683 ! write (iout,'(i5,2e15.5)') i,fric_work(i),fric_work1(i)
5689 ginvfric(i,j)=ginvfric(i,j)+ginv(i,k)*fricmat(k,j)
5693 write (iout,*) "ginvfric"
5695 write (iout,'(i5,100f8.3)') i,(ginvfric(i,j),j=1,dimen)
5697 write (iout,*) "symmetry check"
5700 write (iout,*) i,j,ginvfric(i,j)-ginvfric(j,i)
5706 end subroutine friction_force
5707 !-----------------------------------------------------------------------------
5708 subroutine setup_fricmat
5712 use control_data, only:time_Bcast
5713 use control, only:tcpu
5715 ! implicit real*8 (a-h,o-z)
5719 real(kind=8) :: time00
5721 ! include 'DIMENSIONS'
5722 ! include 'COMMON.VAR'
5723 ! include 'COMMON.CHAIN'
5724 ! include 'COMMON.DERIV'
5725 ! include 'COMMON.GEO'
5726 ! include 'COMMON.LOCAL'
5727 ! include 'COMMON.INTERACT'
5728 ! include 'COMMON.MD'
5729 ! include 'COMMON.SETUP'
5730 ! include 'COMMON.TIME1'
5731 ! integer licznik /0/
5734 ! include 'COMMON.LANGEVIN'
5736 ! include 'COMMON.LANGEVIN.lang0'
5738 ! include 'COMMON.IOUNITS'
5740 integer :: i,j,ind,ind1,m
5741 logical :: lprn = .false.
5742 real(kind=8) :: dtdi !el ,gamvec(2*nres)
5743 !el real(kind=8),dimension(2*nres,2*nres) :: ginvfric,fcopy
5744 ! real(kind=8),allocatable,dimension(:,:) :: fcopy
5745 !el real(kind=8),dimension(2*nres*(2*nres+1)/2) :: Ghalf !(mmaxres2) (mmaxres2=(maxres2*(maxres2+1)/2))
5746 !el common /syfek/ gamvec
5747 real(kind=8) :: work(8*2*nres)
5748 integer :: iwork(2*nres)
5749 !el common /przechowalnia/ ginvfric,Ghalf,fcopy
5750 integer :: ii,iti,k,l,nzero,nres2,nres6,ierr,mnum
5754 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
5755 if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
5756 if (fg_rank.ne.king) goto 10
5761 if(.not.allocated(gamvec)) allocate(gamvec(nres2)) !(MAXRES2)
5762 if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
5763 if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
5764 !el allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
5765 if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !maxres2=2*maxres
5767 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
5768 ! Zeroing out fricmat
5774 ! Load the friction coefficients corresponding to peptide groups
5779 gamvec(ind1)=gamp(mnum)
5782 if (molnum(nct).eq.5) then
5785 gamvec(ind1)=gamp(mnum)
5787 ! Load the friction coefficients corresponding to side chains
5791 gamsc(ntyp1_molec(j),j)=1.0d0
5798 gamvec(ii)=gamsc(iabs(iti),mnum)
5800 if (surfarea) call sdarea(gamvec)
5802 ! write (iout,*) "Matrix A and vector gamma"
5804 ! write (iout,'(i2,$)') i
5806 ! write (iout,'(f4.1,$)') A(i,j)
5808 ! write (iout,'(f8.3)') gamvec(i)
5812 write (iout,*) "Vector gamvec"
5814 write (iout,'(i5,f10.5)') i, gamvec(i)
5818 ! The friction matrix
5823 dtdi=dtdi+A(j,k)*A(j,i)*gamvec(j)
5830 write (iout,'(//a)') "Matrix fricmat"
5831 call matout2(dimen,dimen,nres2,nres2,fricmat)
5833 if (lang.eq.2 .or. lang.eq.3) then
5834 ! Mass-scale the friction matrix if non-direct integration will be performed
5840 Ginvfric(i,j)=Ginvfric(i,j)+ &
5841 Gsqrm(i,k)*Gsqrm(l,j)*fricmat(k,l)
5846 ! Diagonalize the friction matrix
5851 Ghalf(ind)=Ginvfric(i,j)
5854 call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
5857 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
5858 " mass-scaled friction matrix"
5859 call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
5861 ! Precompute matrices for tinker stochastic integrator
5868 mt1(i,j)=mt1(i,j)+fricvec(k,i)*gsqrm(k,j)
5869 mt2(i,j)=mt2(i,j)+fricvec(k,i)*gsqrp(k,j)
5875 else if (lang.eq.4) then
5876 ! Diagonalize the friction matrix
5881 Ghalf(ind)=fricmat(i,j)
5884 call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
5887 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
5889 call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
5891 ! Determine the number of zero eigenvalues of the friction matrix
5892 nzero=max0(dimen-dimen1,0)
5893 ! do while (fricgam(nzero+1).le.1.0d-5 .and. nzero.lt.dimen)
5896 write (iout,*) "Number of zero eigenvalues:",nzero
5901 fricmat(i,j)=fricmat(i,j) &
5902 +fricvec(i,k)*fricvec(j,k)/fricgam(k)
5907 write (iout,'(//a)') "Generalized inverse of fricmat"
5908 call matout(dimen,dimen,nres6,nres6,fricmat)
5913 if (nfgtasks.gt.1) then
5914 if (fg_rank.eq.0) then
5915 ! The matching BROADCAST for fg processors is called in ERGASTULUM
5921 call MPI_Bcast(10,1,MPI_INTEGER,king,FG_COMM,IERROR)
5923 time_Bcast=time_Bcast+MPI_Wtime()-time00
5925 time_Bcast=time_Bcast+tcpu()-time00
5927 ! print *,"Processor",myrank,
5928 ! & " BROADCAST iorder in SETUP_FRICMAT"
5931 write (iout,*) "setup_fricmat licznik"!,licznik !sp
5937 ! Scatter the friction matrix
5938 call MPI_Scatterv(fricmat(1,1),nginv_counts(0),&
5939 nginv_start(0),MPI_DOUBLE_PRECISION,fcopy(1,1),&
5940 myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
5943 time_scatter=time_scatter+MPI_Wtime()-time00
5944 time_scatter_fmat=time_scatter_fmat+MPI_Wtime()-time00
5946 time_scatter=time_scatter+tcpu()-time00
5947 time_scatter_fmat=time_scatter_fmat+tcpu()-time00
5951 do j=1,2*my_ng_count
5952 fricmat(j,i)=fcopy(i,j)
5955 ! write (iout,*) "My chunk of fricmat"
5956 ! call MATOUT2(my_ng_count,dimen,maxres2,maxres2,fcopy)
5960 end subroutine setup_fricmat
5961 !-----------------------------------------------------------------------------
5962 subroutine stochastic_force(stochforcvec)
5965 use random, only:anorm_distr
5966 ! implicit real*8 (a-h,o-z)
5967 ! include 'DIMENSIONS'
5968 use control, only: tcpu
5973 ! include 'COMMON.VAR'
5974 ! include 'COMMON.CHAIN'
5975 ! include 'COMMON.DERIV'
5976 ! include 'COMMON.GEO'
5977 ! include 'COMMON.LOCAL'
5978 ! include 'COMMON.INTERACT'
5979 ! include 'COMMON.MD'
5980 ! include 'COMMON.TIME1'
5982 ! include 'COMMON.LANGEVIN'
5984 ! include 'COMMON.LANGEVIN.lang0'
5986 ! include 'COMMON.IOUNITS'
5988 real(kind=8) :: x,sig,lowb,highb
5989 real(kind=8) :: ff(3),force(3,0:2*nres),zeta2,lowb2
5990 real(kind=8) :: highb2,sig2,forcvec(6*nres),stochforcvec(6*nres)
5991 real(kind=8) :: time00
5992 logical :: lprn = .false.
5993 integer :: i,j,ind,mnum
5995 integer ichain,innt,inct,iposc
6000 stochforc(j,i)=0.0d0
6010 ! Compute the stochastic forces acting on bodies. Store in force.
6016 force(j,i)=anorm_distr(x,sig,lowb,highb)
6024 force(j,i+nres)=anorm_distr(x,sig2,lowb2,highb2)
6028 time_fsample=time_fsample+MPI_Wtime()-time00
6030 time_fsample=time_fsample+tcpu()-time00
6035 innt=chain_border(1,ichain)
6036 inct=chain_border(2,ichain)
6037 iposc=iposd_chain(ichain)
6038 !c for debugging only
6039 !c innt=chain_border(1,1)
6040 !c inct=chain_border(2,1)
6041 !c iposc=iposd_chain(1)
6042 !c write (iout,*)"stochastic_force ichain=",ichain," innt",innt,
6043 !c & " inct",inct," iposc",iposc
6045 stochforcvec(ind+j)=0.5d0*force(j,innt)
6047 if (iabs(itype(innt),molnum(iint)).eq.10) then
6049 stochforcvec(ind+j)=stochforcvec(ind+j)+force(j,innt+nres)
6055 stochforcvec(ind+j)=force(j,innt+nres)
6061 stochforcvec(ind+j)=0.5d0*(force(j,i)+force(j,i-1))
6063 if (iabs(itype(i,molnum(i)).eq.10) then
6065 stochforcvec(ind+j)=stochforcvec(ind+j)+force(j,i+nres)
6071 stochforcvec(ind+j)=force(j,i+nres)
6077 stochforcvec(ind+j)=0.5d0*force(j,inct-1)
6079 if (iabs(itype(inct),molnum(inct)).eq.10) then
6081 stochforcvec(ind+j)=stochforcvec(ind+j)+force(j,inct+nres)
6087 stochforcvec(ind+j)=force(j,inct+nres)
6091 !c write (iout,*) "chain",ichain," ind",ind
6094 write (iout,*) "stochforcvec"
6095 write (iout,'(3f10.5)') (stochforcvec(j),j=1,ind)
6098 ! Compute the stochastic forces acting on virtual-bond vectors.
6104 stochforc(j,i)=ff(j)+0.5d0*force(j,i)
6107 ff(j)=ff(j)+force(j,i)
6109 ! if (itype(i+1,1).ne.ntyp1) then
6111 if (itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
6113 stochforc(j,i)=stochforc(j,i)+force(j,i+nres+1)
6114 ff(j)=ff(j)+force(j,i+nres+1)
6119 stochforc(j,0)=ff(j)+force(j,nnt+nres)
6123 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
6124 .and.(mnum.lt.4)) then
6125 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
6127 stochforc(j,i+nres)=force(j,i+nres)
6133 stochforcvec(j)=stochforc(j,0)
6138 stochforcvec(ind+j)=stochforc(j,i)
6144 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
6145 .and.(mnum.lt.4)) then
6146 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
6148 stochforcvec(ind+j)=stochforc(j,i+nres)
6154 write (iout,*) "stochforcvec"
6156 write(iout,'(i5,e15.5)') i,stochforcvec(i)
6158 write(iout,*) "Stochastic forces backbone"
6160 write(iout,'(i5,3e15.5)') i,(stochforc(j,i),j=1,3)
6162 write(iout,*) "Stochastic forces side chain"
6164 write(iout,'(i5,3e15.5)') &
6165 i,(stochforc(j,i+nres),j=1,3)
6173 write (iout,*) i,ind
6175 forcvec(ind+j)=force(j,i)
6180 write (iout,*) i,ind
6182 forcvec(j+ind)=force(j,i+nres)
6187 write (iout,*) "forcvec"
6191 write (iout,'(2i3,2f10.5)') i,j,force(j,i),&
6198 write (iout,'(2i3,2f10.5)') i,j,force(j,i+nres),&
6207 end subroutine stochastic_force
6208 !-----------------------------------------------------------------------------
6209 subroutine sdarea(gamvec)
6211 ! Scale the friction coefficients according to solvent accessible surface areas
6212 ! Code adapted from TINKER
6216 ! implicit real*8 (a-h,o-z)
6217 ! include 'DIMENSIONS'
6218 ! include 'COMMON.CONTROL'
6219 ! include 'COMMON.VAR'
6220 ! include 'COMMON.MD'
6222 ! include 'COMMON.LANGEVIN'
6224 ! include 'COMMON.LANGEVIN.lang0'
6226 ! include 'COMMON.CHAIN'
6227 ! include 'COMMON.DERIV'
6228 ! include 'COMMON.GEO'
6229 ! include 'COMMON.LOCAL'
6230 ! include 'COMMON.INTERACT'
6231 ! include 'COMMON.IOUNITS'
6232 ! include 'COMMON.NAMES'
6233 real(kind=8),dimension(2*nres) :: radius,gamvec !(maxres2)
6234 real(kind=8),parameter :: twosix = 1.122462048309372981d0
6235 logical :: lprn = .false.
6236 real(kind=8) :: probe,area,ratio
6237 integer :: i,j,ind,iti,mnum
6239 ! determine new friction coefficients every few SD steps
6241 ! set the atomic radii to estimates of sigma values
6243 ! print *,"Entered sdarea"
6249 ! Load peptide group radii
6252 radius(i)=pstok(mnum)
6254 ! Load side chain radii
6258 radius(i+nres)=restok(iti,mnum)
6261 ! write (iout,*) "i",i," radius",radius(i)
6264 radius(i) = radius(i) / twosix
6265 if (radius(i) .ne. 0.0d0) radius(i) = radius(i) + probe
6268 ! scale atomic friction coefficients by accessible area
6270 if (lprn) write (iout,*) &
6271 "Original gammas, surface areas, scaling factors, new gammas, ",&
6272 "std's of stochastic forces"
6275 if (radius(i).gt.0.0d0) then
6276 call surfatom (i,area,radius)
6277 ratio = dmax1(area/(4.0d0*pi*radius(i)**2),1.0d-1)
6278 if (lprn) write (iout,'(i5,3f10.5,$)') &
6279 i,gamvec(ind+1),area,ratio
6282 gamvec(ind) = ratio * gamvec(ind)
6285 stdforcp(i)=stdfp(mnum)*dsqrt(gamvec(ind))
6286 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcp(i)
6290 if (radius(i+nres).gt.0.0d0) then
6291 call surfatom (i+nres,area,radius)
6292 ratio = dmax1(area/(4.0d0*pi*radius(i+nres)**2),1.0d-1)
6293 if (lprn) write (iout,'(i5,3f10.5,$)') &
6294 i,gamvec(ind+1),area,ratio
6297 gamvec(ind) = ratio * gamvec(ind)
6300 stdforcsc(i)=stdfsc(itype(i,mnum),mnum)*dsqrt(gamvec(ind))
6301 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcsc(i)
6306 end subroutine sdarea
6307 !-----------------------------------------------------------------------------
6309 !-----------------------------------------------------------------------------
6312 ! ###################################################
6313 ! ## COPYRIGHT (C) 1996 by Jay William Ponder ##
6314 ! ## All Rights Reserved ##
6315 ! ###################################################
6317 ! ################################################################
6319 ! ## subroutine surfatom -- exposed surface area of an atom ##
6321 ! ################################################################
6324 ! "surfatom" performs an analytical computation of the surface
6325 ! area of a specified atom; a simplified version of "surface"
6327 ! literature references:
6329 ! T. J. Richmond, "Solvent Accessible Surface Area and
6330 ! Excluded Volume in Proteins", Journal of Molecular Biology,
6333 ! L. Wesson and D. Eisenberg, "Atomic Solvation Parameters
6334 ! Applied to Molecular Dynamics of Proteins in Solution",
6335 ! Protein Science, 1, 227-235 (1992)
6337 ! variables and parameters:
6339 ! ir number of atom for which area is desired
6340 ! area accessible surface area of the atom
6341 ! radius radii of each of the individual atoms
6344 subroutine surfatom(ir,area,radius)
6346 ! implicit real*8 (a-h,o-z)
6347 ! include 'DIMENSIONS'
6349 ! include 'COMMON.GEO'
6350 ! include 'COMMON.IOUNITS'
6352 integer :: nsup,nstart_sup
6353 ! double precision c,dc,dc_old,d_c_work,xloc,xrot,dc_norm
6354 ! common /chain/ c(3,maxres2+2),dc(3,0:maxres2),dc_old(3,0:maxres2),
6355 ! & xloc(3,maxres),xrot(3,maxres),dc_norm(3,0:maxres2),
6356 ! & dc_work(MAXRES6),nres,nres0
6357 integer,parameter :: maxarc=300
6361 integer :: mi,ni,narc
6362 integer :: key(maxarc)
6363 integer :: intag(maxarc)
6364 integer :: intag1(maxarc)
6365 real(kind=8) :: area,arcsum
6366 real(kind=8) :: arclen,exang
6367 real(kind=8) :: delta,delta2
6368 real(kind=8) :: eps,rmove
6369 real(kind=8) :: xr,yr,zr
6370 real(kind=8) :: rr,rrsq
6371 real(kind=8) :: rplus,rminus
6372 real(kind=8) :: axx,axy,axz
6373 real(kind=8) :: ayx,ayy
6374 real(kind=8) :: azx,azy,azz
6375 real(kind=8) :: uxj,uyj,uzj
6376 real(kind=8) :: tx,ty,tz
6377 real(kind=8) :: txb,tyb,td
6378 real(kind=8) :: tr2,tr,txr,tyr
6379 real(kind=8) :: tk1,tk2
6380 real(kind=8) :: thec,the,t,tb
6381 real(kind=8) :: txk,tyk,tzk
6382 real(kind=8) :: t1,ti,tf,tt
6383 real(kind=8) :: txj,tyj,tzj
6384 real(kind=8) :: ccsq,cc,xysq
6385 real(kind=8) :: bsqk,bk,cosine
6386 real(kind=8) :: dsqj,gi,pix2
6387 real(kind=8) :: therk,dk,gk
6388 real(kind=8) :: risqk,rik
6389 real(kind=8) :: radius(2*nres) !(maxatm) (maxatm=maxres2)
6390 real(kind=8) :: ri(maxarc),risq(maxarc)
6391 real(kind=8) :: ux(maxarc),uy(maxarc),uz(maxarc)
6392 real(kind=8) :: xc(maxarc),yc(maxarc),zc(maxarc)
6393 real(kind=8) :: xc1(maxarc),yc1(maxarc),zc1(maxarc)
6394 real(kind=8) :: dsq(maxarc),bsq(maxarc)
6395 real(kind=8) :: dsq1(maxarc),bsq1(maxarc)
6396 real(kind=8) :: arci(maxarc),arcf(maxarc)
6397 real(kind=8) :: ex(maxarc),lt(maxarc),gr(maxarc)
6398 real(kind=8) :: b(maxarc),b1(maxarc),bg(maxarc)
6399 real(kind=8) :: kent(maxarc),kout(maxarc)
6400 real(kind=8) :: ther(maxarc)
6401 logical :: moved,top
6402 logical :: omit(maxarc)
6405 ! maxatm = 2*nres !maxres2 maxres2=2*maxres
6406 ! maxlight = 8*maxatm
6409 ! maxtors = 4*maxatm
6411 ! zero out the surface area for the sphere of interest
6414 ! write (2,*) "ir",ir," radius",radius(ir)
6415 if (radius(ir) .eq. 0.0d0) return
6417 ! set the overlap significance and connectivity shift
6421 delta2 = delta * delta
6426 ! store coordinates and radius of the sphere of interest
6434 ! initialize values of some counters and summations
6443 ! test each sphere to see if it overlaps the sphere of interest
6446 if (i.eq.ir .or. radius(i).eq.0.0d0) goto 30
6447 rplus = rr + radius(i)
6449 if (abs(tx) .ge. rplus) goto 30
6451 if (abs(ty) .ge. rplus) goto 30
6453 if (abs(tz) .ge. rplus) goto 30
6455 ! check for sphere overlap by testing distance against radii
6457 xysq = tx*tx + ty*ty
6458 if (xysq .lt. delta2) then
6465 if (rplus-cc .le. delta) goto 30
6466 rminus = rr - radius(i)
6468 ! check to see if sphere of interest is completely buried
6470 if (cc-abs(rminus) .le. delta) then
6471 if (rminus .le. 0.0d0) goto 170
6475 ! check for too many overlaps with sphere of interest
6477 if (io .ge. maxarc) then
6479 20 format (/,' SURFATOM -- Increase the Value of MAXARC')
6483 ! get overlap between current sphere and sphere of interest
6492 gr(io) = (ccsq+rplus*rminus) / (2.0d0*rr*b1(io))
6498 ! case where no other spheres overlap the sphere of interest
6501 area = 4.0d0 * pi * rrsq
6505 ! case where only one sphere overlaps the sphere of interest
6508 area = pix2 * (1.0d0 + gr(1))
6509 area = mod(area,4.0d0*pi) * rrsq
6513 ! case where many spheres intersect the sphere of interest;
6514 ! sort the intersecting spheres by their degree of overlap
6516 call sort2 (io,gr,key)
6519 intag(i) = intag1(k)
6528 ! get radius of each overlap circle on surface of the sphere
6533 risq(i) = rrsq - gi*gi
6534 ri(i) = sqrt(risq(i))
6535 ther(i) = 0.5d0*pi - asin(min(1.0d0,max(-1.0d0,gr(i))))
6538 ! find boundary of inaccessible area on sphere of interest
6541 if (.not. omit(k)) then
6548 ! check to see if J circle is intersecting K circle;
6549 ! get distance between circle centers and sum of radii
6552 if (omit(j)) goto 60
6553 cc = (txk*xc(j)+tyk*yc(j)+tzk*zc(j))/(bk*b(j))
6554 cc = acos(min(1.0d0,max(-1.0d0,cc)))
6555 td = therk + ther(j)
6557 ! check to see if circles enclose separate regions
6559 if (cc .ge. td) goto 60
6561 ! check for circle J completely inside circle K
6563 if (cc+ther(j) .lt. therk) goto 40
6565 ! check for circles that are essentially parallel
6567 if (cc .gt. delta) goto 50
6572 ! check to see if sphere of interest is completely buried
6575 if (pix2-cc .le. td) goto 170
6581 ! find T value of circle intersections
6584 if (omit(k)) goto 110
6599 ! rotation matrix elements
6611 if (.not. omit(j)) then
6616 ! rotate spheres so K vector colinear with z-axis
6618 uxj = txj*axx + tyj*axy - tzj*axz
6619 uyj = tyj*ayy - txj*ayx
6620 uzj = txj*azx + tyj*azy + tzj*azz
6621 cosine = min(1.0d0,max(-1.0d0,uzj/b(j)))
6622 if (acos(cosine) .lt. therk+ther(j)) then
6623 dsqj = uxj*uxj + uyj*uyj
6628 tr2 = risqk*dsqj - tb*tb
6634 ! get T values of intersection for K circle
6637 tb = min(1.0d0,max(-1.0d0,tb))
6639 if (tyb-txr .lt. 0.0d0) tk1 = pix2 - tk1
6641 tb = min(1.0d0,max(-1.0d0,tb))
6643 if (tyb+txr .lt. 0.0d0) tk2 = pix2 - tk2
6644 thec = (rrsq*uzj-gk*bg(j)) / (rik*ri(j)*b(j))
6645 if (abs(thec) .lt. 1.0d0) then
6647 else if (thec .ge. 1.0d0) then
6649 else if (thec .le. -1.0d0) then
6653 ! see if "tk1" is entry or exit point; check t=0 point;
6654 ! "ti" is exit point, "tf" is entry point
6656 cosine = min(1.0d0,max(-1.0d0, &
6657 (uzj*gk-uxj*rik)/(b(j)*rr)))
6658 if ((acos(cosine)-ther(j))*(tk2-tk1) .le. 0.0d0) then
6666 if (narc .ge. maxarc) then
6668 70 format (/,' SURFATOM -- Increase the Value',&
6672 if (tf .le. ti) then
6693 ! special case; K circle without intersections
6695 if (narc .le. 0) goto 90
6697 ! general case; sum up arclength and set connectivity code
6699 call sort2 (narc,arci,key)
6704 if (narc .gt. 1) then
6707 if (t .lt. arci(j)) then
6708 arcsum = arcsum + arci(j) - t
6709 exang = exang + ex(ni)
6711 if (jb .ge. maxarc) then
6713 80 format (/,' SURFATOM -- Increase the Value',&
6718 kent(jb) = maxarc*i + k
6720 kout(jb) = maxarc*k + i
6729 arcsum = arcsum + pix2 - t
6731 exang = exang + ex(ni)
6734 kent(jb) = maxarc*i + k
6736 kout(jb) = maxarc*k + i
6743 arclen = arclen + gr(k)*arcsum
6746 if (arclen .eq. 0.0d0) goto 170
6747 if (jb .eq. 0) goto 150
6749 ! find number of independent boundaries and check connectivity
6753 if (kout(k) .ne. 0) then
6760 if (m .eq. kent(ii)) then
6763 if (j .eq. jb) goto 150
6775 ! attempt to fix connectivity error by moving atom slightly
6779 140 format (/,' SURFATOM -- Connectivity Error at Atom',i6)
6788 ! compute the exposed surface area for the sphere of interest
6791 area = ib*pix2 + exang + arclen
6792 area = mod(area,4.0d0*pi) * rrsq
6794 ! attempt to fix negative area by moving atom slightly
6796 if (area .lt. 0.0d0) then
6799 160 format (/,' SURFATOM -- Negative Area at Atom',i6)
6810 end subroutine surfatom
6811 !----------------------------------------------------------------
6812 !----------------------------------------------------------------
6813 subroutine alloc_MD_arrays
6814 !EL Allocation of arrays used by MD module
6816 integer :: nres2,nres6
6819 !----------------------
6823 allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
6824 allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
6825 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
6826 allocate(fricvec(nres2,nres2))
6827 allocate(pfric_mat(nres2,nres2),vfric_mat(nres2,nres2))
6828 allocate(afric_mat(nres2,nres2),prand_mat(nres2,nres2))
6829 allocate(vrand_mat1(nres2,nres2),vrand_mat2(nres2,nres2)) !(MAXRES2,MAXRES2)
6830 allocate(pfric0_mat(nres2,nres2,0:maxflag_stoch))
6831 allocate(afric0_mat(nres2,nres2,0:maxflag_stoch))
6832 allocate(vfric0_mat(nres2,nres2,0:maxflag_stoch))
6833 allocate(prand0_mat(nres2,nres2,0:maxflag_stoch))
6834 allocate(vrand0_mat1(nres2,nres2,0:maxflag_stoch))
6835 allocate(vrand0_mat2(nres2,nres2,0:maxflag_stoch)) !(MAXRES2,MAXRES2,0:maxflag_stoch)
6836 allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
6838 allocate(mt1(nres2,nres2),mt2(nres2,nres2),mt3(nres2,nres2)) !(maxres2,maxres2)
6839 !----------------------
6841 ! commom.langevin.lang0
6843 allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
6844 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
6845 allocate(fricvec(nres2,nres2)) !(MAXRES2,MAXRES2)
6846 allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
6847 allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
6850 if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2))
6851 !----------------------
6852 ! commom.hairpin in CSA module
6853 !----------------------
6854 ! common.mce in MCM_MD module
6855 !----------------------
6857 ! common /mdgrad/ in module.energy
6858 ! common /back_constr/ in module.energy
6859 ! common /qmeas/ in module.energy
6862 allocate(potEcomp(0:n_ene+4)) !(0:n_ene+4)
6864 allocate(d_t(3,0:nres2),d_a(3,0:nres2),d_t_old(3,0:nres2)) !(3,0:MAXRES2)
6865 allocate(d_a_work(nres6)) !(6*MAXRES)
6867 allocate(DM(nres2),DU1(nres2),DU2(nres2))
6868 allocate(DMorig(nres2),DU1orig(nres2),DU2orig(nres2))
6870 write (iout,*) "Before A Allocation",nfgtasks-1
6872 allocate(Gmat(nres2,nres2),A(nres2,nres2))
6873 if(.not.allocated(Ginv)) allocate(Ginv(nres2,nres2)) !in control: ergastulum
6874 allocate(Gsqrp(nres2,nres2),Gsqrm(nres2,nres2),Gvec(nres2,nres2)) !(maxres2,maxres2)
6876 allocate(Geigen(nres2)) !(maxres2)
6877 if(.not.allocated(vtot)) allocate(vtot(nres2)) !(maxres2)
6878 ! common /inertia/ in io_conf: parmread
6879 ! real(kind=8),dimension(:),allocatable :: ISC,msc !(ntyp+1)
6880 ! common /langevin/in io read_MDpar
6881 ! real(kind=8),dimension(:),allocatable :: gamsc !(ntyp1)
6882 ! real(kind=8),dimension(:),allocatable :: stdfsc !(ntyp)
6883 ! in io_conf: parmread
6884 ! real(kind=8),dimension(:),allocatable :: restok !(ntyp+1)
6885 ! common /mdpmpi/ in control: ergastulum
6886 if(.not.allocated(ng_start)) allocate(ng_start(0:nfgtasks-1))
6887 if(.not.allocated(ng_counts)) allocate(ng_counts(0:nfgtasks-1))
6888 if(.not.allocated(nginv_counts)) allocate(nginv_counts(0:nfgtasks-1)) !(0:MaxProcs-1)
6889 if(.not.allocated(nginv_start)) allocate(nginv_start(0:nfgtasks)) !(0:MaxProcs)
6890 !----------------------
6891 ! common.muca in read_muca
6892 ! common /double_muca/
6893 ! real(kind=8) :: elow,ehigh,factor,hbin,factor_min
6894 ! real(kind=8),dimension(:),allocatable :: emuca,nemuca,&
6895 ! nemuca2,hist !(4*maxres)
6896 ! common /integer_muca/
6897 ! integer :: nmuca,imtime,muca_smooth
6899 ! real(kind=8),dimension(:),allocatable :: elowi,ehighi !(maxprocs)
6900 !----------------------
6902 ! common /mdgrad/ in module.energy
6903 ! common /back_constr/ in module.energy
6904 ! common /qmeas/ in module.energy
6908 allocate(d_t_work(nres6),d_t_work_new(nres6),d_af_work(nres6))
6909 allocate(d_as_work(nres6),kinetic_force(nres6)) !(MAXRES6)
6910 allocate(d_t_new(3,0:nres2),d_a_old(3,0:nres2),d_a_short(3,0:nres2)) !,d_a !(3,0:MAXRES2)
6911 allocate(stdforcp(nres),stdforcsc(nres)) !(MAXRES)
6912 !----------------------
6914 allocate(D_ban(nres6)) !(MAXRES6) maxres6=6*maxres
6915 ! common /stochcalc/ stochforcvec
6916 allocate(stochforcvec(nres6)) !(MAXRES6) maxres6=6*maxres
6917 !----------------------
6919 end subroutine alloc_MD_arrays
6920 !-----------------------------------------------------------------------------
6921 !-----------------------------------------------------------------------------