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(20)
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 !-----------------------------------------------------------------------------
648 !-----------------------------------------------------------------------------
649 subroutine kinetic(KE_total)
650 !----------------------------------------------------------------
651 ! This subroutine calculates the total kinetic energy of the chain
652 !-----------------------------------------------------------------
654 ! implicit real*8 (a-h,o-z)
655 ! include 'DIMENSIONS'
656 ! include 'COMMON.VAR'
657 ! include 'COMMON.CHAIN'
658 ! include 'COMMON.DERIV'
659 ! include 'COMMON.GEO'
660 ! include 'COMMON.LOCAL'
661 ! include 'COMMON.INTERACT'
662 ! include 'COMMON.MD'
663 ! include 'COMMON.IOUNITS'
664 real(kind=8) :: KE_total,mscab
666 integer :: i,j,k,iti,mnum,term
667 real(kind=8) :: KEt_p,KEt_sc,KEr_p,KEr_sc,incr(3),&
670 write (iout,*) "Velocities, kietic"
672 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
673 (d_t(j,i+nres),j=1,3)
678 ! write (iout,*) "ISC",(isc(itype(i,1)),i=1,nres)
679 ! The translational part for peptide virtual bonds
684 ! if (molnum(nct).gt.3) term=nct
687 if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
688 ! write (iout,*) "Kinetic trp:",i,(incr(j),j=1,3),mp(mnum)
691 v(j)=incr(j)+d_t(j,i)
695 v(j)=incr(j)+0.5d0*d_t(j,i)
698 vtot(i)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
699 KEt_p=KEt_p+mp(mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
701 incr(j)=incr(j)+d_t(j,i)
704 ! write(iout,*) 'KEt_p', KEt_p
705 ! The translational part for the side chain virtual bond
706 ! Only now we can initialize incr with zeros. It must be equal
707 ! to the velocities of the first Calpha.
713 iti=iabs(itype(i,mnum))
714 ! if (mnum.ge.4) then
719 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
720 if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
721 .or.(mnum.ge.4)) then
727 v(j)=incr(j)+d_t(j,nres+i)
730 ! write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3)
731 ! write (iout,*) "i",i," msc",msc(iti,mnum)," v",(v(j),j=1,3)
732 KEt_sc=KEt_sc+mscab*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
733 vtot(i+nres)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
735 incr(j)=incr(j)+d_t(j,i)
739 ! write(iout,*) 'KEt_sc', KEt_sc
740 ! The part due to stretching and rotation of the peptide groups
744 ! write (iout,*) "i",i
745 ! write (iout,*) "i",i," mag1",mag1," mag2",mag2
749 ! write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3)
750 KEr_p=KEr_p+Ip(mnum)*(incr(1)*incr(1)+incr(2)*incr(2) &
754 ! write(iout,*) 'KEr_p', KEr_p
755 ! The rotational part of the side chain virtual bond
759 iti=iabs(itype(i,mnum))
760 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
761 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
762 .and.(mnum.lt.4)) then
764 incr(j)=d_t(j,nres+i)
766 ! write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3)
767 KEr_sc=KEr_sc+Isc(iti,mnum)*(incr(1)*incr(1)+incr(2)*incr(2)+ &
771 ! The total kinetic energy
773 ! write(iout,*) 'KEr_sc', KEr_sc
774 KE_total=0.5d0*(KEt_p+KEt_sc+0.25d0*KEr_p+KEr_sc)
775 ! write (iout,*) "KE_total",KE_total
777 end subroutine kinetic
778 !-----------------------------------------------------------------------------
780 !-----------------------------------------------------------------------------
782 !------------------------------------------------
783 ! The driver for molecular dynamics subroutines
784 !------------------------------------------------
787 use control, only:tcpu,ovrtim
788 ! use io_comm, only:ilen
790 use compare, only:secondary2,hairpin
791 use io, only:cartout,statout
792 ! implicit real*8 (a-h,o-z)
793 ! include 'DIMENSIONS'
796 integer :: IERROR,ERRCODE
798 ! include 'COMMON.SETUP'
799 ! include 'COMMON.CONTROL'
800 ! include 'COMMON.VAR'
801 ! include 'COMMON.MD'
803 ! include 'COMMON.LANGEVIN'
805 ! include 'COMMON.LANGEVIN.lang0'
807 ! include 'COMMON.CHAIN'
808 ! include 'COMMON.DERIV'
809 ! include 'COMMON.GEO'
810 ! include 'COMMON.LOCAL'
811 ! include 'COMMON.INTERACT'
812 ! include 'COMMON.IOUNITS'
813 ! include 'COMMON.NAMES'
814 ! include 'COMMON.TIME1'
815 ! include 'COMMON.HAIRPIN'
816 real(kind=8),dimension(3) :: L,vcm
818 real(kind=8),dimension(6*nres) :: v_work,v_transf !(maxres6) maxres6=6*maxres
820 integer :: rstcount !ilen,
822 character(len=50) :: tytul
823 !el common /gucio/ cm
825 integer,dimension(4,nres/3) :: iharp !(4,nres/3)(4,maxres/3)
827 real(kind=8) :: tt0,scalfac
828 integer :: nres2,itime
833 print *,"MY tmpdir",tmpdir,ilen(tmpdir)
834 if (ilen(tmpdir).gt.0) &
835 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_" &
836 //liczba(:ilen(liczba))//'.rst')
838 if (ilen(tmpdir).gt.0) &
839 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"//'.rst')
846 write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
852 print *,"just befor setup matix",nres
853 ! Determine the inverse of the inertia matrix.
854 call setup_MD_matrices
856 print *,"AFTER SETUP MATRICES"
858 print *,"AFTER INIT MD"
861 t_MDsetup = MPI_Wtime()-tt0
863 t_MDsetup = tcpu()-tt0
866 ! Entering the MD loop
872 if (lang.eq.2 .or. lang.eq.3) then
876 call sd_verlet_p_setup
878 call sd_verlet_ciccotti_setup
882 pfric0_mat(i,j,0)=pfric_mat(i,j)
883 afric0_mat(i,j,0)=afric_mat(i,j)
884 vfric0_mat(i,j,0)=vfric_mat(i,j)
885 prand0_mat(i,j,0)=prand_mat(i,j)
886 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
887 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
892 flag_stoch(i)=.false.
896 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
898 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
902 else if (lang.eq.1 .or. lang.eq.4) then
903 print *,"before setup_fricmat"
905 print *,"after setup_fricmat"
908 t_langsetup=MPI_Wtime()-tt0
911 t_langsetup=tcpu()-tt0
914 do itime=1,n_timestep
916 if (large.and. mod(itime,ntwe).eq.0) &
917 write (iout,*) "itime",itime
919 if (lang.gt.0 .and. surfarea .and. &
920 mod(itime,reset_fricmat).eq.0) then
921 if (lang.eq.2 .or. lang.eq.3) then
925 call sd_verlet_p_setup
927 call sd_verlet_ciccotti_setup
931 pfric0_mat(i,j,0)=pfric_mat(i,j)
932 afric0_mat(i,j,0)=afric_mat(i,j)
933 vfric0_mat(i,j,0)=vfric_mat(i,j)
934 prand0_mat(i,j,0)=prand_mat(i,j)
935 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
936 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
941 flag_stoch(i)=.false.
944 else if (lang.eq.1 .or. lang.eq.4) then
947 write (iout,'(a,i10)') &
948 "Friction matrix reset based on surface area, itime",itime
950 if (reset_vel .and. tbf .and. lang.eq.0 &
951 .and. mod(itime,count_reset_vel).eq.0) then
953 write(iout,'(a,f20.2)') &
954 "Velocities reset to random values, time",totT
957 d_t_old(j,i)=d_t(j,i)
961 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
965 d_t(j,0)=d_t(j,0)-vcm(j)
968 kinetic_T=2.0d0/(dimen3*Rb)*EK
969 scalfac=dsqrt(T_bath/kinetic_T)
970 write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
973 d_t_old(j,i)=scalfac*d_t(j,i)
979 ! Time-reversible RESPA algorithm
980 ! (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
981 call RESPA_step(itime)
983 ! Variable time step algorithm.
984 call velverlet_step(itime)
988 call brown_step(itime)
990 print *,"Brown dynamics not here!"
992 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
999 if (mod(itime,ntwe).eq.0) then
1003 ! call check_ecartint
1013 v_work(ind)=d_t(j,i)
1018 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1.and.mnum.lt.4) then
1021 v_work(ind)=d_t(j,i+nres)
1026 write (66,'(80f10.5)') &
1027 ((d_t(j,i),j=1,3),i=0,nres-1),((d_t(j,i+nres),j=1,3),i=1,nres)
1031 v_transf(i)=v_transf(i)+gvec(j,i)*v_work(j)
1033 v_transf(i)= v_transf(i)*dsqrt(geigen(i))
1035 write (67,'(80f10.5)') (v_transf(i),i=1,ind)
1038 if (mod(itime,ntwx).eq.0) then
1040 write (tytul,'("time",f8.2)') totT
1042 call hairpin(.true.,nharp,iharp)
1043 call secondary2(.true.)
1044 call pdbout(potE,tytul,ipdb)
1049 if (rstcount.eq.1000.or.itime.eq.n_timestep) then
1050 open(irest2,file=rest2name,status='unknown')
1051 write(irest2,*) totT,EK,potE,totE,t_bath
1053 ! AL 4/17/17: Now writing d_t(0,:) too
1055 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
1057 ! AL 4/17/17: Now writing d_c(0,:) too
1059 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
1067 t_MD=MPI_Wtime()-tt0
1071 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))') &
1073 'MD calculations setup:',t_MDsetup,&
1074 'Energy & gradient evaluation:',t_enegrad,&
1075 'Stochastic MD setup:',t_langsetup,&
1076 'Stochastic MD step setup:',t_sdsetup,&
1078 write (iout,'(/28(1h=),a25,27(1h=))') &
1079 ' End of MD calculation '
1081 write (iout,*) "time for etotal",t_etotal," elong",t_elong,&
1083 write (iout,*) "time_fric",time_fric," time_stoch",time_stoch,&
1084 " time_fricmatmult",time_fricmatmult," time_fsample ",&
1089 !-----------------------------------------------------------------------------
1090 subroutine velverlet_step(itime)
1091 !-------------------------------------------------------------------------------
1092 ! Perform a single velocity Verlet step; the time step can be rescaled if
1093 ! increments in accelerations exceed the threshold
1094 !-------------------------------------------------------------------------------
1095 ! implicit real*8 (a-h,o-z)
1096 ! include 'DIMENSIONS'
1098 use control, only:tcpu
1100 use minimm, only:minim_dc
1103 integer :: ierror,ierrcode
1104 real(kind=8) :: errcode
1106 ! include 'COMMON.SETUP'
1107 ! include 'COMMON.VAR'
1108 ! include 'COMMON.MD'
1110 ! include 'COMMON.LANGEVIN'
1112 ! include 'COMMON.LANGEVIN.lang0'
1114 ! include 'COMMON.CHAIN'
1115 ! include 'COMMON.DERIV'
1116 ! include 'COMMON.GEO'
1117 ! include 'COMMON.LOCAL'
1118 ! include 'COMMON.INTERACT'
1119 ! include 'COMMON.IOUNITS'
1120 ! include 'COMMON.NAMES'
1121 ! include 'COMMON.TIME1'
1122 ! include 'COMMON.MUCA'
1123 real(kind=8),dimension(3) :: vcm,incr
1124 real(kind=8),dimension(3) :: L
1125 integer :: count,rstcount !ilen,
1127 character(len=50) :: tytul
1128 integer :: maxcount_scale = 30
1129 !el common /gucio/ cm
1130 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1131 !el common /stochcalc/ stochforcvec
1132 integer :: icount_scale,itime_scal,i,j,ifac_time,iretcode,itime
1134 real(kind=8) :: epdrift,tt0,fac_time
1136 if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres)) !(MAXRES6) maxres6=6*maxres
1142 else if (lang.eq.2 .or. lang.eq.3) then
1144 call stochastic_force(stochforcvec)
1147 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
1149 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1156 icount_scale=icount_scale+1
1157 ! write(iout,*) "icount_scale",icount_scale,scalel
1158 if (icount_scale.gt.maxcount_scale) then
1160 "ERROR: too many attempts at scaling down the time step. ",&
1161 "amax=",amax,"epdrift=",epdrift,&
1162 "damax=",damax,"edriftmax=",edriftmax,&
1166 call MPI_Abort(MPI_COMM_WORLD,IERROR,IERRCODE)
1170 ! First step of the velocity Verlet algorithm
1175 else if (lang.eq.3) then
1177 call sd_verlet1_ciccotti
1179 else if (lang.eq.1) then
1184 ! Build the chain from the newly calculated coordinates
1185 call chainbuild_cart
1186 if (rattle) call rattle1
1188 if (large.and. mod(itime,ntwe).eq.0) then
1189 write (iout,*) "Cartesian and internal coordinates: step 1"
1194 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),&
1195 (dc(j,i+nres),j=1,3)
1197 write (iout,*) "Accelerations"
1199 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1200 (d_a(j,i+nres),j=1,3)
1202 write (iout,*) "Velocities, step 1"
1204 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1205 (d_t(j,i+nres),j=1,3)
1214 ! Calculate energy and forces
1216 call etotal(potEcomp)
1217 ! AL 4/17/17: Reduce the steps if NaNs occurred.
1218 if (potEcomp(0).gt.0.99e18 .or. isnan(potEcomp(0)).gt.0) then
1219 call enerprint(potEcomp)
1221 if (icount_scale.gt.15) then
1222 write (iout,*) "Tu jest problem",potEcomp(0),d_time
1223 ! call gen_rand_conf(1,*335)
1224 ! call minim_dc(potEcomp(0),iretcode,100)
1227 ! call etotal(potEcomp)
1228 ! write(iout,*) "needed to repara,",potEcomp
1231 ! 335 write(iout,*) "Failed genrand"
1235 if (large.and. mod(itime,ntwe).eq.0) &
1236 call enerprint(potEcomp)
1239 t_etotal=t_etotal+MPI_Wtime()-tt0
1241 t_etotal=t_etotal+tcpu()-tt0
1244 potE=potEcomp(0)-potEcomp(20)
1246 ! Get the new accelerations
1249 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1251 t_enegrad=t_enegrad+tcpu()-tt0
1253 ! Determine maximum acceleration and scale down the timestep if needed
1255 amax=amax/(itime_scal+1)**2
1256 call predict_edrift(epdrift)
1257 ! write(iout,*) "amax=",amax,damax,epdrift,edriftmax,amax/(itime_scal+1)
1259 ! write (iout,*) "before enter if",scalel,icount_scale
1260 if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then
1261 ! write(iout,*) "I enter if"
1262 ! Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
1264 ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax)) &
1266 itime_scal=itime_scal+ifac_time
1267 ! fac_time=dmin1(damax/amax,0.5d0)
1268 fac_time=0.5d0**ifac_time
1269 d_time=d_time*fac_time
1270 if (lang.eq.2 .or. lang.eq.3) then
1272 ! write (iout,*) "Calling sd_verlet_setup: 1"
1273 ! Rescale the stochastic forces and recalculate or restore
1274 ! the matrices of tinker integrator
1275 if (itime_scal.gt.maxflag_stoch) then
1276 if (large) write (iout,'(a,i5,a)') &
1277 "Calculate matrices for stochastic step;",&
1278 " itime_scal ",itime_scal
1280 call sd_verlet_p_setup
1282 call sd_verlet_ciccotti_setup
1284 write (iout,'(2a,i3,a,i3,1h.)') &
1285 "Warning: cannot store matrices for stochastic",&
1286 " integration because the index",itime_scal,&
1287 " is greater than",maxflag_stoch
1288 write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",&
1289 " integration Langevin algorithm for better efficiency."
1290 else if (flag_stoch(itime_scal)) then
1291 if (large) write (iout,'(a,i5,a,l1)') &
1292 "Restore matrices for stochastic step; itime_scal ",&
1293 itime_scal," flag ",flag_stoch(itime_scal)
1296 pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
1297 afric_mat(i,j)=afric0_mat(i,j,itime_scal)
1298 vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
1299 prand_mat(i,j)=prand0_mat(i,j,itime_scal)
1300 vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
1301 vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
1305 if (large) write (iout,'(2a,i5,a,l1)') &
1306 "Calculate & store matrices for stochastic step;",&
1307 " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
1309 call sd_verlet_p_setup
1311 call sd_verlet_ciccotti_setup
1313 flag_stoch(ifac_time)=.true.
1316 pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
1317 afric0_mat(i,j,itime_scal)=afric_mat(i,j)
1318 vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
1319 prand0_mat(i,j,itime_scal)=prand_mat(i,j)
1320 vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
1321 vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
1325 fac_time=1.0d0/dsqrt(fac_time)
1327 stochforcvec(i)=fac_time*stochforcvec(i)
1330 else if (lang.eq.1) then
1331 ! Rescale the accelerations due to stochastic forces
1332 fac_time=1.0d0/dsqrt(fac_time)
1334 d_as_work(i)=d_as_work(i)*fac_time
1337 if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)') &
1338 "itime",itime," Timestep scaled down to ",&
1339 d_time," ifac_time",ifac_time," itime_scal",itime_scal
1341 ! Second step of the velocity Verlet algorithm
1346 else if (lang.eq.3) then
1348 call sd_verlet2_ciccotti
1350 else if (lang.eq.1) then
1355 if (rattle) call rattle2
1358 if (d_time.ne.d_time0) then
1361 if (lang.eq.2 .or. lang.eq.3) then
1362 if (large) write (iout,'(a)') &
1363 "Restore original matrices for stochastic step"
1364 ! write (iout,*) "Calling sd_verlet_setup: 2"
1365 ! Restore the matrices of tinker integrator if the time step has been restored
1368 pfric_mat(i,j)=pfric0_mat(i,j,0)
1369 afric_mat(i,j)=afric0_mat(i,j,0)
1370 vfric_mat(i,j)=vfric0_mat(i,j,0)
1371 prand_mat(i,j)=prand0_mat(i,j,0)
1372 vrand_mat1(i,j)=vrand0_mat1(i,j,0)
1373 vrand_mat2(i,j)=vrand0_mat2(i,j,0)
1381 ! Calculate the kinetic and the total energy and the kinetic temperature
1385 ! call kinetic1(EK1)
1386 ! write (iout,*) "step",itime," EK",EK," EK1",EK1
1388 ! Couple the system to Berendsen bath if needed
1389 if (tbf .and. lang.eq.0) then
1392 kinetic_T=2.0d0/(dimen3*Rb)*EK
1393 ! Backup the coordinates, velocities, and accelerations
1397 d_t_old(j,i)=d_t(j,i)
1398 d_a_old(j,i)=d_a(j,i)
1402 if (mod(itime,ntwe).eq.0 .and. large) then
1403 write (iout,*) "Velocities, step 2"
1405 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1406 (d_t(j,i+nres),j=1,3)
1411 end subroutine velverlet_step
1412 !-----------------------------------------------------------------------------
1413 subroutine RESPA_step(itime)
1414 !-------------------------------------------------------------------------------
1415 ! Perform a single RESPA step.
1416 !-------------------------------------------------------------------------------
1417 ! implicit real*8 (a-h,o-z)
1418 ! include 'DIMENSIONS'
1422 use control, only:tcpu
1424 ! use io_conf, only:cartprint
1427 integer :: IERROR,ERRCODE
1429 ! include 'COMMON.SETUP'
1430 ! include 'COMMON.CONTROL'
1431 ! include 'COMMON.VAR'
1432 ! include 'COMMON.MD'
1434 ! include 'COMMON.LANGEVIN'
1436 ! include 'COMMON.LANGEVIN.lang0'
1438 ! include 'COMMON.CHAIN'
1439 ! include 'COMMON.DERIV'
1440 ! include 'COMMON.GEO'
1441 ! include 'COMMON.LOCAL'
1442 ! include 'COMMON.INTERACT'
1443 ! include 'COMMON.IOUNITS'
1444 ! include 'COMMON.NAMES'
1445 ! include 'COMMON.TIME1'
1446 real(kind=8),dimension(0:n_ene) :: energia_short,energia_long
1447 real(kind=8),dimension(3) :: L,vcm,incr
1448 real(kind=8),dimension(3,0:2*nres) :: dc_old0,d_t_old0,d_a_old0 !(3,0:maxres2) maxres2=2*maxres
1449 logical :: PRINT_AMTS_MSG = .false.
1450 integer :: count,rstcount !ilen,
1452 character(len=50) :: tytul
1453 integer :: maxcount_scale = 10
1454 !el common /gucio/ cm
1455 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1456 !el common /stochcalc/ stochforcvec
1457 integer :: itt,i,j,itsplit,itime
1459 !el common /cipiszcze/ itt
1461 real(kind=8) :: epdrift,tt0,epdriftmax
1464 if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres)) !(MAXRES6) maxres6=6*maxres
1468 if (large.and. mod(itime,ntwe).eq.0) then
1469 write (iout,*) "***************** RESPA itime",itime
1470 write (iout,*) "Cartesian and internal coordinates: step 0"
1472 call pdbout(0.0d0,"cipiszcze",iout)
1474 write (iout,*) "Accelerations from long-range forces"
1476 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1477 (d_a(j,i+nres),j=1,3)
1479 write (iout,*) "Velocities, step 0"
1481 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1482 (d_t(j,i+nres),j=1,3)
1487 ! Perform the initial RESPA step (increment velocities)
1488 ! write (iout,*) "*********************** RESPA ini"
1491 if (mod(itime,ntwe).eq.0 .and. large) then
1492 write (iout,*) "Velocities, end"
1494 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1495 (d_t(j,i+nres),j=1,3)
1499 ! Compute the short-range forces
1505 ! 7/2/2009 commented out
1507 ! call etotal_short(energia_short)
1510 ! 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
1513 d_a(j,i)=d_a_short(j,i)
1517 if (large.and. mod(itime,ntwe).eq.0) then
1518 write (iout,*) "energia_short",energia_short(0)
1519 write (iout,*) "Accelerations from short-range forces"
1521 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1522 (d_a(j,i+nres),j=1,3)
1527 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1529 t_enegrad=t_enegrad+tcpu()-tt0
1534 d_t_old(j,i)=d_t(j,i)
1535 d_a_old(j,i)=d_a(j,i)
1538 ! 6/30/08 A-MTS: attempt at increasing the split number
1541 dc_old0(j,i)=dc_old(j,i)
1542 d_t_old0(j,i)=d_t_old(j,i)
1543 d_a_old0(j,i)=d_a_old(j,i)
1546 if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
1547 if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
1554 ! write (iout,*) "itime",itime," ntime_split",ntime_split
1555 ! Split the time step
1556 d_time=d_time0/ntime_split
1557 ! Perform the short-range RESPA steps (velocity Verlet increments of
1558 ! positions and velocities using short-range forces)
1559 ! write (iout,*) "*********************** RESPA split"
1560 do itsplit=1,ntime_split
1563 else if (lang.eq.2 .or. lang.eq.3) then
1565 call stochastic_force(stochforcvec)
1568 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
1570 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1575 ! First step of the velocity Verlet algorithm
1580 else if (lang.eq.3) then
1582 call sd_verlet1_ciccotti
1584 else if (lang.eq.1) then
1589 ! Build the chain from the newly calculated coordinates
1590 call chainbuild_cart
1591 if (rattle) call rattle1
1593 if (large.and. mod(itime,ntwe).eq.0) then
1594 write (iout,*) "***** ITSPLIT",itsplit
1595 write (iout,*) "Cartesian and internal coordinates: step 1"
1596 call pdbout(0.0d0,"cipiszcze",iout)
1599 write (iout,*) "Velocities, step 1"
1601 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1602 (d_t(j,i+nres),j=1,3)
1611 ! Calculate energy and forces
1613 call etotal_short(energia_short)
1614 ! AL 4/17/17: Exit itime_split loop when energy goes infinite
1615 if (energia_short(0).gt.0.99e20 .or. isnan(energia_short(0)) ) then
1616 if (PRINT_AMTS_MSG) &
1617 write (iout,*) "Infinities/NaNs in energia_short",energia_short(0),"; increasing ntime_split to",ntime_split
1618 ntime_split=ntime_split*2
1619 if (ntime_split.gt.maxtime_split) then
1622 "Cannot rescue the run; aborting job. Retry with a smaller time step"
1624 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1627 "Cannot rescue the run; terminating. Retry with a smaller time step"
1633 if (large.and. mod(itime,ntwe).eq.0) &
1634 call enerprint(energia_short)
1637 t_eshort=t_eshort+MPI_Wtime()-tt0
1639 t_eshort=t_eshort+tcpu()-tt0
1643 ! Get the new accelerations
1645 ! 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
1648 d_a_short(j,i)=d_a(j,i)
1652 if (large.and. mod(itime,ntwe).eq.0) then
1653 write (iout,*)"energia_short",energia_short(0)
1654 write (iout,*) "Accelerations from short-range forces"
1656 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1657 (d_a(j,i+nres),j=1,3)
1662 ! Determine maximum acceleration and scale down the timestep if needed
1664 amax=amax/ntime_split**2
1665 call predict_edrift(epdrift)
1666 if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0) &
1667 write (iout,*) "amax",amax," damax",damax,&
1668 " epdrift",epdrift," epdriftmax",epdriftmax
1669 ! Exit loop and try with increased split number if the change of
1670 ! acceleration is too big
1671 if (amax.gt.damax .or. epdrift.gt.edriftmax) then
1672 if (ntime_split.lt.maxtime_split) then
1674 ntime_split=ntime_split*2
1675 ! AL 4/17/17: We should exit the itime_split loop when acceleration change is too big
1679 dc_old(j,i)=dc_old0(j,i)
1680 d_t_old(j,i)=d_t_old0(j,i)
1681 d_a_old(j,i)=d_a_old0(j,i)
1684 if (PRINT_AMTS_MSG) then
1685 write (iout,*) "acceleration/energy drift too large",amax,&
1686 epdrift," split increased to ",ntime_split," itime",itime,&
1692 "Uh-hu. Bumpy landscape. Maximum splitting number",&
1694 " already reached!!! Trying to carry on!"
1698 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1700 t_enegrad=t_enegrad+tcpu()-tt0
1702 ! Second step of the velocity Verlet algorithm
1707 else if (lang.eq.3) then
1709 call sd_verlet2_ciccotti
1711 else if (lang.eq.1) then
1716 if (rattle) call rattle2
1717 ! Backup the coordinates, velocities, and accelerations
1721 d_t_old(j,i)=d_t(j,i)
1722 d_a_old(j,i)=d_a(j,i)
1729 ! Restore the time step
1731 ! Compute long-range forces
1738 call etotal_long(energia_long)
1739 if (energia_long(0).gt.0.99e20 .or. isnan(energia_long(0))) then
1742 "Infinitied/NaNs in energia_long, Aborting MPI job."
1744 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1746 write (iout,*) "Infinitied/NaNs in energia_long, terminating."
1750 if (large.and. mod(itime,ntwe).eq.0) &
1751 call enerprint(energia_long)
1754 t_elong=t_elong+MPI_Wtime()-tt0
1756 t_elong=t_elong+tcpu()-tt0
1762 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1764 t_enegrad=t_enegrad+tcpu()-tt0
1766 ! Compute accelerations from long-range forces
1768 if (large.and. mod(itime,ntwe).eq.0) then
1769 write (iout,*) "energia_long",energia_long(0)
1770 write (iout,*) "Cartesian and internal coordinates: step 2"
1772 call pdbout(0.0d0,"cipiszcze",iout)
1774 write (iout,*) "Accelerations from long-range forces"
1776 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1777 (d_a(j,i+nres),j=1,3)
1779 write (iout,*) "Velocities, step 2"
1781 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1782 (d_t(j,i+nres),j=1,3)
1786 ! Compute the final RESPA step (increment velocities)
1787 ! write (iout,*) "*********************** RESPA fin"
1789 ! Compute the complete potential energy
1791 potEcomp(i)=energia_short(i)+energia_long(i)
1793 potE=potEcomp(0)-potEcomp(20)
1794 ! potE=energia_short(0)+energia_long(0)
1797 ! Calculate the kinetic and the total energy and the kinetic temperature
1800 ! Couple the system to Berendsen bath if needed
1801 if (tbf .and. lang.eq.0) then
1804 kinetic_T=2.0d0/(dimen3*Rb)*EK
1805 ! Backup the coordinates, velocities, and accelerations
1807 if (mod(itime,ntwe).eq.0 .and. large) then
1808 write (iout,*) "Velocities, end"
1810 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1811 (d_t(j,i+nres),j=1,3)
1816 end subroutine RESPA_step
1817 !-----------------------------------------------------------------------------
1818 subroutine RESPA_vel
1819 ! First and last RESPA step (incrementing velocities using long-range
1822 ! implicit real*8 (a-h,o-z)
1823 ! include 'DIMENSIONS'
1824 ! include 'COMMON.CONTROL'
1825 ! include 'COMMON.VAR'
1826 ! include 'COMMON.MD'
1827 ! include 'COMMON.CHAIN'
1828 ! include 'COMMON.DERIV'
1829 ! include 'COMMON.GEO'
1830 ! include 'COMMON.LOCAL'
1831 ! include 'COMMON.INTERACT'
1832 ! include 'COMMON.IOUNITS'
1833 ! include 'COMMON.NAMES'
1834 integer :: i,j,inres,mnum
1837 d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
1841 d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
1846 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
1847 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
1848 .and.(mnum.lt.4)) then
1851 d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
1856 end subroutine RESPA_vel
1857 !-----------------------------------------------------------------------------
1859 ! Applying velocity Verlet algorithm - step 1 to coordinates
1861 ! implicit real*8 (a-h,o-z)
1862 ! include 'DIMENSIONS'
1863 ! include 'COMMON.CONTROL'
1864 ! include 'COMMON.VAR'
1865 ! include 'COMMON.MD'
1866 ! include 'COMMON.CHAIN'
1867 ! include 'COMMON.DERIV'
1868 ! include 'COMMON.GEO'
1869 ! include 'COMMON.LOCAL'
1870 ! include 'COMMON.INTERACT'
1871 ! include 'COMMON.IOUNITS'
1872 ! include 'COMMON.NAMES'
1873 real(kind=8) :: adt,adt2
1874 integer :: i,j,inres,mnum
1877 write (iout,*) "VELVERLET1 START: DC"
1879 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),&
1880 (dc(j,i+nres),j=1,3)
1884 adt=d_a_old(j,0)*d_time
1886 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1887 d_t_new(j,0)=d_t_old(j,0)+adt2
1888 d_t(j,0)=d_t_old(j,0)+adt
1892 adt=d_a_old(j,i)*d_time
1894 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1895 d_t_new(j,i)=d_t_old(j,i)+adt2
1896 d_t(j,i)=d_t_old(j,i)+adt
1901 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
1902 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
1903 .and.(mnum.lt.4)) then
1906 adt=d_a_old(j,inres)*d_time
1908 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1909 d_t_new(j,inres)=d_t_old(j,inres)+adt2
1910 d_t(j,inres)=d_t_old(j,inres)+adt
1915 write (iout,*) "VELVERLET1 END: DC"
1917 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),&
1918 (dc(j,i+nres),j=1,3)
1922 end subroutine verlet1
1923 !-----------------------------------------------------------------------------
1925 ! Step 2 of the velocity Verlet algorithm: update velocities
1927 ! implicit real*8 (a-h,o-z)
1928 ! include 'DIMENSIONS'
1929 ! include 'COMMON.CONTROL'
1930 ! include 'COMMON.VAR'
1931 ! include 'COMMON.MD'
1932 ! include 'COMMON.CHAIN'
1933 ! include 'COMMON.DERIV'
1934 ! include 'COMMON.GEO'
1935 ! include 'COMMON.LOCAL'
1936 ! include 'COMMON.INTERACT'
1937 ! include 'COMMON.IOUNITS'
1938 ! include 'COMMON.NAMES'
1939 integer :: i,j,inres,mnum
1942 d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
1946 d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
1951 ! iti=iabs(itype(i,mnum))
1952 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
1953 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
1954 .and.(mnum.lt.4)) then
1957 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
1962 end subroutine verlet2
1963 !-----------------------------------------------------------------------------
1964 subroutine sddir_precalc
1965 ! Applying velocity Verlet algorithm - step 1 to coordinates
1966 ! implicit real*8 (a-h,o-z)
1967 ! include 'DIMENSIONS'
1973 ! include 'COMMON.CONTROL'
1974 ! include 'COMMON.VAR'
1975 ! include 'COMMON.MD'
1977 ! include 'COMMON.LANGEVIN'
1979 ! include 'COMMON.LANGEVIN.lang0'
1981 ! include 'COMMON.CHAIN'
1982 ! include 'COMMON.DERIV'
1983 ! include 'COMMON.GEO'
1984 ! include 'COMMON.LOCAL'
1985 ! include 'COMMON.INTERACT'
1986 ! include 'COMMON.IOUNITS'
1987 ! include 'COMMON.NAMES'
1988 ! include 'COMMON.TIME1'
1989 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1990 !el common /stochcalc/ stochforcvec
1991 real(kind=8) :: time00
1993 ! Compute friction and stochastic forces
1998 time_fric=time_fric+MPI_Wtime()-time00
2000 call stochastic_force(stochforcvec)
2001 time_stoch=time_stoch+MPI_Wtime()-time00
2004 ! Compute the acceleration due to friction forces (d_af_work) and stochastic
2005 ! forces (d_as_work)
2007 call ginv_mult(fric_work, d_af_work)
2008 call ginv_mult(stochforcvec, d_as_work)
2010 end subroutine sddir_precalc
2011 !-----------------------------------------------------------------------------
2012 subroutine sddir_verlet1
2013 ! Applying velocity Verlet algorithm - step 1 to velocities
2016 ! implicit real*8 (a-h,o-z)
2017 ! include 'DIMENSIONS'
2018 ! include 'COMMON.CONTROL'
2019 ! include 'COMMON.VAR'
2020 ! include 'COMMON.MD'
2022 ! include 'COMMON.LANGEVIN'
2024 ! include 'COMMON.LANGEVIN.lang0'
2026 ! include 'COMMON.CHAIN'
2027 ! include 'COMMON.DERIV'
2028 ! include 'COMMON.GEO'
2029 ! include 'COMMON.LOCAL'
2030 ! include 'COMMON.INTERACT'
2031 ! include 'COMMON.IOUNITS'
2032 ! include 'COMMON.NAMES'
2033 ! Revised 3/31/05 AL: correlation between random contributions to
2034 ! position and velocity increments included.
2035 real(kind=8) :: sqrt13 = 0.57735026918962576451d0 ! 1/sqrt(3)
2036 real(kind=8) :: adt,adt2
2037 integer :: i,j,ind,inres,mnum
2039 ! Add the contribution from BOTH friction and stochastic force to the
2040 ! coordinates, but ONLY the contribution from the friction forces to velocities
2043 adt=(d_a_old(j,0)+d_af_work(j))*d_time
2044 adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
2045 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
2046 d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
2047 d_t(j,0)=d_t_old(j,0)+adt
2052 adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
2053 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
2054 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
2055 d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
2056 d_t(j,i)=d_t_old(j,i)+adt
2062 ! iti=iabs(itype(i,mnum))
2063 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2064 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2065 .and.(mnum.lt.4)) then
2068 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
2069 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
2070 ! write(iout,*) "adt",adt,"ads",adt2
2071 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
2072 d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
2073 d_t(j,inres)=d_t_old(j,inres)+adt
2079 end subroutine sddir_verlet1
2080 !-----------------------------------------------------------------------------
2081 subroutine sddir_verlet2
2082 ! Calculating the adjusted velocities for accelerations
2085 ! implicit real*8 (a-h,o-z)
2086 ! include 'DIMENSIONS'
2087 ! include 'COMMON.CONTROL'
2088 ! include 'COMMON.VAR'
2089 ! include 'COMMON.MD'
2091 ! include 'COMMON.LANGEVIN'
2093 ! include 'COMMON.LANGEVIN.lang0'
2095 ! include 'COMMON.CHAIN'
2096 ! include 'COMMON.DERIV'
2097 ! include 'COMMON.GEO'
2098 ! include 'COMMON.LOCAL'
2099 ! include 'COMMON.INTERACT'
2100 ! include 'COMMON.IOUNITS'
2101 ! include 'COMMON.NAMES'
2102 real(kind=8),dimension(6*nres) :: stochforcvec,d_as_work1 !(MAXRES6) maxres6=6*maxres
2103 real(kind=8) :: cos60 = 0.5d0, sin60 = 0.86602540378443864676d0
2104 integer :: i,j,ind,inres,mnum
2105 ! Revised 3/31/05 AL: correlation between random contributions to
2106 ! position and velocity increments included.
2107 ! The correlation coefficients are calculated at low-friction limit.
2108 ! Also, friction forces are now not calculated with new velocities.
2110 ! call friction_force
2111 call stochastic_force(stochforcvec)
2113 ! Compute the acceleration due to friction forces (d_af_work) and stochastic
2114 ! forces (d_as_work)
2116 call ginv_mult(stochforcvec, d_as_work1)
2122 d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j)) &
2123 +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
2128 d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j)) &
2129 +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
2135 ! iti=iabs(itype(i,mnum))
2136 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2137 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2138 .and.(mnum.lt.4)) then
2141 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres) &
2142 +d_af_work(ind+j))+sin60*d_as_work(ind+j) &
2143 +cos60*d_as_work1(ind+j))*d_time
2149 end subroutine sddir_verlet2
2150 !-----------------------------------------------------------------------------
2151 subroutine max_accel
2153 ! Find the maximum difference in the accelerations of the the sites
2154 ! at the beginning and the end of the time step.
2157 ! implicit real*8 (a-h,o-z)
2158 ! include 'DIMENSIONS'
2159 ! include 'COMMON.CONTROL'
2160 ! include 'COMMON.VAR'
2161 ! include 'COMMON.MD'
2162 ! include 'COMMON.CHAIN'
2163 ! include 'COMMON.DERIV'
2164 ! include 'COMMON.GEO'
2165 ! include 'COMMON.LOCAL'
2166 ! include 'COMMON.INTERACT'
2167 ! include 'COMMON.IOUNITS'
2168 real(kind=8),dimension(3) :: aux,accel,accel_old
2169 real(kind=8) :: dacc
2173 ! aux(j)=d_a(j,0)-d_a_old(j,0)
2174 accel_old(j)=d_a_old(j,0)
2181 ! 7/3/08 changed to asymmetric difference
2183 ! accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
2184 accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
2185 accel(j)=accel(j)+0.5d0*d_a(j,i)
2186 ! if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
2187 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
2188 dacc=dabs(accel(j)-accel_old(j))
2189 ! write (iout,*) i,dacc
2190 if (dacc.gt.amax) amax=dacc
2198 accel_old(j)=d_a_old(j,0)
2203 accel_old(j)=accel_old(j)+d_a_old(j,1)
2204 accel(j)=accel(j)+d_a(j,1)
2209 ! iti=iabs(itype(i,mnum))
2210 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2211 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2212 .and.(mnum.lt.4)) then
2214 ! accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
2215 accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
2216 accel(j)=accel(j)+d_a(j,i+nres)
2220 ! if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
2221 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
2222 dacc=dabs(accel(j)-accel_old(j))
2223 ! write (iout,*) "side-chain",i,dacc
2224 if (dacc.gt.amax) amax=dacc
2228 accel_old(j)=accel_old(j)+d_a_old(j,i)
2229 accel(j)=accel(j)+d_a(j,i)
2230 ! aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
2234 end subroutine max_accel
2235 !-----------------------------------------------------------------------------
2236 subroutine predict_edrift(epdrift)
2238 ! Predict the drift of the potential energy
2241 use control_data, only: lmuca
2242 ! implicit real*8 (a-h,o-z)
2243 ! include 'DIMENSIONS'
2244 ! include 'COMMON.CONTROL'
2245 ! include 'COMMON.VAR'
2246 ! include 'COMMON.MD'
2247 ! include 'COMMON.CHAIN'
2248 ! include 'COMMON.DERIV'
2249 ! include 'COMMON.GEO'
2250 ! include 'COMMON.LOCAL'
2251 ! include 'COMMON.INTERACT'
2252 ! include 'COMMON.IOUNITS'
2253 ! include 'COMMON.MUCA'
2254 real(kind=8) :: epdrift,epdriftij
2256 ! Drift of the potential energy
2262 epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
2263 if (lmuca) epdriftij=epdriftij*factor
2264 ! write (iout,*) "back",i,j,epdriftij
2265 if (epdriftij.gt.epdrift) epdrift=epdriftij
2269 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1.and.&
2270 molnum(i).lt.4) then
2273 dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
2274 if (lmuca) epdriftij=epdriftij*factor
2275 ! write (iout,*) "side",i,j,epdriftij
2276 if (epdriftij.gt.epdrift) epdrift=epdriftij
2280 epdrift=0.5d0*epdrift*d_time*d_time
2281 ! write (iout,*) "epdrift",epdrift
2283 end subroutine predict_edrift
2284 !-----------------------------------------------------------------------------
2285 subroutine verlet_bath
2287 ! Coupling to the thermostat by using the Berendsen algorithm
2290 ! implicit real*8 (a-h,o-z)
2291 ! include 'DIMENSIONS'
2292 ! include 'COMMON.CONTROL'
2293 ! include 'COMMON.VAR'
2294 ! include 'COMMON.MD'
2295 ! include 'COMMON.CHAIN'
2296 ! include 'COMMON.DERIV'
2297 ! include 'COMMON.GEO'
2298 ! include 'COMMON.LOCAL'
2299 ! include 'COMMON.INTERACT'
2300 ! include 'COMMON.IOUNITS'
2301 ! include 'COMMON.NAMES'
2302 real(kind=8) :: T_half,fact
2303 integer :: i,j,inres,mnum
2305 T_half=2.0d0/(dimen3*Rb)*EK
2306 fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
2307 ! write(iout,*) "T_half", T_half
2308 ! write(iout,*) "EK", EK
2309 ! write(iout,*) "fact", fact
2311 d_t(j,0)=fact*d_t(j,0)
2315 d_t(j,i)=fact*d_t(j,i)
2320 ! iti=iabs(itype(i,mnum))
2321 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2322 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2323 .and.(mnum.lt.4)) then
2326 d_t(j,inres)=fact*d_t(j,inres)
2331 end subroutine verlet_bath
2332 !-----------------------------------------------------------------------------
2334 ! Set up the initial conditions of a MD simulation
2337 use control, only:tcpu
2338 !el use io_basic, only:ilen
2341 use minimm, only:minim_dc,minimize,sc_move
2342 use io_config, only:readrst
2343 use io, only:statout
2344 ! implicit real*8 (a-h,o-z)
2345 ! include 'DIMENSIONS'
2348 character(len=16) :: form
2349 integer :: IERROR,ERRCODE
2351 integer :: iranmin,itrial,itmp
2352 ! include 'COMMON.SETUP'
2353 ! include 'COMMON.CONTROL'
2354 ! include 'COMMON.VAR'
2355 ! include 'COMMON.MD'
2357 ! include 'COMMON.LANGEVIN'
2359 ! include 'COMMON.LANGEVIN.lang0'
2361 ! include 'COMMON.CHAIN'
2362 ! include 'COMMON.DERIV'
2363 ! include 'COMMON.GEO'
2364 ! include 'COMMON.LOCAL'
2365 ! include 'COMMON.INTERACT'
2366 ! include 'COMMON.IOUNITS'
2367 ! include 'COMMON.NAMES'
2368 ! include 'COMMON.REMD'
2369 real(kind=8),dimension(0:n_ene) :: energia_long,energia_short
2370 real(kind=8),dimension(3) :: vcm,incr,L
2371 real(kind=8) :: xv,sigv,lowb,highb
2372 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
2373 character(len=256) :: qstr
2376 character(len=50) :: tytul
2377 logical :: file_exist
2378 !el common /gucio/ cm
2379 integer :: i,j,ipos,iq,iw,nft_sc,iretcode,nfun,ierr,mnum,itime
2380 real(kind=8) :: etot,tt0
2384 ! write(iout,*) "d_time", d_time
2385 ! Compute the standard deviations of stochastic forces for Langevin dynamics
2386 ! if the friction coefficients do not depend on surface area
2387 if (lang.gt.0 .and. .not.surfarea) then
2390 stdforcp(i)=stdfp(mnum)*dsqrt(gamp(mnum))
2394 stdforcsc(i)=stdfsc(iabs(itype(i,mnum)),mnum) &
2395 *dsqrt(gamsc(iabs(itype(i,mnum)),mnum))
2399 ! Open the pdb file for snapshotshots
2402 if (ilen(tmpdir).gt.0) &
2403 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2404 liczba(:ilen(liczba))//".pdb")
2406 file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2410 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) &
2411 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2412 liczba(:ilen(liczba))//".x")
2413 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2416 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) &
2417 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2418 liczba(:ilen(liczba))//".cx")
2419 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2425 if (ilen(tmpdir).gt.0) &
2426 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
2427 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
2429 if (ilen(tmpdir).gt.0) &
2430 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
2431 cartname=prefix(:ilen(prefix))//"_MD.cx"
2435 write (qstr,'(256(1h ))')
2438 iq = qinfrag(i,iset)*10
2439 iw = wfrag(i,iset)/100
2441 if(me.eq.king.or..not.out1file) &
2442 write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
2443 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
2448 iq = qinpair(i,iset)*10
2449 iw = wpair(i,iset)/100
2451 if(me.eq.king.or..not.out1file) &
2452 write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
2453 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
2457 ! pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
2459 ! cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
2461 ! cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
2463 ! statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
2467 if (restart1file) then
2469 inquire(file=mremd_rst_name,exist=file_exist)
2470 write (*,*) me," Before broadcast: file_exist",file_exist
2472 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,&
2475 write (*,*) me," After broadcast: file_exist",file_exist
2476 ! inquire(file=mremd_rst_name,exist=file_exist)
2477 if(me.eq.king.or..not.out1file) &
2478 write(iout,*) "Initial state read by master and distributed"
2480 if (ilen(tmpdir).gt.0) &
2481 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_' &
2482 //liczba(:ilen(liczba))//'.rst')
2483 inquire(file=rest2name,exist=file_exist)
2486 if(.not.restart1file) then
2487 if(me.eq.king.or..not.out1file) &
2488 write(iout,*) "Initial state will be read from file ",&
2489 rest2name(:ilen(rest2name))
2492 call rescale_weights(t_bath)
2494 if(me.eq.king.or..not.out1file)then
2495 if (restart1file) then
2496 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),&
2499 write(iout,*) "File ",rest2name(:ilen(rest2name)),&
2502 write(iout,*) "Initial velocities randomly generated"
2509 ! Generate initial velocities
2510 if(me.eq.king.or..not.out1file) &
2511 write(iout,*) "Initial velocities randomly generated"
2516 ! rest2name = prefix(:ilen(prefix))//'.rst'
2517 if(me.eq.king.or..not.out1file)then
2518 write (iout,*) "Initial velocities"
2520 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2521 (d_t(j,i+nres),j=1,3)
2523 ! Zeroing the total angular momentum of the system
2524 write(iout,*) "Calling the zero-angular momentum subroutine"
2527 ! Getting the potential energy and forces and velocities and accelerations
2529 ! write (iout,*) "velocity of the center of the mass:"
2530 ! write (iout,*) (vcm(j),j=1,3)
2532 d_t(j,0)=d_t(j,0)-vcm(j)
2534 ! Removing the velocity of the center of mass
2536 if(me.eq.king.or..not.out1file)then
2537 write (iout,*) "vcm right after adjustment:"
2538 write (iout,*) (vcm(j),j=1,3)
2540 if ((.not.rest).or.(forceminim)) then
2541 if (forceminim) call chainbuild_cart
2542 if(iranconf.ne.0 .or.indpdb.gt.0.and..not.unres_pdb .or.preminim) then
2544 print *, 'Calling OVERLAP_SC'
2545 call overlap_sc(fail)
2546 print *,'after OVERLAP'
2549 print *,'call SC_MOVE'
2550 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
2551 print *,'SC_move',nft_sc,etot
2552 if(me.eq.king.or..not.out1file) &
2553 write(iout,*) 'SC_move',nft_sc,etot
2557 print *, 'Calling MINIM_DC'
2558 call minim_dc(etot,iretcode,nfun)
2560 call geom_to_var(nvar,varia)
2561 print *,'Calling MINIMIZE.'
2562 call minimize(etot,varia,iretcode,nfun)
2563 call var_to_geom(nvar,varia)
2565 write(iout,*) "just before minimin"
2567 if(me.eq.king.or..not.out1file) &
2568 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2570 write(iout,*) "just after minimin"
2572 if(iranconf.ne.0) then
2573 !c 8/22/17 AL Loop to produce a low-energy random conformation
2576 if(me.eq.king.or..not.out1file) &
2577 write (iout,*) 'Calling OVERLAP_SC'
2578 call overlap_sc(fail)
2579 endif !endif overlap
2582 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
2583 print *,'SC_move',nft_sc,etot
2584 if(me.eq.king.or..not.out1file) &
2585 write(iout,*) 'SC_move',nft_sc,etot
2589 print *, 'Calling MINIM_DC'
2590 call minim_dc(etot,iretcode,nfun)
2591 call int_from_cart1(.false.)
2593 call geom_to_var(nvar,varia)
2594 print *,'Calling MINIMIZE.'
2595 call minimize(etot,varia,iretcode,nfun)
2596 call var_to_geom(nvar,varia)
2598 if(me.eq.king.or..not.out1file) &
2599 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2600 write(iout,*) "just after minimin"
2602 if (isnan(etot) .or. etot.gt.4.0d6) then
2603 write (iout,*) "Energy too large",etot, &
2604 " trying another random conformation"
2607 call gen_rand_conf(itmp,*30)
2609 30 write (iout,*) 'Failed to generate random conformation', &
2611 write (*,*) 'Processor:',me, &
2612 ' Failed to generate random conformation',&
2621 write (iout,'(a,i3,a)') 'Processor:',me, &
2622 ' error in generating random conformation.'
2623 write (*,'(a,i3,a)') 'Processor:',me, &
2624 ' error in generating random conformation.'
2627 ! call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
2628 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
2638 write (iout,'(a,i3,a)') 'Processor:',me, &
2639 ' failed to generate a low-energy random conformation.'
2640 write (*,'(a,i3,a,f10.3)') 'Processor:',me, &
2641 ' failed to generate a low-energy random conformation.',etot
2645 ! call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
2646 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
2653 call chainbuild_cart
2655 write(iout,*) "just after kinetic"
2660 kinetic_T=2.0d0/(dimen3*Rb)*EK
2661 if(me.eq.king.or..not.out1file)then
2662 write(iout,*) "just after verlet_bath"
2672 write(iout,*) "before ETOTAL"
2673 call etotal(potEcomp)
2674 if (large) call enerprint(potEcomp)
2677 t_etotal=t_etotal+MPI_Wtime()-tt0
2679 t_etotal=t_etotal+tcpu()-tt0
2686 if (amax*d_time .gt. dvmax) then
2687 d_time=d_time*dvmax/amax
2688 if(me.eq.king.or..not.out1file) write (iout,*) &
2689 "Time step reduced to",d_time,&
2690 " because of too large initial acceleration."
2692 if(me.eq.king.or..not.out1file)then
2693 write(iout,*) "Potential energy and its components"
2694 call enerprint(potEcomp)
2695 ! write(iout,*) (potEcomp(i),i=0,n_ene)
2697 potE=potEcomp(0)-potEcomp(20)
2701 if (ntwe.ne.0) call statout(itime)
2702 if(me.eq.king.or..not.out1file) &
2703 write (iout,'(/a/3(a25,1pe14.5/))') "Initial:", &
2704 " Kinetic energy",EK," Potential energy",potE, &
2705 " Total energy",totE," Maximum acceleration ", &
2708 write (iout,*) "Initial coordinates"
2710 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),&
2713 write (iout,*) "Initial dC"
2715 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),&
2716 (dc(j,i+nres),j=1,3)
2718 write (iout,*) "Initial velocities"
2719 write (iout,"(13x,' backbone ',23x,' side chain')")
2721 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2722 (d_t(j,i+nres),j=1,3)
2724 write (iout,*) "Initial accelerations"
2726 ! write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2727 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),&
2728 (d_a(j,i+nres),j=1,3)
2734 d_t_old(j,i)=d_t(j,i)
2735 d_a_old(j,i)=d_a(j,i)
2737 ! write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
2746 call etotal_short(energia_short)
2747 if (large) call enerprint(potEcomp)
2750 t_eshort=t_eshort+MPI_Wtime()-tt0
2752 t_eshort=t_eshort+tcpu()-tt0
2757 if(.not.out1file .and. large) then
2758 write (iout,*) "energia_long",energia_long(0),&
2759 " energia_short",energia_short(0),&
2760 " total",energia_long(0)+energia_short(0)
2761 write (iout,*) "Initial fast-force accelerations"
2763 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
2764 (d_a(j,i+nres),j=1,3)
2767 ! 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
2770 d_a_short(j,i)=d_a(j,i)
2779 call etotal_long(energia_long)
2780 if (large) call enerprint(potEcomp)
2783 t_elong=t_elong+MPI_Wtime()-tt0
2785 t_elong=t_elong+tcpu()-tt0
2790 if(.not.out1file .and. large) then
2791 write (iout,*) "energia_long",energia_long(0)
2792 write (iout,*) "Initial slow-force accelerations"
2794 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
2795 (d_a(j,i+nres),j=1,3)
2799 t_enegrad=t_enegrad+MPI_Wtime()-tt0
2801 t_enegrad=t_enegrad+tcpu()-tt0
2805 end subroutine init_MD
2806 !-----------------------------------------------------------------------------
2807 subroutine random_vel
2809 ! implicit real*8 (a-h,o-z)
2811 use random, only:anorm_distr
2813 ! include 'DIMENSIONS'
2814 ! include 'COMMON.CONTROL'
2815 ! include 'COMMON.VAR'
2816 ! include 'COMMON.MD'
2818 ! include 'COMMON.LANGEVIN'
2820 ! include 'COMMON.LANGEVIN.lang0'
2822 ! include 'COMMON.CHAIN'
2823 ! include 'COMMON.DERIV'
2824 ! include 'COMMON.GEO'
2825 ! include 'COMMON.LOCAL'
2826 ! include 'COMMON.INTERACT'
2827 ! include 'COMMON.IOUNITS'
2828 ! include 'COMMON.NAMES'
2829 ! include 'COMMON.TIME1'
2830 real(kind=8) :: xv,sigv,lowb,highb ,Ek1
2833 real(kind=8) ,allocatable, dimension(:) :: DDU1,DDU2,DL2,DL1,xsolv,DML,rs
2834 real(kind=8) :: sumx
2836 real(kind=8) ,allocatable, dimension(:) :: rsold
2837 real (kind=8),allocatable,dimension(:,:) :: matold
2841 integer :: i,j,ii,k,ind,mark,imark,mnum
2842 ! Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
2843 ! First generate velocities in the eigenspace of the G matrix
2844 ! write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
2851 sigv=dsqrt((Rb*t_bath)/geigen(i))
2854 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
2856 write (iout,*) "i",i," ii",ii," geigen",geigen(i),&
2857 " d_t_work_new",d_t_work_new(ii)
2868 Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
2871 write (iout,*) "Ek from eigenvectors",Ek1
2872 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
2876 ! Transform velocities to UNRES coordinate space
2877 allocate (DL1(2*nres))
2878 allocate (DDU1(2*nres))
2879 allocate (DL2(2*nres))
2880 allocate (DDU2(2*nres))
2881 allocate (xsolv(2*nres))
2882 allocate (DML(2*nres))
2883 allocate (rs(2*nres))
2885 allocate (rsold(2*nres))
2886 allocate (matold(2*nres,2*nres))
2888 matold(1,1)=DMorig(1)
2889 matold(1,2)=DU1orig(1)
2890 matold(1,3)=DU2orig(1)
2891 write (*,*) DMorig(1),DU1orig(1),DU2orig(1)
2896 matold(i,j)=DMorig(i)
2897 matold(i,j-1)=DU1orig(i-1)
2899 matold(i,j-2)=DU2orig(i-2)
2907 matold(i,j+1)=DU1orig(i)
2913 matold(i,j+2)=DU2orig(i)
2917 matold(dimen,dimen)=DMorig(dimen)
2918 matold(dimen,dimen-1)=DU1orig(dimen-1)
2919 matold(dimen,dimen-2)=DU2orig(dimen-2)
2920 write(iout,*) "old gmatrix"
2921 call matout(dimen,dimen,2*nres,2*nres,matold)
2925 ! Find the ith eigenvector of the pentadiagonal inertiq matrix
2929 DML(j)=DMorig(j)-geigen(i)
2932 DML(j-1)=DMorig(j)-geigen(i)
2937 DDU1(imark-1)=DU2orig(imark-1)
2938 do j=imark+1,dimen-1
2939 DDU1(j-1)=DU1orig(j)
2947 DDU2(j)=DU2orig(j+1)
2956 write (iout,*) "DL2,DL1,DML,DDU1,DDU2"
2957 write(iout,'(10f10.5)') (DL2(k),k=3,dimen-1)
2958 write(iout,'(10f10.5)') (DL1(k),k=2,dimen-1)
2959 write(iout,'(10f10.5)') (DML(k),k=1,dimen-1)
2960 write(iout,'(10f10.5)') (DDU1(k),k=1,dimen-2)
2961 write(iout,'(10f10.5)') (DDU2(k),k=1,dimen-3)
2964 if (imark.gt.2) rs(imark-2)=-DU2orig(imark-2)
2965 if (imark.gt.1) rs(imark-1)=-DU1orig(imark-1)
2966 if (imark.lt.dimen) rs(imark)=-DU1orig(imark)
2967 if (imark.lt.dimen-1) rs(imark+1)=-DU2orig(imark)
2971 ! write (iout,*) "Vector rs"
2973 ! write (iout,*) j,rs(j)
2976 call FDIAG(dimen-1,DL2,DL1,DML,DDU1,DDU2,rs,xsolv,mark)
2983 sumx=-geigen(i)*xsolv(j)
2985 sumx=sumx+matold(j,k)*xsolv(k)
2988 sumx=sumx+matold(j,k)*xsolv(k-1)
2990 write(iout,'(i5,3f10.5)') j,sumx,rsold(j),sumx-rsold(j)
2993 sumx=-geigen(i)*xsolv(j-1)
2995 sumx=sumx+matold(j,k)*xsolv(k)
2998 sumx=sumx+matold(j,k)*xsolv(k-1)
3000 write(iout,'(i5,3f10.5)') j-1,sumx,rsold(j-1),sumx-rsold(j-1)
3004 "Solution of equations system",i," for eigenvalue",geigen(i)
3006 write(iout,'(i5,f10.5)') j,xsolv(j)
3009 do j=dimen-1,imark,-1
3014 write (iout,*) "Un-normalized eigenvector",i," for eigenvalue",geigen(i)
3016 write(iout,'(i5,f10.5)') j,xsolv(j)
3019 ! Normalize ith eigenvector
3022 sumx=sumx+xsolv(j)**2
3026 xsolv(j)=xsolv(j)/sumx
3029 write (iout,*) "Eigenvector",i," for eigenvalue",geigen(i)
3031 write(iout,'(i5,3f10.5)') j,xsolv(j)
3034 ! All done at this point for eigenvector i, exit loop
3042 write (iout,*) "Unable to find eigenvector",i
3045 ! write (iout,*) "i=",i
3047 ! write (iout,*) "k=",k
3050 ! write(iout,*) "ind",ind," ind1",3*(i-1)+k
3051 d_t_work(ind)=d_t_work(ind) &
3052 +xsolv(j)*d_t_work_new(3*(i-1)+k)
3055 enddo ! i (loop over the eigenvectors)
3058 write (iout,*) "d_t_work"
3060 write (iout,"(i5,f10.5)") i,d_t_work(i)
3065 ! if (itype(i,1).eq.10) then
3067 if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
3068 iti=iabs(itype(i,mnum))
3069 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3070 if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
3071 .or.(mnum.ge.4)) then
3078 ! write (iout,*) "k",k," ii+k",ii+k," ii+j+k",ii+j+k,"EK1",Ek1
3079 Ek1=Ek1+0.5d0*mp(mnum)*((d_t_work(ii+j+k)+d_t_work_new(ii+k))/2)**2+&
3080 0.5d0*0.25d0*IP(mnum)*(d_t_work(ii+j+k)-d_t_work_new(ii+k))**2
3083 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3084 .and.(mnum.lt.4)) ii=ii+3
3085 write (iout,*) "i",i," itype",itype(i,mnum)," mass",msc(itype(i,mnum),mnum)
3086 write (iout,*) "ii",ii
3089 write (iout,*) "k",k," ii",ii,"EK1",EK1
3090 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3092 Ek1=Ek1+0.5d0*Isc(iabs(itype(i,mnum)),mnum)*(d_t_work(ii)-d_t_work(ii-3))**2
3093 Ek1=Ek1+0.5d0*msc(iabs(itype(i,mnum)),mnum)*d_t_work(ii)**2
3095 write (iout,*) "i",i," ii",ii
3097 write (iout,*) "Ek from d_t_work",Ek1
3098 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
3100 if(allocated(DDU1)) deallocate(DDU1)
3101 if(allocated(DDU2)) deallocate(DDU2)
3102 if(allocated(DL2)) deallocate(DL2)
3103 if(allocated(DL1)) deallocate(DL1)
3104 if(allocated(xsolv)) deallocate(xsolv)
3105 if(allocated(DML)) deallocate(DML)
3106 if(allocated(rs)) deallocate(rs)
3108 if(allocated(matold)) deallocate(matold)
3109 if(allocated(rsold)) deallocate(rsold)
3114 d_t(k,j)=d_t_work(ind)
3118 if (itype(j,1).ne.10 .and. itype(j,mnum).ne.ntyp1_molec(mnum)&
3119 .and.(mnum.lt.4)) then
3121 d_t(k,j+nres)=d_t_work(ind)
3127 write (iout,*) "Random velocities in the Calpha,SC space"
3129 write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3)
3132 write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3)
3139 ! if (itype(i,1).eq.10) then
3141 if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
3142 .or.(mnum.ge.4)) then
3144 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
3148 d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
3149 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
3154 write (iout,*)"Random velocities in the virtual-bond-vector space"
3156 write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3)
3159 write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3)
3162 write (iout,*) "Ek from d_t_work",Ek1
3163 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
3171 d_t_work(ind)=d_t_work(ind) &
3172 +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
3174 ! write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
3178 ! Transfer to the d_t vector
3180 d_t(j,0)=d_t_work(j)
3186 d_t(j,i)=d_t_work(ind)
3191 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3192 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3193 .and.(mnum.lt.4)) then
3196 d_t(j,i+nres)=d_t_work(ind)
3202 ! write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",&
3203 ! 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
3205 ! write(iout,*) "end init MD"
3207 end subroutine random_vel
3208 !-----------------------------------------------------------------------------
3210 subroutine sd_verlet_p_setup
3211 ! Sets up the parameters of stochastic Verlet algorithm
3212 ! implicit real*8 (a-h,o-z)
3213 ! include 'DIMENSIONS'
3214 use control, only: tcpu
3219 ! include 'COMMON.CONTROL'
3220 ! include 'COMMON.VAR'
3221 ! include 'COMMON.MD'
3223 ! include 'COMMON.LANGEVIN'
3225 ! include 'COMMON.LANGEVIN.lang0'
3227 ! include 'COMMON.CHAIN'
3228 ! include 'COMMON.DERIV'
3229 ! include 'COMMON.GEO'
3230 ! include 'COMMON.LOCAL'
3231 ! include 'COMMON.INTERACT'
3232 ! include 'COMMON.IOUNITS'
3233 ! include 'COMMON.NAMES'
3234 ! include 'COMMON.TIME1'
3235 real(kind=8),dimension(6*nres) :: emgdt !(MAXRES6) maxres6=6*maxres
3236 real(kind=8) :: pterm,vterm,rho,rhoc,vsig
3237 real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
3238 prand_vec,vrand_vec1,vrand_vec2 !(MAXRES6) maxres6=6*maxres
3239 logical :: lprn = .false.
3240 real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
3241 real(kind=8) :: ktm,gdt,egdt,gdt2,gdt3,gdt4,gdt5,gdt6,gdt7,gdt8,&
3243 integer :: i,maxres2
3250 ! AL 8/17/04 Code adapted from tinker
3252 ! Get the frictional and random terms for stochastic dynamics in the
3253 ! eigenspace of mass-scaled UNRES friction matrix
3257 gdt = fricgam(i) * d_time
3259 ! Stochastic dynamics reduces to simple MD for zero friction
3261 if (gdt .le. zero) then
3262 pfric_vec(i) = 1.0d0
3263 vfric_vec(i) = d_time
3264 afric_vec(i) = 0.5d0 * d_time * d_time
3265 prand_vec(i) = 0.0d0
3266 vrand_vec1(i) = 0.0d0
3267 vrand_vec2(i) = 0.0d0
3269 ! Analytical expressions when friction coefficient is large
3272 if (gdt .ge. gdt_radius) then
3275 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
3276 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
3277 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
3278 vterm = 1.0d0 - egdt**2
3279 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
3281 ! Use series expansions when friction coefficient is small
3292 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0 &
3293 - gdt5/120.0d0 + gdt6/720.0d0 &
3294 - gdt7/5040.0d0 + gdt8/40320.0d0 &
3295 - gdt9/362880.0d0) / fricgam(i)**2
3296 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
3297 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
3298 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0 &
3299 + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0 &
3300 + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0 &
3301 + 127.0d0*gdt9/90720.0d0
3302 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0 &
3303 - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0 &
3304 - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0 &
3305 - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
3306 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0 &
3307 - 17.0d0*gdt2/1280.0d0 &
3308 + 17.0d0*gdt3/6144.0d0 &
3309 + 40967.0d0*gdt4/34406400.0d0 &
3310 - 57203.0d0*gdt5/275251200.0d0 &
3311 - 1429487.0d0*gdt6/13212057600.0d0)
3314 ! Compute the scaling factors of random terms for the nonzero friction case
3316 ktm = 0.5d0*d_time/fricgam(i)
3317 psig = dsqrt(ktm*pterm) / fricgam(i)
3318 vsig = dsqrt(ktm*vterm)
3319 rhoc = dsqrt(1.0d0 - rho*rho)
3321 vrand_vec1(i) = vsig * rho
3322 vrand_vec2(i) = vsig * rhoc
3327 "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",&
3330 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),&
3331 afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
3335 ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
3338 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
3339 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
3340 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
3341 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
3342 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
3343 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
3346 t_sdsetup=t_sdsetup+MPI_Wtime()
3348 t_sdsetup=t_sdsetup+tcpu()-tt0
3351 end subroutine sd_verlet_p_setup
3352 !-----------------------------------------------------------------------------
3353 subroutine eigtransf1(n,ndim,ab,d,c)
3357 real(kind=8) :: ab(ndim,ndim,n),c(ndim,n),d(ndim)
3363 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
3368 end subroutine eigtransf1
3369 !-----------------------------------------------------------------------------
3370 subroutine eigtransf(n,ndim,a,b,d,c)
3374 real(kind=8) :: a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
3380 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
3385 end subroutine eigtransf
3386 !-----------------------------------------------------------------------------
3387 subroutine sd_verlet1
3389 ! Applying stochastic velocity Verlet algorithm - step 1 to velocities
3391 ! implicit real*8 (a-h,o-z)
3392 ! include 'DIMENSIONS'
3393 ! include 'COMMON.CONTROL'
3394 ! include 'COMMON.VAR'
3395 ! include 'COMMON.MD'
3397 ! include 'COMMON.LANGEVIN'
3399 ! include 'COMMON.LANGEVIN.lang0'
3401 ! include 'COMMON.CHAIN'
3402 ! include 'COMMON.DERIV'
3403 ! include 'COMMON.GEO'
3404 ! include 'COMMON.LOCAL'
3405 ! include 'COMMON.INTERACT'
3406 ! include 'COMMON.IOUNITS'
3407 ! include 'COMMON.NAMES'
3408 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
3409 !el common /stochcalc/ stochforcvec
3410 logical :: lprn = .false.
3411 real(kind=8) :: ddt1,ddt2
3412 integer :: i,j,ind,inres
3414 ! write (iout,*) "dc_old"
3416 ! write (iout,'(i5,3f10.5,5x,3f10.5)')
3417 ! & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3420 dc_work(j)=dc_old(j,0)
3421 d_t_work(j)=d_t_old(j,0)
3422 d_a_work(j)=d_a_old(j,0)
3427 dc_work(ind+j)=dc_old(j,i)
3428 d_t_work(ind+j)=d_t_old(j,i)
3429 d_a_work(ind+j)=d_a_old(j,i)
3435 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3436 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3437 .and.(mnum.lt.4)) then
3439 dc_work(ind+j)=dc_old(j,i+nres)
3440 d_t_work(ind+j)=d_t_old(j,i+nres)
3441 d_a_work(ind+j)=d_a_old(j,i+nres)
3449 "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",&
3453 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
3454 vfric_mat(i,j),afric_mat(i,j),&
3455 prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3463 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) &
3464 +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3465 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3466 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3468 d_t_work_new(i)=ddt1+0.5d0*ddt2
3469 d_t_work(i)=ddt1+ddt2
3474 d_t(j,0)=d_t_work(j)
3479 dc(j,i)=dc_work(ind+j)
3480 d_t(j,i)=d_t_work(ind+j)
3486 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3487 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3488 .and.(mnum.lt.4)) then
3491 dc(j,inres)=dc_work(ind+j)
3492 d_t(j,inres)=d_t_work(ind+j)
3498 end subroutine sd_verlet1
3499 !-----------------------------------------------------------------------------
3500 subroutine sd_verlet2
3502 ! Calculating the adjusted velocities for accelerations
3504 ! implicit real*8 (a-h,o-z)
3505 ! include 'DIMENSIONS'
3506 ! include 'COMMON.CONTROL'
3507 ! include 'COMMON.VAR'
3508 ! include 'COMMON.MD'
3510 ! include 'COMMON.LANGEVIN'
3512 ! include 'COMMON.LANGEVIN.lang0'
3514 ! include 'COMMON.CHAIN'
3515 ! include 'COMMON.DERIV'
3516 ! include 'COMMON.GEO'
3517 ! include 'COMMON.LOCAL'
3518 ! include 'COMMON.INTERACT'
3519 ! include 'COMMON.IOUNITS'
3520 ! include 'COMMON.NAMES'
3521 !el real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV !(MAXRES6) maxres6=6*maxres
3522 real(kind=8),dimension(6*nres) :: stochforcvecV !(MAXRES6) maxres6=6*maxres
3523 !el common /stochcalc/ stochforcvec
3525 real(kind=8) :: ddt1,ddt2
3526 integer :: i,j,ind,inres
3527 ! Compute the stochastic forces which contribute to velocity change
3529 call stochastic_force(stochforcvecV)
3536 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3537 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+ &
3538 vrand_mat2(i,j)*stochforcvecV(j)
3540 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
3544 d_t(j,0)=d_t_work(j)
3549 d_t(j,i)=d_t_work(ind+j)
3555 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3556 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3557 .and.(mnum.lt.4)) then
3560 d_t(j,inres)=d_t_work(ind+j)
3566 end subroutine sd_verlet2
3567 !-----------------------------------------------------------------------------
3568 subroutine sd_verlet_ciccotti_setup
3570 ! Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
3572 ! implicit real*8 (a-h,o-z)
3573 ! include 'DIMENSIONS'
3574 use control, only: tcpu
3579 ! include 'COMMON.CONTROL'
3580 ! include 'COMMON.VAR'
3581 ! include 'COMMON.MD'
3583 ! include 'COMMON.LANGEVIN'
3585 ! include 'COMMON.LANGEVIN.lang0'
3587 ! include 'COMMON.CHAIN'
3588 ! include 'COMMON.DERIV'
3589 ! include 'COMMON.GEO'
3590 ! include 'COMMON.LOCAL'
3591 ! include 'COMMON.INTERACT'
3592 ! include 'COMMON.IOUNITS'
3593 ! include 'COMMON.NAMES'
3594 ! include 'COMMON.TIME1'
3595 real(kind=8),dimension(6*nres) :: emgdt !(MAXRES6) maxres6=6*maxres
3596 real(kind=8) :: pterm,vterm,rho,rhoc,vsig
3597 real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
3598 prand_vec,vrand_vec1,vrand_vec2 !(MAXRES6) maxres6=6*maxres
3599 logical :: lprn = .false.
3600 real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
3601 real(kind=8) :: ktm,gdt,egdt,tt0
3602 integer :: i,maxres2
3609 ! AL 8/17/04 Code adapted from tinker
3611 ! Get the frictional and random terms for stochastic dynamics in the
3612 ! eigenspace of mass-scaled UNRES friction matrix
3616 write (iout,*) "i",i," fricgam",fricgam(i)
3617 gdt = fricgam(i) * d_time
3619 ! Stochastic dynamics reduces to simple MD for zero friction
3621 if (gdt .le. zero) then
3622 pfric_vec(i) = 1.0d0
3623 vfric_vec(i) = d_time
3624 afric_vec(i) = 0.5d0*d_time*d_time
3625 prand_vec(i) = afric_vec(i)
3626 vrand_vec2(i) = vfric_vec(i)
3628 ! Analytical expressions when friction coefficient is large
3633 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
3634 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
3635 prand_vec(i) = afric_vec(i)
3636 vrand_vec2(i) = vfric_vec(i)
3638 ! Compute the scaling factors of random terms for the nonzero friction case
3640 ! ktm = 0.5d0*d_time/fricgam(i)
3641 ! psig = dsqrt(ktm*pterm) / fricgam(i)
3642 ! vsig = dsqrt(ktm*vterm)
3643 ! prand_vec(i) = psig*afric_vec(i)
3644 ! vrand_vec2(i) = vsig*vfric_vec(i)
3649 "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",&
3652 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),&
3653 afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
3657 ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
3659 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
3660 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
3661 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
3662 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
3663 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
3665 t_sdsetup=t_sdsetup+MPI_Wtime()
3667 t_sdsetup=t_sdsetup+tcpu()-tt0
3670 end subroutine sd_verlet_ciccotti_setup
3671 !-----------------------------------------------------------------------------
3672 subroutine sd_verlet1_ciccotti
3674 ! Applying stochastic velocity Verlet algorithm - step 1 to velocities
3675 ! implicit real*8 (a-h,o-z)
3677 ! include 'DIMENSIONS'
3681 ! include 'COMMON.CONTROL'
3682 ! include 'COMMON.VAR'
3683 ! include 'COMMON.MD'
3685 ! include 'COMMON.LANGEVIN'
3687 ! include 'COMMON.LANGEVIN.lang0'
3689 ! include 'COMMON.CHAIN'
3690 ! include 'COMMON.DERIV'
3691 ! include 'COMMON.GEO'
3692 ! include 'COMMON.LOCAL'
3693 ! include 'COMMON.INTERACT'
3694 ! include 'COMMON.IOUNITS'
3695 ! include 'COMMON.NAMES'
3696 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
3697 !el common /stochcalc/ stochforcvec
3698 logical :: lprn = .false.
3699 real(kind=8) :: ddt1,ddt2
3700 integer :: i,j,ind,inres
3701 ! write (iout,*) "dc_old"
3703 ! write (iout,'(i5,3f10.5,5x,3f10.5)')
3704 ! & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3707 dc_work(j)=dc_old(j,0)
3708 d_t_work(j)=d_t_old(j,0)
3709 d_a_work(j)=d_a_old(j,0)
3714 dc_work(ind+j)=dc_old(j,i)
3715 d_t_work(ind+j)=d_t_old(j,i)
3716 d_a_work(ind+j)=d_a_old(j,i)
3721 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3723 dc_work(ind+j)=dc_old(j,i+nres)
3724 d_t_work(ind+j)=d_t_old(j,i+nres)
3725 d_a_work(ind+j)=d_a_old(j,i+nres)
3734 "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",&
3738 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
3739 vfric_mat(i,j),afric_mat(i,j),&
3740 prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3748 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) &
3749 +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3750 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3751 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3753 d_t_work_new(i)=ddt1+0.5d0*ddt2
3754 d_t_work(i)=ddt1+ddt2
3759 d_t(j,0)=d_t_work(j)
3764 dc(j,i)=dc_work(ind+j)
3765 d_t(j,i)=d_t_work(ind+j)
3771 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3772 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3773 .and.(mnum.lt.4)) then
3776 dc(j,inres)=dc_work(ind+j)
3777 d_t(j,inres)=d_t_work(ind+j)
3783 end subroutine sd_verlet1_ciccotti
3784 !-----------------------------------------------------------------------------
3785 subroutine sd_verlet2_ciccotti
3787 ! Calculating the adjusted velocities for accelerations
3789 ! implicit real*8 (a-h,o-z)
3790 ! include 'DIMENSIONS'
3791 ! include 'COMMON.CONTROL'
3792 ! include 'COMMON.VAR'
3793 ! include 'COMMON.MD'
3795 ! include 'COMMON.LANGEVIN'
3797 ! include 'COMMON.LANGEVIN.lang0'
3799 ! include 'COMMON.CHAIN'
3800 ! include 'COMMON.DERIV'
3801 ! include 'COMMON.GEO'
3802 ! include 'COMMON.LOCAL'
3803 ! include 'COMMON.INTERACT'
3804 ! include 'COMMON.IOUNITS'
3805 ! include 'COMMON.NAMES'
3806 !el real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV !(MAXRES6) maxres6=6*maxres
3807 real(kind=8),dimension(6*nres) :: stochforcvecV !(MAXRES6) maxres6=6*maxres
3808 !el common /stochcalc/ stochforcvec
3809 real(kind=8) :: ddt1,ddt2
3810 integer :: i,j,ind,inres
3812 ! Compute the stochastic forces which contribute to velocity change
3814 call stochastic_force(stochforcvecV)
3821 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3822 ! ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
3823 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
3825 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
3829 d_t(j,0)=d_t_work(j)
3834 d_t(j,i)=d_t_work(ind+j)
3840 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3842 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3845 d_t(j,inres)=d_t_work(ind+j)
3851 end subroutine sd_verlet2_ciccotti
3853 !-----------------------------------------------------------------------------
3855 !-----------------------------------------------------------------------------
3856 subroutine inertia_tensor
3858 ! Calculating the intertia tensor for the entire protein in order to
3859 ! remove the perpendicular components of velocity matrix which cause
3860 ! the molecule to rotate.
3863 ! implicit real*8 (a-h,o-z)
3864 ! include 'DIMENSIONS'
3865 ! include 'COMMON.CONTROL'
3866 ! include 'COMMON.VAR'
3867 ! include 'COMMON.MD'
3868 ! include 'COMMON.CHAIN'
3869 ! include 'COMMON.DERIV'
3870 ! include 'COMMON.GEO'
3871 ! include 'COMMON.LOCAL'
3872 ! include 'COMMON.INTERACT'
3873 ! include 'COMMON.IOUNITS'
3874 ! include 'COMMON.NAMES'
3876 real(kind=8),dimension(3,3) :: Im,Imcp,eigvec,Id
3877 real(kind=8),dimension(3) :: pr,eigval,L,vp,vrot
3878 real(kind=8) :: M_SC,mag,mag2,M_PEP
3879 real(kind=8),dimension(3,0:nres) :: vpp !(3,0:MAXRES)
3880 real(kind=8),dimension(3) :: vs_p,pp,incr,v
3881 real(kind=8),dimension(3,3) :: pr1,pr2
3883 !el common /gucio/ cm
3884 integer :: iti,inres,i,j,k,mnum
3895 ! calculating the center of the mass of the protein
3899 if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
3900 ! if (itype(i,mnum).eq.ntyp1_molec(mnum)) cycle
3901 M_PEP=M_PEP+mp(mnum)
3904 cm(j)=cm(j)+(c(j,i)+0.5d0*dc(j,i))*mp(mnum)
3913 if (mnum.ge.5) cycle
3914 iti=iabs(itype(i,mnum))
3915 M_SC=M_SC+msc(iabs(iti),mnum)
3918 cm(j)=cm(j)+msc(iabs(iti),mnum)*c(j,inres)
3922 cm(j)=cm(j)/(M_SC+M_PEP)
3927 if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
3929 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
3931 Im(1,1)=Im(1,1)+mp(mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
3932 Im(1,2)=Im(1,2)-mp(mnum)*pr(1)*pr(2)
3933 Im(1,3)=Im(1,3)-mp(mnum)*pr(1)*pr(3)
3934 Im(2,3)=Im(2,3)-mp(mnum)*pr(2)*pr(3)
3935 Im(2,2)=Im(2,2)+mp(mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
3936 Im(3,3)=Im(3,3)+mp(mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
3941 iti=iabs(itype(i,mnum))
3942 if (mnum.ge.5) cycle
3945 pr(j)=c(j,inres)-cm(j)
3947 Im(1,1)=Im(1,1)+msc(iabs(iti),mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
3948 Im(1,2)=Im(1,2)-msc(iabs(iti),mnum)*pr(1)*pr(2)
3949 Im(1,3)=Im(1,3)-msc(iabs(iti),mnum)*pr(1)*pr(3)
3950 Im(2,3)=Im(2,3)-msc(iabs(iti),mnum)*pr(2)*pr(3)
3951 Im(2,2)=Im(2,2)+msc(iabs(iti),mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
3952 Im(3,3)=Im(3,3)+msc(iabs(iti),mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
3957 Im(1,1)=Im(1,1)+Ip(mnum)*(1-dc_norm(1,i)*dc_norm(1,i))* &
3958 vbld(i+1)*vbld(i+1)*0.25d0
3959 Im(1,2)=Im(1,2)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(2,i))* &
3960 vbld(i+1)*vbld(i+1)*0.25d0
3961 Im(1,3)=Im(1,3)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(3,i))* &
3962 vbld(i+1)*vbld(i+1)*0.25d0
3963 Im(2,3)=Im(2,3)+Ip(mnum)*(-dc_norm(2,i)*dc_norm(3,i))* &
3964 vbld(i+1)*vbld(i+1)*0.25d0
3965 Im(2,2)=Im(2,2)+Ip(mnum)*(1-dc_norm(2,i)*dc_norm(2,i))* &
3966 vbld(i+1)*vbld(i+1)*0.25d0
3967 Im(3,3)=Im(3,3)+Ip(mnum)*(1-dc_norm(3,i)*dc_norm(3,i))* &
3968 vbld(i+1)*vbld(i+1)*0.25d0
3974 ! if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
3975 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3976 .and.(mnum.lt.4)) then
3977 iti=iabs(itype(i,mnum))
3979 Im(1,1)=Im(1,1)+Isc(iti,mnum)*(1-dc_norm(1,inres)* &
3980 dc_norm(1,inres))*vbld(inres)*vbld(inres)
3981 Im(1,2)=Im(1,2)-Isc(iti,mnum)*(dc_norm(1,inres)* &
3982 dc_norm(2,inres))*vbld(inres)*vbld(inres)
3983 Im(1,3)=Im(1,3)-Isc(iti,mnum)*(dc_norm(1,inres)* &
3984 dc_norm(3,inres))*vbld(inres)*vbld(inres)
3985 Im(2,3)=Im(2,3)-Isc(iti,mnum)*(dc_norm(2,inres)* &
3986 dc_norm(3,inres))*vbld(inres)*vbld(inres)
3987 Im(2,2)=Im(2,2)+Isc(iti,mnum)*(1-dc_norm(2,inres)* &
3988 dc_norm(2,inres))*vbld(inres)*vbld(inres)
3989 Im(3,3)=Im(3,3)+Isc(iti,mnum)*(1-dc_norm(3,inres)* &
3990 dc_norm(3,inres))*vbld(inres)*vbld(inres)
3995 ! write(iout,*) "The angular momentum before adjustment:"
3996 ! write(iout,*) (L(j),j=1,3)
4002 ! Copying the Im matrix for the djacob subroutine
4010 ! Finding the eigenvectors and eignvalues of the inertia tensor
4011 call djacob(3,3,10000,1.0d-10,Imcp,eigvec,eigval)
4012 ! write (iout,*) "Eigenvalues & Eigenvectors"
4013 ! write (iout,'(5x,3f10.5)') (eigval(i),i=1,3)
4016 ! write (iout,'(i5,3f10.5)') i,(eigvec(i,j),j=1,3)
4018 ! Constructing the diagonalized matrix
4020 if (dabs(eigval(i)).gt.1.0d-15) then
4021 Id(i,i)=1.0d0/eigval(i)
4028 Imcp(i,j)=eigvec(j,i)
4034 pr1(i,j)=pr1(i,j)+Id(i,k)*Imcp(k,j)
4041 pr2(i,j)=pr2(i,j)+eigvec(i,k)*pr1(k,j)
4045 ! Calculating the total rotational velocity of the molecule
4048 vrot(i)=vrot(i)+pr2(i,j)*L(j)
4051 ! Resetting the velocities
4053 call vecpr(vrot(1),dc(1,i),vp)
4055 ! print *,"HERE2",d_t(j,i),vp(j)
4056 d_t(j,i)=d_t(j,i)-vp(j)
4057 ! print *,"HERE2",d_t(j,i),vp(j)
4062 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4063 .and.(mnum.lt.4)) then
4064 ! if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
4065 ! if(itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
4067 call vecpr(vrot(1),dc(1,inres),vp)
4069 d_t(j,inres)=d_t(j,inres)-vp(j)
4074 ! write(iout,*) "The angular momentum after adjustment:"
4075 ! write(iout,*) (L(j),j=1,3)
4078 end subroutine inertia_tensor
4079 !-----------------------------------------------------------------------------
4080 subroutine angmom(cm,L)
4083 ! implicit real*8 (a-h,o-z)
4084 ! include 'DIMENSIONS'
4085 ! include 'COMMON.CONTROL'
4086 ! include 'COMMON.VAR'
4087 ! include 'COMMON.MD'
4088 ! include 'COMMON.CHAIN'
4089 ! include 'COMMON.DERIV'
4090 ! include 'COMMON.GEO'
4091 ! include 'COMMON.LOCAL'
4092 ! include 'COMMON.INTERACT'
4093 ! include 'COMMON.IOUNITS'
4094 ! include 'COMMON.NAMES'
4095 real(kind=8) :: mscab
4096 real(kind=8),dimension(3) :: L,cm,pr,vp,vrot,incr,v,pp
4097 integer :: iti,inres,i,j,mnum
4098 ! Calculate the angular momentum
4107 if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
4109 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
4112 v(j)=incr(j)+0.5d0*d_t(j,i)
4115 incr(j)=incr(j)+d_t(j,i)
4117 call vecpr(pr(1),v(1),vp)
4119 L(j)=L(j)+mp(mnum)*vp(j)
4120 ! print *,"HERE3",J,i,L(j),mp(mnum),Ip(mnum),mnum
4124 pp(j)=0.5d0*d_t(j,i)
4126 call vecpr(pr(1),pp(1),vp)
4129 L(j)=L(j)+Ip(mnum)*vp(j)
4137 iti=iabs(itype(i,mnum))
4145 pr(j)=c(j,inres)-cm(j)
4148 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4149 .and.(mnum.lt.4)) then
4151 v(j)=incr(j)+d_t(j,inres)
4158 ! print *,i,pr,"pr",v
4159 call vecpr(pr(1),v(1),vp)
4160 ! write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),&
4161 ! " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
4163 L(j)=L(j)+mscab*vp(j)
4165 ! write (iout,*) "L",(l(j),j=1,3)
4166 ! print *,"L",(l(j),j=1,3),i,vp(1)
4168 if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4169 .and.(mnum.lt.4)) then
4171 v(j)=incr(j)+d_t(j,inres)
4173 call vecpr(dc(1,inres),d_t(1,inres),vp)
4175 L(j)=L(j)+Isc(iti,mnum)*vp(j)
4179 incr(j)=incr(j)+d_t(j,i)
4183 end subroutine angmom
4184 !-----------------------------------------------------------------------------
4185 subroutine vcm_vel(vcm)
4188 ! implicit real*8 (a-h,o-z)
4189 ! include 'DIMENSIONS'
4190 ! include 'COMMON.VAR'
4191 ! include 'COMMON.MD'
4192 ! include 'COMMON.CHAIN'
4193 ! include 'COMMON.DERIV'
4194 ! include 'COMMON.GEO'
4195 ! include 'COMMON.LOCAL'
4196 ! include 'COMMON.INTERACT'
4197 ! include 'COMMON.IOUNITS'
4198 real(kind=8),dimension(3) :: vcm,vv
4199 real(kind=8) :: summas,amas
4209 if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
4211 summas=summas+mp(mnum)
4213 vcm(j)=vcm(j)+mp(mnum)*(vv(j)+0.5d0*d_t(j,i))
4214 ! print *,i,j,vv(j),d_t(j,i)
4218 amas=msc(iabs(itype(i,mnum)),mnum)
4223 if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4224 .and.(mnum.lt.4)) then
4226 vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres))
4230 vcm(j)=vcm(j)+amas*vv(j)
4234 vv(j)=vv(j)+d_t(j,i)
4237 ! write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas
4239 vcm(j)=vcm(j)/summas
4242 end subroutine vcm_vel
4243 !-----------------------------------------------------------------------------
4245 !-----------------------------------------------------------------------------
4247 ! RATTLE algorithm for velocity Verlet - step 1, UNRES
4251 ! implicit real*8 (a-h,o-z)
4252 ! include 'DIMENSIONS'
4254 ! include 'COMMON.CONTROL'
4255 ! include 'COMMON.VAR'
4256 ! include 'COMMON.MD'
4258 ! include 'COMMON.LANGEVIN'
4260 ! include 'COMMON.LANGEVIN.lang0'
4262 ! include 'COMMON.CHAIN'
4263 ! include 'COMMON.DERIV'
4264 ! include 'COMMON.GEO'
4265 ! include 'COMMON.LOCAL'
4266 ! include 'COMMON.INTERACT'
4267 ! include 'COMMON.IOUNITS'
4268 ! include 'COMMON.NAMES'
4269 ! include 'COMMON.TIME1'
4270 !el real(kind=8) :: gginv(2*nres,2*nres),&
4271 !el gdc(3,2*nres,2*nres)
4272 real(kind=8) :: dC_uncor(3,2*nres) !,&
4273 !el real(kind=8) :: Cmat(2*nres,2*nres)
4274 real(kind=8) :: x(2*nres),xcorr(3,2*nres) !maxres2=2*maxres
4275 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
4276 !el common /przechowalnia/ nbond
4277 integer :: max_rattle = 5
4278 logical :: lprn = .false., lprn1 = .false., not_done
4279 real(kind=8) :: tol_rattle = 1.0d-5
4281 integer :: ii,i,j,jj,l,ind,ind1,nres2
4284 !el /common/ przechowalnia
4286 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4287 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4288 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4290 if (lprn) write (iout,*) "RATTLE1"
4294 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4295 .and.(mnum.lt.4)) nbond=nbond+1
4297 ! Make a folded form of the Ginv-matrix
4310 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4315 if (itype(k,1).ne.10 .and. itype(k,mnum).ne.ntyp1_molec(mnum)&
4316 .and.(mnum.lt.4)) then
4320 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4328 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4339 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4343 if (itype(k,1).ne.10) then
4347 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4355 write (iout,*) "Matrix GGinv"
4356 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
4362 if (iter.gt.max_rattle) then
4363 write (iout,*) "Error - too many iterations in RATTLE."
4366 ! Calculate the matrix C = GG**(-1) dC_old o dC
4371 dC_uncor(j,ind1)=dC(j,i)
4375 if (itype(i,1).ne.10) then
4378 dC_uncor(j,ind1)=dC(j,i+nres)
4387 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
4391 if (itype(k,1).ne.10) then
4394 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
4399 ! Calculate deviations from standard virtual-bond lengths
4403 x(ind)=vbld(i+1)**2-vbl**2
4406 if (itype(i,1).ne.10) then
4408 x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4412 write (iout,*) "Coordinates and violations"
4414 write(iout,'(i5,3f10.5,5x,e15.5)') &
4415 i,(dC_uncor(j,i),j=1,3),x(i)
4417 write (iout,*) "Velocities and violations"
4421 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4422 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4426 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4427 .and.(mnum.lt.4)) then
4430 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4431 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4432 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4435 ! write (iout,*) "gdc"
4437 ! write (iout,*) "i",i
4439 ! write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4445 if (dabs(x(i)).gt.xmax) then
4449 if (xmax.lt.tol_rattle) then
4453 ! Calculate the matrix of the system of equations
4458 Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
4463 write (iout,*) "Matrix Cmat"
4464 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4466 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4467 ! Add constraint term to positions
4474 xx = xx+x(ii)*gdc(j,ind,ii)
4478 d_t_new(j,i)=d_t_new(j,i)-xx/d_time
4482 if (itype(i,1).ne.10) then
4487 xx = xx+x(ii)*gdc(j,ind,ii)
4490 dC(j,i+nres)=dC(j,i+nres)-xx
4491 d_t_new(j,i+nres)=d_t_new(j,i+nres)-xx/d_time
4495 ! Rebuild the chain using the new coordinates
4496 call chainbuild_cart
4498 write (iout,*) "New coordinates, Lagrange multipliers,",&
4499 " and differences between actual and standard bond lengths"
4503 xx=vbld(i+1)**2-vbl**2
4504 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4505 i,(dC(j,i),j=1,3),x(ind),xx
4509 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4512 xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4513 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4514 i,(dC(j,i+nres),j=1,3),x(ind),xx
4517 write (iout,*) "Velocities and violations"
4521 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4522 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4525 if (itype(i,1).ne.10) then
4527 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4528 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4529 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4536 10 write (iout,*) "Error - singularity in solving the system",&
4537 " of equations for Lagrange multipliers."
4541 "RATTLE inactive; use -DRATTLE switch at compile time."
4544 end subroutine rattle1
4545 !-----------------------------------------------------------------------------
4547 ! RATTLE algorithm for velocity Verlet - step 2, UNRES
4551 ! implicit real*8 (a-h,o-z)
4552 ! include 'DIMENSIONS'
4554 ! include 'COMMON.CONTROL'
4555 ! include 'COMMON.VAR'
4556 ! include 'COMMON.MD'
4558 ! include 'COMMON.LANGEVIN'
4560 ! include 'COMMON.LANGEVIN.lang0'
4562 ! include 'COMMON.CHAIN'
4563 ! include 'COMMON.DERIV'
4564 ! include 'COMMON.GEO'
4565 ! include 'COMMON.LOCAL'
4566 ! include 'COMMON.INTERACT'
4567 ! include 'COMMON.IOUNITS'
4568 ! include 'COMMON.NAMES'
4569 ! include 'COMMON.TIME1'
4570 !el real(kind=8) :: gginv(2*nres,2*nres),&
4571 !el gdc(3,2*nres,2*nres)
4572 real(kind=8) :: dC_uncor(3,2*nres) !,&
4573 !el Cmat(2*nres,2*nres)
4574 real(kind=8) :: x(2*nres) !maxres2=2*maxres
4575 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
4576 !el common /przechowalnia/ nbond
4577 integer :: max_rattle = 5
4578 logical :: lprn = .false., lprn1 = .false., not_done
4579 real(kind=8) :: tol_rattle = 1.0d-5
4583 !el /common/ przechowalnia
4584 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4585 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4586 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4588 if (lprn) write (iout,*) "RATTLE2"
4589 if (lprn) write (iout,*) "Velocity correction"
4590 ! Calculate the matrix G dC
4596 gdc(j,i,ind)=GGinv(i,ind)*dC(j,k)
4601 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4602 .and.(mnum.lt.4)) then
4605 gdc(j,i,ind)=GGinv(i,ind)*dC(j,k+nres)
4611 ! write (iout,*) "gdc"
4613 ! write (iout,*) "i",i
4615 ! write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4619 ! Calculate the matrix of the system of equations
4626 Cmat(ind,j)=Cmat(ind,j)+dC(k,i)*gdc(k,ind,j)
4632 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4633 .and.(mnum.lt.4)) then
4638 Cmat(ind,j)=Cmat(ind,j)+dC(k,i+nres)*gdc(k,ind,j)
4643 ! Calculate the scalar product dC o d_t_new
4647 x(ind)=scalar(d_t(1,i),dC(1,i))
4651 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4652 .and.(mnum.lt.4)) then
4654 x(ind)=scalar(d_t(1,i+nres),dC(1,i+nres))
4658 write (iout,*) "Velocities and violations"
4662 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4663 i,ind,(d_t(j,i),j=1,3),x(ind)
4667 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4668 .and.(mnum.lt.4)) then
4670 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4671 i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind)
4677 if (dabs(x(i)).gt.xmax) then
4681 if (xmax.lt.tol_rattle) then
4686 write (iout,*) "Matrix Cmat"
4687 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4689 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4690 ! Add constraint term to velocities
4697 xx = xx+x(ii)*gdc(j,ind,ii)
4699 d_t(j,i)=d_t(j,i)-xx
4704 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4705 .and.(mnum.lt.4)) then
4710 xx = xx+x(ii)*gdc(j,ind,ii)
4712 d_t(j,i+nres)=d_t(j,i+nres)-xx
4718 "New velocities, Lagrange multipliers violations"
4722 if (lprn) write (iout,'(2i5,3f10.5,5x,2e15.5)') &
4723 i,ind,(d_t(j,i),j=1,3),x(ind),scalar(d_t(1,i),dC(1,i))
4727 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4730 write (iout,'(2i5,3f10.5,5x,2e15.5)') &
4731 i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind),&
4732 scalar(d_t(1,i+nres),dC(1,i+nres))
4738 10 write (iout,*) "Error - singularity in solving the system",&
4739 " of equations for Lagrange multipliers."
4743 "RATTLE inactive; use -DRATTLE option at compile time."
4746 end subroutine rattle2
4747 !-----------------------------------------------------------------------------
4748 subroutine rattle_brown
4749 ! RATTLE/LINCS algorithm for Brownian dynamics, UNRES
4753 ! implicit real*8 (a-h,o-z)
4754 ! include 'DIMENSIONS'
4756 ! include 'COMMON.CONTROL'
4757 ! include 'COMMON.VAR'
4758 ! include 'COMMON.MD'
4760 ! include 'COMMON.LANGEVIN'
4762 ! include 'COMMON.LANGEVIN.lang0'
4764 ! include 'COMMON.CHAIN'
4765 ! include 'COMMON.DERIV'
4766 ! include 'COMMON.GEO'
4767 ! include 'COMMON.LOCAL'
4768 ! include 'COMMON.INTERACT'
4769 ! include 'COMMON.IOUNITS'
4770 ! include 'COMMON.NAMES'
4771 ! include 'COMMON.TIME1'
4772 !el real(kind=8) :: gginv(2*nres,2*nres),&
4773 !el gdc(3,2*nres,2*nres)
4774 real(kind=8) :: dC_uncor(3,2*nres) !,&
4775 !el real(kind=8) :: Cmat(2*nres,2*nres)
4776 real(kind=8) :: x(2*nres) !maxres2=2*maxres
4777 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
4778 !el common /przechowalnia/ nbond
4779 integer :: max_rattle = 5
4780 logical :: lprn = .false., lprn1 = .false., not_done
4781 real(kind=8) :: tol_rattle = 1.0d-5
4785 !el /common/ przechowalnia
4786 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4787 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4788 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4791 if (lprn) write (iout,*) "RATTLE_BROWN"
4794 if (itype(i,1).ne.10) nbond=nbond+1
4796 ! Make a folded form of the Ginv-matrix
4809 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4813 if (itype(k,1).ne.10) then
4817 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4824 if (itype(i,1).ne.10) then
4834 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4838 if (itype(k,1).ne.10) then
4842 if (j.eq.1 .and. l.eq.1)GGinv(ii,jj)=fricmat(ind,ind1)
4850 write (iout,*) "Matrix GGinv"
4851 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
4857 if (iter.gt.max_rattle) then
4858 write (iout,*) "Error - too many iterations in RATTLE."
4861 ! Calculate the matrix C = GG**(-1) dC_old o dC
4866 dC_uncor(j,ind1)=dC(j,i)
4870 if (itype(i,1).ne.10) then
4873 dC_uncor(j,ind1)=dC(j,i+nres)
4882 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
4886 if (itype(k,1).ne.10) then
4889 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
4894 ! Calculate deviations from standard virtual-bond lengths
4898 x(ind)=vbld(i+1)**2-vbl**2
4901 if (itype(i,1).ne.10) then
4903 x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4907 write (iout,*) "Coordinates and violations"
4909 write(iout,'(i5,3f10.5,5x,e15.5)') &
4910 i,(dC_uncor(j,i),j=1,3),x(i)
4912 write (iout,*) "Velocities and violations"
4916 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4917 i,ind,(d_t(j,i),j=1,3),scalar(d_t(1,i),dC_old(1,i))
4920 if (itype(i,1).ne.10) then
4922 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4923 i+nres,ind,(d_t(j,i+nres),j=1,3),&
4924 scalar(d_t(1,i+nres),dC_old(1,i+nres))
4927 write (iout,*) "gdc"
4929 write (iout,*) "i",i
4931 write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4937 if (dabs(x(i)).gt.xmax) then
4941 if (xmax.lt.tol_rattle) then
4945 ! Calculate the matrix of the system of equations
4950 Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
4955 write (iout,*) "Matrix Cmat"
4956 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4958 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4959 ! Add constraint term to positions
4966 xx = xx+x(ii)*gdc(j,ind,ii)
4969 d_t(j,i)=d_t(j,i)+xx/d_time
4974 if (itype(i,1).ne.10) then
4979 xx = xx+x(ii)*gdc(j,ind,ii)
4982 d_t(j,i+nres)=d_t(j,i+nres)+xx/d_time
4983 dC(j,i+nres)=dC(j,i+nres)+xx
4987 ! Rebuild the chain using the new coordinates
4988 call chainbuild_cart
4990 write (iout,*) "New coordinates, Lagrange multipliers,",&
4991 " and differences between actual and standard bond lengths"
4995 xx=vbld(i+1)**2-vbl**2
4996 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4997 i,(dC(j,i),j=1,3),x(ind),xx
5000 if (itype(i,1).ne.10) then
5002 xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
5003 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
5004 i,(dC(j,i+nres),j=1,3),x(ind),xx
5007 write (iout,*) "Velocities and violations"
5011 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5012 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
5015 if (itype(i,1).ne.10) then
5017 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5018 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
5019 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
5026 10 write (iout,*) "Error - singularity in solving the system",&
5027 " of equations for Lagrange multipliers."
5031 "RATTLE inactive; use -DRATTLE option at compile time"
5034 end subroutine rattle_brown
5035 !-----------------------------------------------------------------------------
5037 !-----------------------------------------------------------------------------
5038 subroutine friction_force
5043 ! implicit real*8 (a-h,o-z)
5044 ! include 'DIMENSIONS'
5045 ! include 'COMMON.VAR'
5046 ! include 'COMMON.CHAIN'
5047 ! include 'COMMON.DERIV'
5048 ! include 'COMMON.GEO'
5049 ! include 'COMMON.LOCAL'
5050 ! include 'COMMON.INTERACT'
5051 ! include 'COMMON.MD'
5053 ! include 'COMMON.LANGEVIN'
5055 ! include 'COMMON.LANGEVIN.lang0'
5057 ! include 'COMMON.IOUNITS'
5058 !el real(kind=8),dimension(6*nres) :: gamvec !(MAXRES6) maxres6=6*maxres
5059 !el common /syfek/ gamvec
5060 real(kind=8) :: vv(3),vvtot(3,nres),v_work(6*nres) !,&
5061 !el ginvfric(2*nres,2*nres) !maxres2=2*maxres
5062 !el common /przechowalnia/ ginvfric
5064 logical :: lprn, checkmode
5065 integer :: i,j,ind,k,nres2,nres6,mnum
5070 ! if (large) lprn=.true.
5071 ! if (large) checkmode=.true.
5072 if(.not.allocated(gamvec)) allocate(gamvec(nres6)) !(MAXRES6)
5073 if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
5081 d_t_work(j)=d_t(j,0)
5086 d_t_work(ind+j)=d_t(j,i)
5092 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5093 .and.(mnum.lt.4)) then
5095 d_t_work(ind+j)=d_t(j,i+nres)
5101 call fricmat_mult(d_t_work,fric_work)
5103 if (.not.checkmode) return
5106 write (iout,*) "d_t_work and fric_work"
5108 write (iout,'(i3,2e15.5)') i,d_t_work(i),fric_work(i)
5112 friction(j,0)=fric_work(j)
5117 friction(j,i)=fric_work(ind+j)
5123 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5124 .and.(mnum.lt.4)) then
5125 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
5127 friction(j,i+nres)=fric_work(ind+j)
5133 write(iout,*) "Friction backbone"
5135 write(iout,'(i5,3e15.5,5x,3e15.5)') &
5136 i,(friction(j,i),j=1,3),(d_t(j,i),j=1,3)
5138 write(iout,*) "Friction side chain"
5140 write(iout,'(i5,3e15.5,5x,3e15.5)') &
5141 i,(friction(j,i+nres),j=1,3),(d_t(j,i+nres),j=1,3)
5150 vvtot(j,i)=vv(j)+0.5d0*d_t(j,i)
5151 vvtot(j,i+nres)=vv(j)+d_t(j,i+nres)
5152 vv(j)=vv(j)+d_t(j,i)
5155 write (iout,*) "vvtot backbone and sidechain"
5157 write (iout,'(i5,3e15.5,5x,3e15.5)') i,(vvtot(j,i),j=1,3),&
5158 (vvtot(j,i+nres),j=1,3)
5163 v_work(ind+j)=vvtot(j,i)
5169 v_work(ind+j)=vvtot(j,i+nres)
5173 write (iout,*) "v_work gamvec and site-based friction forces"
5175 write (iout,'(i5,3e15.5)') i,v_work(i),gamvec(i),&
5179 ! fric_work1(i)=0.0d0
5181 ! fric_work1(i)=fric_work1(i)-A(j,i)*gamvec(j)*v_work(j)
5184 ! write (iout,*) "fric_work and fric_work1"
5186 ! write (iout,'(i5,2e15.5)') i,fric_work(i),fric_work1(i)
5192 ginvfric(i,j)=ginvfric(i,j)+ginv(i,k)*fricmat(k,j)
5196 write (iout,*) "ginvfric"
5198 write (iout,'(i5,100f8.3)') i,(ginvfric(i,j),j=1,dimen)
5200 write (iout,*) "symmetry check"
5203 write (iout,*) i,j,ginvfric(i,j)-ginvfric(j,i)
5208 end subroutine friction_force
5209 !-----------------------------------------------------------------------------
5210 subroutine setup_fricmat
5214 use control_data, only:time_Bcast
5215 use control, only:tcpu
5217 ! implicit real*8 (a-h,o-z)
5221 real(kind=8) :: time00
5223 ! include 'DIMENSIONS'
5224 ! include 'COMMON.VAR'
5225 ! include 'COMMON.CHAIN'
5226 ! include 'COMMON.DERIV'
5227 ! include 'COMMON.GEO'
5228 ! include 'COMMON.LOCAL'
5229 ! include 'COMMON.INTERACT'
5230 ! include 'COMMON.MD'
5231 ! include 'COMMON.SETUP'
5232 ! include 'COMMON.TIME1'
5233 ! integer licznik /0/
5236 ! include 'COMMON.LANGEVIN'
5238 ! include 'COMMON.LANGEVIN.lang0'
5240 ! include 'COMMON.IOUNITS'
5242 integer :: i,j,ind,ind1,m
5243 logical :: lprn = .false.
5244 real(kind=8) :: dtdi !el ,gamvec(2*nres)
5245 !el real(kind=8),dimension(2*nres,2*nres) :: ginvfric,fcopy
5246 ! real(kind=8),allocatable,dimension(:,:) :: fcopy
5247 !el real(kind=8),dimension(2*nres*(2*nres+1)/2) :: Ghalf !(mmaxres2) (mmaxres2=(maxres2*(maxres2+1)/2))
5248 !el common /syfek/ gamvec
5249 real(kind=8) :: work(8*2*nres)
5250 integer :: iwork(2*nres)
5251 !el common /przechowalnia/ ginvfric,Ghalf,fcopy
5252 integer :: ii,iti,k,l,nzero,nres2,nres6,ierr,mnum
5256 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
5257 if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
5258 if (fg_rank.ne.king) goto 10
5263 if(.not.allocated(gamvec)) allocate(gamvec(nres2)) !(MAXRES2)
5264 if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
5265 if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
5266 !el allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
5267 if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !maxres2=2*maxres
5269 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
5270 ! Zeroing out fricmat
5276 ! Load the friction coefficients corresponding to peptide groups
5281 gamvec(ind1)=gamp(mnum)
5284 if (molnum(nct).eq.5) then
5287 gamvec(ind1)=gamp(mnum)
5289 ! Load the friction coefficients corresponding to side chains
5293 gamsc(ntyp1_molec(j),j)=1.0d0
5300 gamvec(ii)=gamsc(iabs(iti),mnum)
5302 if (surfarea) call sdarea(gamvec)
5304 ! write (iout,*) "Matrix A and vector gamma"
5306 ! write (iout,'(i2,$)') i
5308 ! write (iout,'(f4.1,$)') A(i,j)
5310 ! write (iout,'(f8.3)') gamvec(i)
5314 write (iout,*) "Vector gamvec"
5316 write (iout,'(i5,f10.5)') i, gamvec(i)
5320 ! The friction matrix
5325 dtdi=dtdi+A(j,k)*A(j,i)*gamvec(j)
5332 write (iout,'(//a)') "Matrix fricmat"
5333 call matout2(dimen,dimen,nres2,nres2,fricmat)
5335 if (lang.eq.2 .or. lang.eq.3) then
5336 ! Mass-scale the friction matrix if non-direct integration will be performed
5342 Ginvfric(i,j)=Ginvfric(i,j)+ &
5343 Gsqrm(i,k)*Gsqrm(l,j)*fricmat(k,l)
5348 ! Diagonalize the friction matrix
5353 Ghalf(ind)=Ginvfric(i,j)
5356 call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
5359 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
5360 " mass-scaled friction matrix"
5361 call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
5363 ! Precompute matrices for tinker stochastic integrator
5370 mt1(i,j)=mt1(i,j)+fricvec(k,i)*gsqrm(k,j)
5371 mt2(i,j)=mt2(i,j)+fricvec(k,i)*gsqrp(k,j)
5377 else if (lang.eq.4) then
5378 ! Diagonalize the friction matrix
5383 Ghalf(ind)=fricmat(i,j)
5386 call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
5389 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
5391 call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
5393 ! Determine the number of zero eigenvalues of the friction matrix
5394 nzero=max0(dimen-dimen1,0)
5395 ! do while (fricgam(nzero+1).le.1.0d-5 .and. nzero.lt.dimen)
5398 write (iout,*) "Number of zero eigenvalues:",nzero
5403 fricmat(i,j)=fricmat(i,j) &
5404 +fricvec(i,k)*fricvec(j,k)/fricgam(k)
5409 write (iout,'(//a)') "Generalized inverse of fricmat"
5410 call matout(dimen,dimen,nres6,nres6,fricmat)
5415 if (nfgtasks.gt.1) then
5416 if (fg_rank.eq.0) then
5417 ! The matching BROADCAST for fg processors is called in ERGASTULUM
5423 call MPI_Bcast(10,1,MPI_INTEGER,king,FG_COMM,IERROR)
5425 time_Bcast=time_Bcast+MPI_Wtime()-time00
5427 time_Bcast=time_Bcast+tcpu()-time00
5429 ! print *,"Processor",myrank,
5430 ! & " BROADCAST iorder in SETUP_FRICMAT"
5433 write (iout,*) "setup_fricmat licznik"!,licznik !sp
5439 ! Scatter the friction matrix
5440 call MPI_Scatterv(fricmat(1,1),nginv_counts(0),&
5441 nginv_start(0),MPI_DOUBLE_PRECISION,fcopy(1,1),&
5442 myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
5445 time_scatter=time_scatter+MPI_Wtime()-time00
5446 time_scatter_fmat=time_scatter_fmat+MPI_Wtime()-time00
5448 time_scatter=time_scatter+tcpu()-time00
5449 time_scatter_fmat=time_scatter_fmat+tcpu()-time00
5453 do j=1,2*my_ng_count
5454 fricmat(j,i)=fcopy(i,j)
5457 ! write (iout,*) "My chunk of fricmat"
5458 ! call MATOUT2(my_ng_count,dimen,maxres2,maxres2,fcopy)
5462 end subroutine setup_fricmat
5463 !-----------------------------------------------------------------------------
5464 subroutine stochastic_force(stochforcvec)
5467 use random, only:anorm_distr
5468 ! implicit real*8 (a-h,o-z)
5469 ! include 'DIMENSIONS'
5470 use control, only: tcpu
5475 ! include 'COMMON.VAR'
5476 ! include 'COMMON.CHAIN'
5477 ! include 'COMMON.DERIV'
5478 ! include 'COMMON.GEO'
5479 ! include 'COMMON.LOCAL'
5480 ! include 'COMMON.INTERACT'
5481 ! include 'COMMON.MD'
5482 ! include 'COMMON.TIME1'
5484 ! include 'COMMON.LANGEVIN'
5486 ! include 'COMMON.LANGEVIN.lang0'
5488 ! include 'COMMON.IOUNITS'
5490 real(kind=8) :: x,sig,lowb,highb
5491 real(kind=8) :: ff(3),force(3,0:2*nres),zeta2,lowb2
5492 real(kind=8) :: highb2,sig2,forcvec(6*nres),stochforcvec(6*nres)
5493 real(kind=8) :: time00
5494 logical :: lprn = .false.
5495 integer :: i,j,ind,mnum
5499 stochforc(j,i)=0.0d0
5509 ! Compute the stochastic forces acting on bodies. Store in force.
5515 force(j,i)=anorm_distr(x,sig,lowb,highb)
5523 force(j,i+nres)=anorm_distr(x,sig2,lowb2,highb2)
5527 time_fsample=time_fsample+MPI_Wtime()-time00
5529 time_fsample=time_fsample+tcpu()-time00
5531 ! Compute the stochastic forces acting on virtual-bond vectors.
5537 stochforc(j,i)=ff(j)+0.5d0*force(j,i)
5540 ff(j)=ff(j)+force(j,i)
5542 ! if (itype(i+1,1).ne.ntyp1) then
5544 if (itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
5546 stochforc(j,i)=stochforc(j,i)+force(j,i+nres+1)
5547 ff(j)=ff(j)+force(j,i+nres+1)
5552 stochforc(j,0)=ff(j)+force(j,nnt+nres)
5556 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5557 .and.(mnum.lt.4)) then
5558 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
5560 stochforc(j,i+nres)=force(j,i+nres)
5566 stochforcvec(j)=stochforc(j,0)
5571 stochforcvec(ind+j)=stochforc(j,i)
5577 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5578 .and.(mnum.lt.4)) then
5579 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
5581 stochforcvec(ind+j)=stochforc(j,i+nres)
5587 write (iout,*) "stochforcvec"
5589 write(iout,'(i5,e15.5)') i,stochforcvec(i)
5591 write(iout,*) "Stochastic forces backbone"
5593 write(iout,'(i5,3e15.5)') i,(stochforc(j,i),j=1,3)
5595 write(iout,*) "Stochastic forces side chain"
5597 write(iout,'(i5,3e15.5)') &
5598 i,(stochforc(j,i+nres),j=1,3)
5606 write (iout,*) i,ind
5608 forcvec(ind+j)=force(j,i)
5613 write (iout,*) i,ind
5615 forcvec(j+ind)=force(j,i+nres)
5620 write (iout,*) "forcvec"
5624 write (iout,'(2i3,2f10.5)') i,j,force(j,i),&
5631 write (iout,'(2i3,2f10.5)') i,j,force(j,i+nres),&
5640 end subroutine stochastic_force
5641 !-----------------------------------------------------------------------------
5642 subroutine sdarea(gamvec)
5644 ! Scale the friction coefficients according to solvent accessible surface areas
5645 ! Code adapted from TINKER
5649 ! implicit real*8 (a-h,o-z)
5650 ! include 'DIMENSIONS'
5651 ! include 'COMMON.CONTROL'
5652 ! include 'COMMON.VAR'
5653 ! include 'COMMON.MD'
5655 ! include 'COMMON.LANGEVIN'
5657 ! include 'COMMON.LANGEVIN.lang0'
5659 ! include 'COMMON.CHAIN'
5660 ! include 'COMMON.DERIV'
5661 ! include 'COMMON.GEO'
5662 ! include 'COMMON.LOCAL'
5663 ! include 'COMMON.INTERACT'
5664 ! include 'COMMON.IOUNITS'
5665 ! include 'COMMON.NAMES'
5666 real(kind=8),dimension(2*nres) :: radius,gamvec !(maxres2)
5667 real(kind=8),parameter :: twosix = 1.122462048309372981d0
5668 logical :: lprn = .false.
5669 real(kind=8) :: probe,area,ratio
5670 integer :: i,j,ind,iti,mnum
5672 ! determine new friction coefficients every few SD steps
5674 ! set the atomic radii to estimates of sigma values
5676 ! print *,"Entered sdarea"
5682 ! Load peptide group radii
5685 radius(i)=pstok(mnum)
5687 ! Load side chain radii
5691 radius(i+nres)=restok(iti,mnum)
5694 ! write (iout,*) "i",i," radius",radius(i)
5697 radius(i) = radius(i) / twosix
5698 if (radius(i) .ne. 0.0d0) radius(i) = radius(i) + probe
5701 ! scale atomic friction coefficients by accessible area
5703 if (lprn) write (iout,*) &
5704 "Original gammas, surface areas, scaling factors, new gammas, ",&
5705 "std's of stochastic forces"
5708 if (radius(i).gt.0.0d0) then
5709 call surfatom (i,area,radius)
5710 ratio = dmax1(area/(4.0d0*pi*radius(i)**2),1.0d-1)
5711 if (lprn) write (iout,'(i5,3f10.5,$)') &
5712 i,gamvec(ind+1),area,ratio
5715 gamvec(ind) = ratio * gamvec(ind)
5718 stdforcp(i)=stdfp(mnum)*dsqrt(gamvec(ind))
5719 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcp(i)
5723 if (radius(i+nres).gt.0.0d0) then
5724 call surfatom (i+nres,area,radius)
5725 ratio = dmax1(area/(4.0d0*pi*radius(i+nres)**2),1.0d-1)
5726 if (lprn) write (iout,'(i5,3f10.5,$)') &
5727 i,gamvec(ind+1),area,ratio
5730 gamvec(ind) = ratio * gamvec(ind)
5733 stdforcsc(i)=stdfsc(itype(i,mnum),mnum)*dsqrt(gamvec(ind))
5734 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcsc(i)
5739 end subroutine sdarea
5740 !-----------------------------------------------------------------------------
5742 !-----------------------------------------------------------------------------
5745 ! ###################################################
5746 ! ## COPYRIGHT (C) 1996 by Jay William Ponder ##
5747 ! ## All Rights Reserved ##
5748 ! ###################################################
5750 ! ################################################################
5752 ! ## subroutine surfatom -- exposed surface area of an atom ##
5754 ! ################################################################
5757 ! "surfatom" performs an analytical computation of the surface
5758 ! area of a specified atom; a simplified version of "surface"
5760 ! literature references:
5762 ! T. J. Richmond, "Solvent Accessible Surface Area and
5763 ! Excluded Volume in Proteins", Journal of Molecular Biology,
5766 ! L. Wesson and D. Eisenberg, "Atomic Solvation Parameters
5767 ! Applied to Molecular Dynamics of Proteins in Solution",
5768 ! Protein Science, 1, 227-235 (1992)
5770 ! variables and parameters:
5772 ! ir number of atom for which area is desired
5773 ! area accessible surface area of the atom
5774 ! radius radii of each of the individual atoms
5777 subroutine surfatom(ir,area,radius)
5779 ! implicit real*8 (a-h,o-z)
5780 ! include 'DIMENSIONS'
5782 ! include 'COMMON.GEO'
5783 ! include 'COMMON.IOUNITS'
5785 integer :: nsup,nstart_sup
5786 ! double precision c,dc,dc_old,d_c_work,xloc,xrot,dc_norm
5787 ! common /chain/ c(3,maxres2+2),dc(3,0:maxres2),dc_old(3,0:maxres2),
5788 ! & xloc(3,maxres),xrot(3,maxres),dc_norm(3,0:maxres2),
5789 ! & dc_work(MAXRES6),nres,nres0
5790 integer,parameter :: maxarc=300
5794 integer :: mi,ni,narc
5795 integer :: key(maxarc)
5796 integer :: intag(maxarc)
5797 integer :: intag1(maxarc)
5798 real(kind=8) :: area,arcsum
5799 real(kind=8) :: arclen,exang
5800 real(kind=8) :: delta,delta2
5801 real(kind=8) :: eps,rmove
5802 real(kind=8) :: xr,yr,zr
5803 real(kind=8) :: rr,rrsq
5804 real(kind=8) :: rplus,rminus
5805 real(kind=8) :: axx,axy,axz
5806 real(kind=8) :: ayx,ayy
5807 real(kind=8) :: azx,azy,azz
5808 real(kind=8) :: uxj,uyj,uzj
5809 real(kind=8) :: tx,ty,tz
5810 real(kind=8) :: txb,tyb,td
5811 real(kind=8) :: tr2,tr,txr,tyr
5812 real(kind=8) :: tk1,tk2
5813 real(kind=8) :: thec,the,t,tb
5814 real(kind=8) :: txk,tyk,tzk
5815 real(kind=8) :: t1,ti,tf,tt
5816 real(kind=8) :: txj,tyj,tzj
5817 real(kind=8) :: ccsq,cc,xysq
5818 real(kind=8) :: bsqk,bk,cosine
5819 real(kind=8) :: dsqj,gi,pix2
5820 real(kind=8) :: therk,dk,gk
5821 real(kind=8) :: risqk,rik
5822 real(kind=8) :: radius(2*nres) !(maxatm) (maxatm=maxres2)
5823 real(kind=8) :: ri(maxarc),risq(maxarc)
5824 real(kind=8) :: ux(maxarc),uy(maxarc),uz(maxarc)
5825 real(kind=8) :: xc(maxarc),yc(maxarc),zc(maxarc)
5826 real(kind=8) :: xc1(maxarc),yc1(maxarc),zc1(maxarc)
5827 real(kind=8) :: dsq(maxarc),bsq(maxarc)
5828 real(kind=8) :: dsq1(maxarc),bsq1(maxarc)
5829 real(kind=8) :: arci(maxarc),arcf(maxarc)
5830 real(kind=8) :: ex(maxarc),lt(maxarc),gr(maxarc)
5831 real(kind=8) :: b(maxarc),b1(maxarc),bg(maxarc)
5832 real(kind=8) :: kent(maxarc),kout(maxarc)
5833 real(kind=8) :: ther(maxarc)
5834 logical :: moved,top
5835 logical :: omit(maxarc)
5838 ! maxatm = 2*nres !maxres2 maxres2=2*maxres
5839 ! maxlight = 8*maxatm
5842 ! maxtors = 4*maxatm
5844 ! zero out the surface area for the sphere of interest
5847 ! write (2,*) "ir",ir," radius",radius(ir)
5848 if (radius(ir) .eq. 0.0d0) return
5850 ! set the overlap significance and connectivity shift
5854 delta2 = delta * delta
5859 ! store coordinates and radius of the sphere of interest
5867 ! initialize values of some counters and summations
5876 ! test each sphere to see if it overlaps the sphere of interest
5879 if (i.eq.ir .or. radius(i).eq.0.0d0) goto 30
5880 rplus = rr + radius(i)
5882 if (abs(tx) .ge. rplus) goto 30
5884 if (abs(ty) .ge. rplus) goto 30
5886 if (abs(tz) .ge. rplus) goto 30
5888 ! check for sphere overlap by testing distance against radii
5890 xysq = tx*tx + ty*ty
5891 if (xysq .lt. delta2) then
5898 if (rplus-cc .le. delta) goto 30
5899 rminus = rr - radius(i)
5901 ! check to see if sphere of interest is completely buried
5903 if (cc-abs(rminus) .le. delta) then
5904 if (rminus .le. 0.0d0) goto 170
5908 ! check for too many overlaps with sphere of interest
5910 if (io .ge. maxarc) then
5912 20 format (/,' SURFATOM -- Increase the Value of MAXARC')
5916 ! get overlap between current sphere and sphere of interest
5925 gr(io) = (ccsq+rplus*rminus) / (2.0d0*rr*b1(io))
5931 ! case where no other spheres overlap the sphere of interest
5934 area = 4.0d0 * pi * rrsq
5938 ! case where only one sphere overlaps the sphere of interest
5941 area = pix2 * (1.0d0 + gr(1))
5942 area = mod(area,4.0d0*pi) * rrsq
5946 ! case where many spheres intersect the sphere of interest;
5947 ! sort the intersecting spheres by their degree of overlap
5949 call sort2 (io,gr,key)
5952 intag(i) = intag1(k)
5961 ! get radius of each overlap circle on surface of the sphere
5966 risq(i) = rrsq - gi*gi
5967 ri(i) = sqrt(risq(i))
5968 ther(i) = 0.5d0*pi - asin(min(1.0d0,max(-1.0d0,gr(i))))
5971 ! find boundary of inaccessible area on sphere of interest
5974 if (.not. omit(k)) then
5981 ! check to see if J circle is intersecting K circle;
5982 ! get distance between circle centers and sum of radii
5985 if (omit(j)) goto 60
5986 cc = (txk*xc(j)+tyk*yc(j)+tzk*zc(j))/(bk*b(j))
5987 cc = acos(min(1.0d0,max(-1.0d0,cc)))
5988 td = therk + ther(j)
5990 ! check to see if circles enclose separate regions
5992 if (cc .ge. td) goto 60
5994 ! check for circle J completely inside circle K
5996 if (cc+ther(j) .lt. therk) goto 40
5998 ! check for circles that are essentially parallel
6000 if (cc .gt. delta) goto 50
6005 ! check to see if sphere of interest is completely buried
6008 if (pix2-cc .le. td) goto 170
6014 ! find T value of circle intersections
6017 if (omit(k)) goto 110
6032 ! rotation matrix elements
6044 if (.not. omit(j)) then
6049 ! rotate spheres so K vector colinear with z-axis
6051 uxj = txj*axx + tyj*axy - tzj*axz
6052 uyj = tyj*ayy - txj*ayx
6053 uzj = txj*azx + tyj*azy + tzj*azz
6054 cosine = min(1.0d0,max(-1.0d0,uzj/b(j)))
6055 if (acos(cosine) .lt. therk+ther(j)) then
6056 dsqj = uxj*uxj + uyj*uyj
6061 tr2 = risqk*dsqj - tb*tb
6067 ! get T values of intersection for K circle
6070 tb = min(1.0d0,max(-1.0d0,tb))
6072 if (tyb-txr .lt. 0.0d0) tk1 = pix2 - tk1
6074 tb = min(1.0d0,max(-1.0d0,tb))
6076 if (tyb+txr .lt. 0.0d0) tk2 = pix2 - tk2
6077 thec = (rrsq*uzj-gk*bg(j)) / (rik*ri(j)*b(j))
6078 if (abs(thec) .lt. 1.0d0) then
6080 else if (thec .ge. 1.0d0) then
6082 else if (thec .le. -1.0d0) then
6086 ! see if "tk1" is entry or exit point; check t=0 point;
6087 ! "ti" is exit point, "tf" is entry point
6089 cosine = min(1.0d0,max(-1.0d0, &
6090 (uzj*gk-uxj*rik)/(b(j)*rr)))
6091 if ((acos(cosine)-ther(j))*(tk2-tk1) .le. 0.0d0) then
6099 if (narc .ge. maxarc) then
6101 70 format (/,' SURFATOM -- Increase the Value',&
6105 if (tf .le. ti) then
6126 ! special case; K circle without intersections
6128 if (narc .le. 0) goto 90
6130 ! general case; sum up arclength and set connectivity code
6132 call sort2 (narc,arci,key)
6137 if (narc .gt. 1) then
6140 if (t .lt. arci(j)) then
6141 arcsum = arcsum + arci(j) - t
6142 exang = exang + ex(ni)
6144 if (jb .ge. maxarc) then
6146 80 format (/,' SURFATOM -- Increase the Value',&
6151 kent(jb) = maxarc*i + k
6153 kout(jb) = maxarc*k + i
6162 arcsum = arcsum + pix2 - t
6164 exang = exang + ex(ni)
6167 kent(jb) = maxarc*i + k
6169 kout(jb) = maxarc*k + i
6176 arclen = arclen + gr(k)*arcsum
6179 if (arclen .eq. 0.0d0) goto 170
6180 if (jb .eq. 0) goto 150
6182 ! find number of independent boundaries and check connectivity
6186 if (kout(k) .ne. 0) then
6193 if (m .eq. kent(ii)) then
6196 if (j .eq. jb) goto 150
6208 ! attempt to fix connectivity error by moving atom slightly
6212 140 format (/,' SURFATOM -- Connectivity Error at Atom',i6)
6221 ! compute the exposed surface area for the sphere of interest
6224 area = ib*pix2 + exang + arclen
6225 area = mod(area,4.0d0*pi) * rrsq
6227 ! attempt to fix negative area by moving atom slightly
6229 if (area .lt. 0.0d0) then
6232 160 format (/,' SURFATOM -- Negative Area at Atom',i6)
6243 end subroutine surfatom
6244 !----------------------------------------------------------------
6245 !----------------------------------------------------------------
6246 subroutine alloc_MD_arrays
6247 !EL Allocation of arrays used by MD module
6249 integer :: nres2,nres6
6252 !----------------------
6256 allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
6257 allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
6258 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
6259 allocate(fricvec(nres2,nres2))
6260 allocate(pfric_mat(nres2,nres2),vfric_mat(nres2,nres2))
6261 allocate(afric_mat(nres2,nres2),prand_mat(nres2,nres2))
6262 allocate(vrand_mat1(nres2,nres2),vrand_mat2(nres2,nres2)) !(MAXRES2,MAXRES2)
6263 allocate(pfric0_mat(nres2,nres2,0:maxflag_stoch))
6264 allocate(afric0_mat(nres2,nres2,0:maxflag_stoch))
6265 allocate(vfric0_mat(nres2,nres2,0:maxflag_stoch))
6266 allocate(prand0_mat(nres2,nres2,0:maxflag_stoch))
6267 allocate(vrand0_mat1(nres2,nres2,0:maxflag_stoch))
6268 allocate(vrand0_mat2(nres2,nres2,0:maxflag_stoch)) !(MAXRES2,MAXRES2,0:maxflag_stoch)
6269 allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
6271 allocate(mt1(nres2,nres2),mt2(nres2,nres2),mt3(nres2,nres2)) !(maxres2,maxres2)
6272 !----------------------
6274 ! commom.langevin.lang0
6276 allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
6277 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
6278 allocate(fricvec(nres2,nres2)) !(MAXRES2,MAXRES2)
6279 allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
6280 allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
6283 if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2))
6284 !----------------------
6285 ! commom.hairpin in CSA module
6286 !----------------------
6287 ! common.mce in MCM_MD module
6288 !----------------------
6290 ! common /mdgrad/ in module.energy
6291 ! common /back_constr/ in module.energy
6292 ! common /qmeas/ in module.energy
6295 allocate(potEcomp(0:n_ene+4)) !(0:n_ene+4)
6297 allocate(d_t(3,0:nres2),d_a(3,0:nres2),d_t_old(3,0:nres2)) !(3,0:MAXRES2)
6298 allocate(d_a_work(nres6)) !(6*MAXRES)
6300 allocate(DM(nres2),DU1(nres2),DU2(nres2))
6301 allocate(DMorig(nres2),DU1orig(nres2),DU2orig(nres2))
6303 write (iout,*) "Before A Allocation",nfgtasks-1
6305 allocate(Gmat(nres2,nres2),A(nres2,nres2))
6306 if(.not.allocated(Ginv)) allocate(Ginv(nres2,nres2)) !in control: ergastulum
6307 allocate(Gsqrp(nres2,nres2),Gsqrm(nres2,nres2),Gvec(nres2,nres2)) !(maxres2,maxres2)
6309 allocate(Geigen(nres2)) !(maxres2)
6310 if(.not.allocated(vtot)) allocate(vtot(nres2)) !(maxres2)
6311 ! common /inertia/ in io_conf: parmread
6312 ! real(kind=8),dimension(:),allocatable :: ISC,msc !(ntyp+1)
6313 ! common /langevin/in io read_MDpar
6314 ! real(kind=8),dimension(:),allocatable :: gamsc !(ntyp1)
6315 ! real(kind=8),dimension(:),allocatable :: stdfsc !(ntyp)
6316 ! in io_conf: parmread
6317 ! real(kind=8),dimension(:),allocatable :: restok !(ntyp+1)
6318 ! common /mdpmpi/ in control: ergastulum
6319 if(.not.allocated(ng_start)) allocate(ng_start(0:nfgtasks-1))
6320 if(.not.allocated(ng_counts)) allocate(ng_counts(0:nfgtasks-1))
6321 if(.not.allocated(nginv_counts)) allocate(nginv_counts(0:nfgtasks-1)) !(0:MaxProcs-1)
6322 if(.not.allocated(nginv_start)) allocate(nginv_start(0:nfgtasks)) !(0:MaxProcs)
6323 !----------------------
6324 ! common.muca in read_muca
6325 ! common /double_muca/
6326 ! real(kind=8) :: elow,ehigh,factor,hbin,factor_min
6327 ! real(kind=8),dimension(:),allocatable :: emuca,nemuca,&
6328 ! nemuca2,hist !(4*maxres)
6329 ! common /integer_muca/
6330 ! integer :: nmuca,imtime,muca_smooth
6332 ! real(kind=8),dimension(:),allocatable :: elowi,ehighi !(maxprocs)
6333 !----------------------
6335 ! common /mdgrad/ in module.energy
6336 ! common /back_constr/ in module.energy
6337 ! common /qmeas/ in module.energy
6341 allocate(d_t_work(nres6),d_t_work_new(nres6),d_af_work(nres6))
6342 allocate(d_as_work(nres6),kinetic_force(nres6)) !(MAXRES6)
6343 allocate(d_t_new(3,0:nres2),d_a_old(3,0:nres2),d_a_short(3,0:nres2)) !,d_a !(3,0:MAXRES2)
6344 allocate(stdforcp(nres),stdforcsc(nres)) !(MAXRES)
6345 !----------------------
6347 allocate(D_ban(nres6)) !(MAXRES6) maxres6=6*maxres
6348 ! common /stochcalc/ stochforcvec
6349 allocate(stochforcvec(nres6)) !(MAXRES6) maxres6=6*maxres
6350 !----------------------
6352 end subroutine alloc_MD_arrays
6353 !-----------------------------------------------------------------------------
6354 !-----------------------------------------------------------------------------