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 !el 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(6*nres) :: Td !(maxres6) maxres6=6*maxres
191 real(kind=8),dimension(2*nres) :: ppvec !(maxres2) maxres2=2*maxres
192 !el common /stochcalc/ stochforcvec
193 !el real(kind=8),dimension(3) :: cm !el
194 !el common /gucio/ cm
196 logical :: lprn = .false.,lprn1 = .false.
197 integer :: maxiter = 5
198 real(kind=8) :: difftol = 1.0d-5
199 real(kind=8) :: xx,diffmax,blen2,diffbond,tt0
200 integer :: i,j,nbond,k,ind,ind1,iter
201 integer :: nres2,nres6
206 if (.not.allocated(stochforcvec)) allocate(stochforcvec(nres6)) !(MAXRES6) maxres6=6*maxres
210 if (itype(i).ne.10) nbond=nbond+1
214 write (iout,*) "Generalized inverse of fricmat"
215 call matout(dimen,dimen,nres6,nres6,fricmat)
227 Bmat(ind+j,ind1)=dC_norm(j,i)
232 if (itype(i).ne.10) then
235 Bmat(ind+j,ind1)=dC_norm(j,i+nres)
241 write (iout,*) "Matrix Bmat"
242 call MATOUT(nbond,dimen,nres6,nres6,Bmat)
248 GBmat(i,j)=GBmat(i,j)+fricmat(i,k)*Bmat(k,j)
253 write (iout,*) "Matrix GBmat"
254 call MATOUT(nbond,dimen,nres6,nres2,Gbmat)
260 Cmat_(i,j)=Cmat_(i,j)+Bmat(k,i)*GBmat(k,j)
265 write (iout,*) "Matrix Cmat"
266 call MATOUT(nbond,nbond,nres2,nres2,Cmat_)
268 call matinvert(nbond,nres2,Cmat_,Cinv,osob)
270 write (iout,*) "Matrix Cinv"
271 call MATOUT(nbond,nbond,nres2,nres2,Cinv)
277 Tmat(i,j)=Tmat(i,j)+GBmat(i,k)*Cinv(k,j)
282 write (iout,*) "Matrix Tmat"
283 call MATOUT(nbond,dimen,nres6,nres2,Tmat)
293 Pmat(i,j)=Pmat(i,j)-Tmat(i,k)*Bmat(j,k)
298 write (iout,*) "Matrix Pmat"
299 call MATOUT(dimen,dimen,nres6,nres6,Pmat)
306 Td(i)=Td(i)+vbl*Tmat(i,ind)
309 if (itype(k).ne.10) then
311 Td(i)=Td(i)+vbldsc0(1,itype(k))*Tmat(i,ind)
316 write (iout,*) "Vector Td"
318 write (iout,'(i5,f10.5)') i,Td(i)
321 call stochastic_force(stochforcvec)
323 write (iout,*) "stochforcvec"
325 write (iout,*) i,stochforcvec(i)
329 zapas(j)=-gcart(j,0)+stochforcvec(j)
331 dC_work(j)=dC_old(j,0)
337 zapas(ind)=-gcart(j,i)+stochforcvec(ind)
338 dC_work(ind)=dC_old(j,i)
342 if (itype(i).ne.10) then
345 zapas(ind)=-gxcart(j,i)+stochforcvec(ind)
346 dC_work(ind)=dC_old(j,i+nres)
352 write (iout,*) "Initial d_t_work"
354 write (iout,*) i,d_t_work(i)
361 d_t_work(i)=d_t_work(i)+fricmat(i,j)*zapas(j)
368 zapas(i)=zapas(i)+Pmat(i,j)*(dC_work(j)+d_t_work(j)*d_time)
372 write (iout,*) "Final d_t_work and zapas"
374 write (iout,*) i,d_t_work(i),zapas(i)
388 dc_work(ind+j)=dc(j,i)
394 d_t(j,i+nres)=d_t_work(ind+j)
395 dc(j,i+nres)=zapas(ind+j)
396 dc_work(ind+j)=dc(j,i+nres)
402 write (iout,*) "Before correction for rotational lengthening"
403 write (iout,*) "New coordinates",&
404 " and differences between actual and standard bond lengths"
409 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
413 if (itype(i).ne.10) then
415 xx=vbld(i+nres)-vbldsc0(1,itype(i))
416 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
417 i,(dC(j,i+nres),j=1,3),xx
421 ! Second correction (rotational lengthening)
427 blen2 = scalar(dc(1,i),dc(1,i))
428 ppvec(ind)=2*vbl**2-blen2
429 diffbond=dabs(vbl-dsqrt(blen2))
430 if (diffbond.gt.diffmax) diffmax=diffbond
431 if (ppvec(ind).gt.0.0d0) then
432 ppvec(ind)=dsqrt(ppvec(ind))
437 write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind)
441 if (itype(i).ne.10) then
443 blen2 = scalar(dc(1,i+nres),dc(1,i+nres))
444 ppvec(ind)=2*vbldsc0(1,itype(i))**2-blen2
445 diffbond=dabs(vbldsc0(1,itype(i))-dsqrt(blen2))
446 if (diffbond.gt.diffmax) diffmax=diffbond
447 if (ppvec(ind).gt.0.0d0) then
448 ppvec(ind)=dsqrt(ppvec(ind))
453 write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind)
457 if (lprn) write (iout,*) "iter",iter," diffmax",diffmax
458 if (diffmax.lt.difftol) goto 10
462 Td(i)=Td(i)+ppvec(j)*Tmat(i,j)
468 zapas(i)=zapas(i)+Pmat(i,j)*dc_work(j)
479 dc_work(ind+j)=zapas(ind+j)
484 if (itype(i).ne.10) then
486 dc(j,i+nres)=zapas(ind+j)
487 dc_work(ind+j)=zapas(ind+j)
492 ! Building the chain from the newly calculated coordinates
495 if (large.and. mod(itime,ntwe).eq.0) then
496 write (iout,*) "Cartesian and internal coordinates: step 1"
499 write (iout,'(a)') "Potential forces"
501 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(-gcart(j,i),j=1,3),&
504 write (iout,'(a)') "Stochastic forces"
506 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(stochforc(j,i),j=1,3),&
507 (stochforc(j,i+nres),j=1,3)
509 write (iout,'(a)') "Velocities"
511 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
512 (d_t(j,i+nres),j=1,3)
517 write (iout,*) "After correction for rotational lengthening"
518 write (iout,*) "New coordinates",&
519 " and differences between actual and standard bond lengths"
524 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
528 if (itype(i).ne.10) then
530 xx=vbld(i+nres)-vbldsc0(1,itype(i))
531 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
532 i,(dC(j,i+nres),j=1,3),xx
537 ! write (iout,*) "Too many attempts at correcting the bonds"
545 ! Calculate energy and forces
547 call etotal(potEcomp)
548 potE=potEcomp(0)-potEcomp(20)
551 ! Calculate the kinetic and total energy and the kinetic temperature
554 t_enegrad=t_enegrad+MPI_Wtime()-tt0
556 t_enegrad=t_enegrad+tcpu()-tt0
559 kinetic_T=2.0d0/(dimen*Rb)*EK
561 end subroutine brown_step
562 !-----------------------------------------------------------------------------
564 !-----------------------------------------------------------------------------
565 subroutine gauss(RO,AP,MT,M,N,*)
567 ! CALCULATES (RO**(-1))*AP BY GAUSS ELIMINATION
568 ! RO IS A SQUARE MATRIX
569 ! THE CALCULATED PRODUCT IS STORED IN AP
570 ! ABNORMAL EXIT IF RO IS SINGULAR
572 integer :: MT, M, N, M1,I,J,IM,&
574 real(kind=8) :: RO(MT,M),AP(MT,N),X,RM,PR,Y
580 if(dabs(X).le.1.0D-13) return 1
592 if(DABS(RO(J,I)).LE.RM) goto 2
606 if(dabs(X).le.1.0E-13) return 1
615 8 AP(J,K)=AP(J,K)-Y*AP(I,K)
617 9 RO(J,K)=RO(J,K)-Y*RO(I,K)
621 if(dabs(X).le.1.0E-13) return 1
631 15 X=X-AP(K,J)*RO(MI,K)
636 !-----------------------------------------------------------------------------
638 !-----------------------------------------------------------------------------
639 subroutine kinetic(KE_total)
640 !----------------------------------------------------------------
641 ! This subroutine calculates the total kinetic energy of the chain
642 !-----------------------------------------------------------------
644 ! implicit real*8 (a-h,o-z)
645 ! include 'DIMENSIONS'
646 ! include 'COMMON.VAR'
647 ! include 'COMMON.CHAIN'
648 ! include 'COMMON.DERIV'
649 ! include 'COMMON.GEO'
650 ! include 'COMMON.LOCAL'
651 ! include 'COMMON.INTERACT'
652 ! include 'COMMON.MD'
653 ! include 'COMMON.IOUNITS'
654 real(kind=8) :: KE_total
657 real(kind=8) :: KEt_p,KEt_sc,KEr_p,KEr_sc,incr(3),&
662 ! write (iout,*) "ISC",(isc(itype(i)),i=1,nres)
663 ! The translational part for peptide virtual bonds
668 ! write (iout,*) "Kinetic trp:",i,(incr(j),j=1,3)
670 v(j)=incr(j)+0.5d0*d_t(j,i)
672 vtot(i)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
673 KEt_p=KEt_p+(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
675 incr(j)=incr(j)+d_t(j,i)
678 ! write(iout,*) 'KEt_p', KEt_p
679 ! The translational part for the side chain virtual bond
680 ! Only now we can initialize incr with zeros. It must be equal
681 ! to the velocities of the first Calpha.
687 if (itype(i).eq.10) then
693 v(j)=incr(j)+d_t(j,nres+i)
696 ! write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3)
697 ! write (iout,*) "i",i," msc",msc(iti)," v",(v(j),j=1,3)
698 KEt_sc=KEt_sc+msc(iti)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
699 vtot(i+nres)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
701 incr(j)=incr(j)+d_t(j,i)
705 ! write(iout,*) 'KEt_sc', KEt_sc
706 ! The part due to stretching and rotation of the peptide groups
709 ! write (iout,*) "i",i
710 ! write (iout,*) "i",i," mag1",mag1," mag2",mag2
714 ! write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3)
715 KEr_p=KEr_p+(incr(1)*incr(1)+incr(2)*incr(2) &
719 ! write(iout,*) 'KEr_p', KEr_p
720 ! The rotational part of the side chain virtual bond
724 if (itype(i).ne.10) then
726 incr(j)=d_t(j,nres+i)
728 ! write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3)
729 KEr_sc=KEr_sc+Isc(iti)*(incr(1)*incr(1)+incr(2)*incr(2)+ &
733 ! The total kinetic energy
735 ! write(iout,*) 'KEr_sc', KEr_sc
736 KE_total=0.5d0*(mp*KEt_p+KEt_sc+0.25d0*Ip*KEr_p+KEr_sc)
737 ! write (iout,*) "KE_total",KE_total
739 end subroutine kinetic
740 !-----------------------------------------------------------------------------
742 !-----------------------------------------------------------------------------
744 !------------------------------------------------
745 ! The driver for molecular dynamics subroutines
746 !------------------------------------------------
749 use control, only:tcpu,ovrtim
750 ! use io_comm, only:ilen
752 use compare, only:secondary2,hairpin
753 use io, only:cartout,statout
754 ! implicit real*8 (a-h,o-z)
755 ! include 'DIMENSIONS'
758 integer :: IERROR,ERRCODE
760 ! include 'COMMON.SETUP'
761 ! include 'COMMON.CONTROL'
762 ! include 'COMMON.VAR'
763 ! include 'COMMON.MD'
765 ! include 'COMMON.LANGEVIN'
767 ! include 'COMMON.LANGEVIN.lang0'
769 ! include 'COMMON.CHAIN'
770 ! include 'COMMON.DERIV'
771 ! include 'COMMON.GEO'
772 ! include 'COMMON.LOCAL'
773 ! include 'COMMON.INTERACT'
774 ! include 'COMMON.IOUNITS'
775 ! include 'COMMON.NAMES'
776 ! include 'COMMON.TIME1'
777 ! include 'COMMON.HAIRPIN'
778 real(kind=8),dimension(3) :: L,vcm
780 real(kind=8),dimension(6*nres) :: v_work,v_transf !(maxres6) maxres6=6*maxres
782 integer :: rstcount !ilen,
784 character(len=50) :: tytul
785 !el common /gucio/ cm
786 integer :: itime,i,j,nharp
787 integer,dimension(4,nres/3) :: iharp !(4,nres/3)(4,maxres/3)
789 real(kind=8) :: tt0,scalfac
794 if (ilen(tmpdir).gt.0) &
795 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_" &
796 //liczba(:ilen(liczba))//'.rst')
798 if (ilen(tmpdir).gt.0) &
799 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"//'.rst')
806 write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
812 ! Determine the inverse of the inertia matrix.
813 call setup_MD_matrices
817 t_MDsetup = MPI_Wtime()-tt0
819 t_MDsetup = tcpu()-tt0
822 ! Entering the MD loop
828 if (lang.eq.2 .or. lang.eq.3) then
832 call sd_verlet_p_setup
834 call sd_verlet_ciccotti_setup
838 pfric0_mat(i,j,0)=pfric_mat(i,j)
839 afric0_mat(i,j,0)=afric_mat(i,j)
840 vfric0_mat(i,j,0)=vfric_mat(i,j)
841 prand0_mat(i,j,0)=prand_mat(i,j)
842 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
843 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
848 flag_stoch(i)=.false.
852 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
854 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
858 else if (lang.eq.1 .or. lang.eq.4) then
862 t_langsetup=MPI_Wtime()-tt0
865 t_langsetup=tcpu()-tt0
868 do itime=1,n_timestep
870 if (large.and. mod(itime,ntwe).eq.0) &
871 write (iout,*) "itime",itime
873 if (lang.gt.0 .and. surfarea .and. &
874 mod(itime,reset_fricmat).eq.0) then
875 if (lang.eq.2 .or. lang.eq.3) then
879 call sd_verlet_p_setup
881 call sd_verlet_ciccotti_setup
885 pfric0_mat(i,j,0)=pfric_mat(i,j)
886 afric0_mat(i,j,0)=afric_mat(i,j)
887 vfric0_mat(i,j,0)=vfric_mat(i,j)
888 prand0_mat(i,j,0)=prand_mat(i,j)
889 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
890 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
895 flag_stoch(i)=.false.
898 else if (lang.eq.1 .or. lang.eq.4) then
901 write (iout,'(a,i10)') &
902 "Friction matrix reset based on surface area, itime",itime
904 if (reset_vel .and. tbf .and. lang.eq.0 &
905 .and. mod(itime,count_reset_vel).eq.0) then
907 write(iout,'(a,f20.2)') &
908 "Velocities reset to random values, time",totT
911 d_t_old(j,i)=d_t(j,i)
915 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
919 d_t(j,0)=d_t(j,0)-vcm(j)
922 kinetic_T=2.0d0/(dimen3*Rb)*EK
923 scalfac=dsqrt(T_bath/kinetic_T)
924 write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
927 d_t_old(j,i)=scalfac*d_t(j,i)
933 ! Time-reversible RESPA algorithm
934 ! (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
935 call RESPA_step(itime)
937 ! Variable time step algorithm.
938 call velverlet_step(itime)
942 call brown_step(itime)
944 print *,"Brown dynamics not here!"
946 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
952 if (mod(itime,ntwe).eq.0) call statout(itime)
965 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
968 v_work(ind)=d_t(j,i+nres)
973 write (66,'(80f10.5)') &
974 ((d_t(j,i),j=1,3),i=0,nres-1),((d_t(j,i+nres),j=1,3),i=1,nres)
978 v_transf(i)=v_transf(i)+gvec(j,i)*v_work(j)
980 v_transf(i)= v_transf(i)*dsqrt(geigen(i))
982 write (67,'(80f10.5)') (v_transf(i),i=1,ind)
985 if (mod(itime,ntwx).eq.0) then
986 write (tytul,'("time",f8.2)') totT
988 call hairpin(.true.,nharp,iharp)
989 call secondary2(.true.)
990 call pdbout(potE,tytul,ipdb)
995 if (rstcount.eq.1000.or.itime.eq.n_timestep) then
996 open(irest2,file=rest2name,status='unknown')
997 write(irest2,*) totT,EK,potE,totE,t_bath
999 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
1002 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
1010 t_MD=MPI_Wtime()-tt0
1014 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))') &
1016 'MD calculations setup:',t_MDsetup,&
1017 'Energy & gradient evaluation:',t_enegrad,&
1018 'Stochastic MD setup:',t_langsetup,&
1019 'Stochastic MD step setup:',t_sdsetup,&
1021 write (iout,'(/28(1h=),a25,27(1h=))') &
1022 ' End of MD calculation '
1024 write (iout,*) "time for etotal",t_etotal," elong",t_elong,&
1026 write (iout,*) "time_fric",time_fric," time_stoch",time_stoch,&
1027 " time_fricmatmult",time_fricmatmult," time_fsample ",&
1032 !-----------------------------------------------------------------------------
1033 subroutine velverlet_step(itime)
1034 !-------------------------------------------------------------------------------
1035 ! Perform a single velocity Verlet step; the time step can be rescaled if
1036 ! increments in accelerations exceed the threshold
1037 !-------------------------------------------------------------------------------
1038 ! implicit real*8 (a-h,o-z)
1039 ! include 'DIMENSIONS'
1041 use control, only:tcpu
1045 integer :: ierror,ierrcode
1046 real(kind=8) :: errcode
1048 ! include 'COMMON.SETUP'
1049 ! include 'COMMON.VAR'
1050 ! include 'COMMON.MD'
1052 ! include 'COMMON.LANGEVIN'
1054 ! include 'COMMON.LANGEVIN.lang0'
1056 ! include 'COMMON.CHAIN'
1057 ! include 'COMMON.DERIV'
1058 ! include 'COMMON.GEO'
1059 ! include 'COMMON.LOCAL'
1060 ! include 'COMMON.INTERACT'
1061 ! include 'COMMON.IOUNITS'
1062 ! include 'COMMON.NAMES'
1063 ! include 'COMMON.TIME1'
1064 ! include 'COMMON.MUCA'
1065 real(kind=8),dimension(3) :: vcm,incr
1066 real(kind=8),dimension(3) :: L
1067 integer :: count,rstcount !ilen,
1069 character(len=50) :: tytul
1070 integer :: maxcount_scale = 20
1071 !el common /gucio/ cm
1072 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1073 !el common /stochcalc/ stochforcvec
1074 integer :: itime,icount_scale,itime_scal,i,j,ifac_time
1076 real(kind=8) :: epdrift,tt0,fac_time
1078 if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres)) !(MAXRES6) maxres6=6*maxres
1084 else if (lang.eq.2 .or. lang.eq.3) then
1086 call stochastic_force(stochforcvec)
1089 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
1091 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1098 icount_scale=icount_scale+1
1099 if (icount_scale.gt.maxcount_scale) then
1101 "ERROR: too many attempts at scaling down the time step. ",&
1102 "amax=",amax,"epdrift=",epdrift,&
1103 "damax=",damax,"edriftmax=",edriftmax,&
1107 call MPI_Abort(MPI_COMM_WORLD,IERROR,IERRCODE)
1111 ! First step of the velocity Verlet algorithm
1116 else if (lang.eq.3) then
1118 call sd_verlet1_ciccotti
1120 else if (lang.eq.1) then
1125 ! Build the chain from the newly calculated coordinates
1126 call chainbuild_cart
1127 if (rattle) call rattle1
1129 if (large.and. mod(itime,ntwe).eq.0) then
1130 write (iout,*) "Cartesian and internal coordinates: step 1"
1135 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),&
1136 (dc(j,i+nres),j=1,3)
1138 write (iout,*) "Accelerations"
1140 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1141 (d_a(j,i+nres),j=1,3)
1143 write (iout,*) "Velocities, step 1"
1145 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1146 (d_t(j,i+nres),j=1,3)
1155 ! Calculate energy and forces
1157 call etotal(potEcomp)
1158 if (large.and. mod(itime,ntwe).eq.0) &
1159 call enerprint(potEcomp)
1162 t_etotal=t_etotal+MPI_Wtime()-tt0
1164 t_etotal=t_etotal+tcpu()-tt0
1167 potE=potEcomp(0)-potEcomp(20)
1169 ! Get the new accelerations
1172 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1174 t_enegrad=t_enegrad+tcpu()-tt0
1176 ! Determine maximum acceleration and scale down the timestep if needed
1178 amax=amax/(itime_scal+1)**2
1179 call predict_edrift(epdrift)
1180 if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then
1181 ! Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
1183 ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax)) &
1185 itime_scal=itime_scal+ifac_time
1186 ! fac_time=dmin1(damax/amax,0.5d0)
1187 fac_time=0.5d0**ifac_time
1188 d_time=d_time*fac_time
1189 if (lang.eq.2 .or. lang.eq.3) then
1191 ! write (iout,*) "Calling sd_verlet_setup: 1"
1192 ! Rescale the stochastic forces and recalculate or restore
1193 ! the matrices of tinker integrator
1194 if (itime_scal.gt.maxflag_stoch) then
1195 if (large) write (iout,'(a,i5,a)') &
1196 "Calculate matrices for stochastic step;",&
1197 " itime_scal ",itime_scal
1199 call sd_verlet_p_setup
1201 call sd_verlet_ciccotti_setup
1203 write (iout,'(2a,i3,a,i3,1h.)') &
1204 "Warning: cannot store matrices for stochastic",&
1205 " integration because the index",itime_scal,&
1206 " is greater than",maxflag_stoch
1207 write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",&
1208 " integration Langevin algorithm for better efficiency."
1209 else if (flag_stoch(itime_scal)) then
1210 if (large) write (iout,'(a,i5,a,l1)') &
1211 "Restore matrices for stochastic step; itime_scal ",&
1212 itime_scal," flag ",flag_stoch(itime_scal)
1215 pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
1216 afric_mat(i,j)=afric0_mat(i,j,itime_scal)
1217 vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
1218 prand_mat(i,j)=prand0_mat(i,j,itime_scal)
1219 vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
1220 vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
1224 if (large) write (iout,'(2a,i5,a,l1)') &
1225 "Calculate & store matrices for stochastic step;",&
1226 " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
1228 call sd_verlet_p_setup
1230 call sd_verlet_ciccotti_setup
1232 flag_stoch(ifac_time)=.true.
1235 pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
1236 afric0_mat(i,j,itime_scal)=afric_mat(i,j)
1237 vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
1238 prand0_mat(i,j,itime_scal)=prand_mat(i,j)
1239 vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
1240 vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
1244 fac_time=1.0d0/dsqrt(fac_time)
1246 stochforcvec(i)=fac_time*stochforcvec(i)
1249 else if (lang.eq.1) then
1250 ! Rescale the accelerations due to stochastic forces
1251 fac_time=1.0d0/dsqrt(fac_time)
1253 d_as_work(i)=d_as_work(i)*fac_time
1256 if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)') &
1257 "itime",itime," Timestep scaled down to ",&
1258 d_time," ifac_time",ifac_time," itime_scal",itime_scal
1260 ! Second step of the velocity Verlet algorithm
1265 else if (lang.eq.3) then
1267 call sd_verlet2_ciccotti
1269 else if (lang.eq.1) then
1274 if (rattle) call rattle2
1276 if (d_time.ne.d_time0) then
1279 if (lang.eq.2 .or. lang.eq.3) then
1280 if (large) write (iout,'(a)') &
1281 "Restore original matrices for stochastic step"
1282 ! write (iout,*) "Calling sd_verlet_setup: 2"
1283 ! Restore the matrices of tinker integrator if the time step has been restored
1286 pfric_mat(i,j)=pfric0_mat(i,j,0)
1287 afric_mat(i,j)=afric0_mat(i,j,0)
1288 vfric_mat(i,j)=vfric0_mat(i,j,0)
1289 prand_mat(i,j)=prand0_mat(i,j,0)
1290 vrand_mat1(i,j)=vrand0_mat1(i,j,0)
1291 vrand_mat2(i,j)=vrand0_mat2(i,j,0)
1300 ! Calculate the kinetic and the total energy and the kinetic temperature
1304 ! call kinetic1(EK1)
1305 ! write (iout,*) "step",itime," EK",EK," EK1",EK1
1307 ! Couple the system to Berendsen bath if needed
1308 if (tbf .and. lang.eq.0) then
1311 kinetic_T=2.0d0/(dimen3*Rb)*EK
1312 ! Backup the coordinates, velocities, and accelerations
1316 d_t_old(j,i)=d_t(j,i)
1317 d_a_old(j,i)=d_a(j,i)
1321 if (mod(itime,ntwe).eq.0 .and. large) then
1322 write (iout,*) "Velocities, step 2"
1324 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1325 (d_t(j,i+nres),j=1,3)
1330 end subroutine velverlet_step
1331 !-----------------------------------------------------------------------------
1332 subroutine RESPA_step(itime)
1333 !-------------------------------------------------------------------------------
1334 ! Perform a single RESPA step.
1335 !-------------------------------------------------------------------------------
1336 ! implicit real*8 (a-h,o-z)
1337 ! include 'DIMENSIONS'
1341 use control, only:tcpu
1343 ! use io_conf, only:cartprint
1346 integer :: IERROR,ERRCODE
1348 ! include 'COMMON.SETUP'
1349 ! include 'COMMON.CONTROL'
1350 ! include 'COMMON.VAR'
1351 ! include 'COMMON.MD'
1353 ! include 'COMMON.LANGEVIN'
1355 ! include 'COMMON.LANGEVIN.lang0'
1357 ! include 'COMMON.CHAIN'
1358 ! include 'COMMON.DERIV'
1359 ! include 'COMMON.GEO'
1360 ! include 'COMMON.LOCAL'
1361 ! include 'COMMON.INTERACT'
1362 ! include 'COMMON.IOUNITS'
1363 ! include 'COMMON.NAMES'
1364 ! include 'COMMON.TIME1'
1365 real(kind=8),dimension(0:n_ene) :: energia_short,energia_long
1366 real(kind=8),dimension(3) :: L,vcm,incr
1367 real(kind=8),dimension(3,0:2*nres) :: dc_old0,d_t_old0,d_a_old0 !(3,0:maxres2) maxres2=2*maxres
1368 logical :: PRINT_AMTS_MSG = .false.
1369 integer :: count,rstcount !ilen,
1371 character(len=50) :: tytul
1372 integer :: maxcount_scale = 10
1373 !el common /gucio/ cm
1374 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1375 !el common /stochcalc/ stochforcvec
1376 integer :: itime,itt,i,j,itsplit
1378 !el common /cipiszcze/ itt
1380 real(kind=8) :: epdrift,tt0,epdriftmax
1383 if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres)) !(MAXRES6) maxres6=6*maxres
1387 if (large.and. mod(itime,ntwe).eq.0) then
1388 write (iout,*) "***************** RESPA itime",itime
1389 write (iout,*) "Cartesian and internal coordinates: step 0"
1391 call pdbout(0.0d0,"cipiszcze",iout)
1393 write (iout,*) "Accelerations from long-range forces"
1395 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1396 (d_a(j,i+nres),j=1,3)
1398 write (iout,*) "Velocities, step 0"
1400 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1401 (d_t(j,i+nres),j=1,3)
1406 ! Perform the initial RESPA step (increment velocities)
1407 ! write (iout,*) "*********************** RESPA ini"
1410 if (mod(itime,ntwe).eq.0 .and. large) then
1411 write (iout,*) "Velocities, end"
1413 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1414 (d_t(j,i+nres),j=1,3)
1418 ! Compute the short-range forces
1424 ! 7/2/2009 commented out
1426 ! call etotal_short(energia_short)
1429 ! 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
1432 d_a(j,i)=d_a_short(j,i)
1436 if (large.and. mod(itime,ntwe).eq.0) then
1437 write (iout,*) "energia_short",energia_short(0)
1438 write (iout,*) "Accelerations from short-range forces"
1440 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1441 (d_a(j,i+nres),j=1,3)
1446 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1448 t_enegrad=t_enegrad+tcpu()-tt0
1453 d_t_old(j,i)=d_t(j,i)
1454 d_a_old(j,i)=d_a(j,i)
1457 ! 6/30/08 A-MTS: attempt at increasing the split number
1460 dc_old0(j,i)=dc_old(j,i)
1461 d_t_old0(j,i)=d_t_old(j,i)
1462 d_a_old0(j,i)=d_a_old(j,i)
1465 if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
1466 if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
1473 ! write (iout,*) "itime",itime," ntime_split",ntime_split
1474 ! Split the time step
1475 d_time=d_time0/ntime_split
1476 ! Perform the short-range RESPA steps (velocity Verlet increments of
1477 ! positions and velocities using short-range forces)
1478 ! write (iout,*) "*********************** RESPA split"
1479 do itsplit=1,ntime_split
1482 else if (lang.eq.2 .or. lang.eq.3) then
1484 call stochastic_force(stochforcvec)
1487 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
1489 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1494 ! First step of the velocity Verlet algorithm
1499 else if (lang.eq.3) then
1501 call sd_verlet1_ciccotti
1503 else if (lang.eq.1) then
1508 ! Build the chain from the newly calculated coordinates
1509 call chainbuild_cart
1510 if (rattle) call rattle1
1512 if (large.and. mod(itime,ntwe).eq.0) then
1513 write (iout,*) "***** ITSPLIT",itsplit
1514 write (iout,*) "Cartesian and internal coordinates: step 1"
1515 call pdbout(0.0d0,"cipiszcze",iout)
1518 write (iout,*) "Velocities, step 1"
1520 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1521 (d_t(j,i+nres),j=1,3)
1530 ! Calculate energy and forces
1532 call etotal_short(energia_short)
1533 if (large.and. mod(itime,ntwe).eq.0) &
1534 call enerprint(energia_short)
1537 t_eshort=t_eshort+MPI_Wtime()-tt0
1539 t_eshort=t_eshort+tcpu()-tt0
1543 ! Get the new accelerations
1545 ! 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
1548 d_a_short(j,i)=d_a(j,i)
1552 if (large.and. mod(itime,ntwe).eq.0) then
1553 write (iout,*)"energia_short",energia_short(0)
1554 write (iout,*) "Accelerations from short-range forces"
1556 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1557 (d_a(j,i+nres),j=1,3)
1562 ! Determine maximum acceleration and scale down the timestep if needed
1564 amax=amax/ntime_split**2
1565 call predict_edrift(epdrift)
1566 if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0) &
1567 write (iout,*) "amax",amax," damax",damax,&
1568 " epdrift",epdrift," epdriftmax",epdriftmax
1569 ! Exit loop and try with increased split number if the change of
1570 ! acceleration is too big
1571 if (amax.gt.damax .or. epdrift.gt.edriftmax) then
1572 if (ntime_split.lt.maxtime_split) then
1574 ntime_split=ntime_split*2
1577 dc_old(j,i)=dc_old0(j,i)
1578 d_t_old(j,i)=d_t_old0(j,i)
1579 d_a_old(j,i)=d_a_old0(j,i)
1582 if (PRINT_AMTS_MSG) then
1583 write (iout,*) "acceleration/energy drift too large",amax,&
1584 epdrift," split increased to ",ntime_split," itime",itime,&
1590 "Uh-hu. Bumpy landscape. Maximum splitting number",&
1592 " already reached!!! Trying to carry on!"
1596 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1598 t_enegrad=t_enegrad+tcpu()-tt0
1600 ! Second step of the velocity Verlet algorithm
1605 else if (lang.eq.3) then
1607 call sd_verlet2_ciccotti
1609 else if (lang.eq.1) then
1614 if (rattle) call rattle2
1615 ! Backup the coordinates, velocities, and accelerations
1619 d_t_old(j,i)=d_t(j,i)
1620 d_a_old(j,i)=d_a(j,i)
1627 ! Restore the time step
1629 ! Compute long-range forces
1636 call etotal_long(energia_long)
1637 if (large.and. mod(itime,ntwe).eq.0) &
1638 call enerprint(energia_long)
1641 t_elong=t_elong+MPI_Wtime()-tt0
1643 t_elong=t_elong+tcpu()-tt0
1649 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1651 t_enegrad=t_enegrad+tcpu()-tt0
1653 ! Compute accelerations from long-range forces
1655 if (large.and. mod(itime,ntwe).eq.0) then
1656 write (iout,*) "energia_long",energia_long(0)
1657 write (iout,*) "Cartesian and internal coordinates: step 2"
1659 call pdbout(0.0d0,"cipiszcze",iout)
1661 write (iout,*) "Accelerations from long-range forces"
1663 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1664 (d_a(j,i+nres),j=1,3)
1666 write (iout,*) "Velocities, step 2"
1668 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1669 (d_t(j,i+nres),j=1,3)
1673 ! Compute the final RESPA step (increment velocities)
1674 ! write (iout,*) "*********************** RESPA fin"
1676 ! Compute the complete potential energy
1678 potEcomp(i)=energia_short(i)+energia_long(i)
1680 potE=potEcomp(0)-potEcomp(20)
1681 ! potE=energia_short(0)+energia_long(0)
1683 ! Calculate the kinetic and the total energy and the kinetic temperature
1686 ! Couple the system to Berendsen bath if needed
1687 if (tbf .and. lang.eq.0) then
1690 kinetic_T=2.0d0/(dimen3*Rb)*EK
1691 ! Backup the coordinates, velocities, and accelerations
1693 if (mod(itime,ntwe).eq.0 .and. large) then
1694 write (iout,*) "Velocities, end"
1696 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1697 (d_t(j,i+nres),j=1,3)
1702 end subroutine RESPA_step
1703 !-----------------------------------------------------------------------------
1704 subroutine RESPA_vel
1705 ! First and last RESPA step (incrementing velocities using long-range
1708 ! implicit real*8 (a-h,o-z)
1709 ! include 'DIMENSIONS'
1710 ! include 'COMMON.CONTROL'
1711 ! include 'COMMON.VAR'
1712 ! include 'COMMON.MD'
1713 ! include 'COMMON.CHAIN'
1714 ! include 'COMMON.DERIV'
1715 ! include 'COMMON.GEO'
1716 ! include 'COMMON.LOCAL'
1717 ! include 'COMMON.INTERACT'
1718 ! include 'COMMON.IOUNITS'
1719 ! include 'COMMON.NAMES'
1720 integer :: i,j,inres
1723 d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
1727 d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
1731 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1734 d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
1739 end subroutine RESPA_vel
1740 !-----------------------------------------------------------------------------
1742 ! Applying velocity Verlet algorithm - step 1 to coordinates
1744 ! implicit real*8 (a-h,o-z)
1745 ! include 'DIMENSIONS'
1746 ! include 'COMMON.CONTROL'
1747 ! include 'COMMON.VAR'
1748 ! include 'COMMON.MD'
1749 ! include 'COMMON.CHAIN'
1750 ! include 'COMMON.DERIV'
1751 ! include 'COMMON.GEO'
1752 ! include 'COMMON.LOCAL'
1753 ! include 'COMMON.INTERACT'
1754 ! include 'COMMON.IOUNITS'
1755 ! include 'COMMON.NAMES'
1756 real(kind=8) :: adt,adt2
1757 integer :: i,j,inres
1760 write (iout,*) "VELVERLET1 START: DC"
1762 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),&
1763 (dc(j,i+nres),j=1,3)
1767 adt=d_a_old(j,0)*d_time
1769 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1770 d_t_new(j,0)=d_t_old(j,0)+adt2
1771 d_t(j,0)=d_t_old(j,0)+adt
1775 adt=d_a_old(j,i)*d_time
1777 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1778 d_t_new(j,i)=d_t_old(j,i)+adt2
1779 d_t(j,i)=d_t_old(j,i)+adt
1783 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1786 adt=d_a_old(j,inres)*d_time
1788 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1789 d_t_new(j,inres)=d_t_old(j,inres)+adt2
1790 d_t(j,inres)=d_t_old(j,inres)+adt
1795 write (iout,*) "VELVERLET1 END: DC"
1797 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),&
1798 (dc(j,i+nres),j=1,3)
1802 end subroutine verlet1
1803 !-----------------------------------------------------------------------------
1805 ! Step 2 of the velocity Verlet algorithm: update velocities
1807 ! implicit real*8 (a-h,o-z)
1808 ! include 'DIMENSIONS'
1809 ! include 'COMMON.CONTROL'
1810 ! include 'COMMON.VAR'
1811 ! include 'COMMON.MD'
1812 ! include 'COMMON.CHAIN'
1813 ! include 'COMMON.DERIV'
1814 ! include 'COMMON.GEO'
1815 ! include 'COMMON.LOCAL'
1816 ! include 'COMMON.INTERACT'
1817 ! include 'COMMON.IOUNITS'
1818 ! include 'COMMON.NAMES'
1819 integer :: i,j,inres
1822 d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
1826 d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
1830 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1833 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
1838 end subroutine verlet2
1839 !-----------------------------------------------------------------------------
1840 subroutine sddir_precalc
1841 ! Applying velocity Verlet algorithm - step 1 to coordinates
1842 ! implicit real*8 (a-h,o-z)
1843 ! include 'DIMENSIONS'
1849 ! include 'COMMON.CONTROL'
1850 ! include 'COMMON.VAR'
1851 ! include 'COMMON.MD'
1853 ! include 'COMMON.LANGEVIN'
1855 ! include 'COMMON.LANGEVIN.lang0'
1857 ! include 'COMMON.CHAIN'
1858 ! include 'COMMON.DERIV'
1859 ! include 'COMMON.GEO'
1860 ! include 'COMMON.LOCAL'
1861 ! include 'COMMON.INTERACT'
1862 ! include 'COMMON.IOUNITS'
1863 ! include 'COMMON.NAMES'
1864 ! include 'COMMON.TIME1'
1865 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1866 !el common /stochcalc/ stochforcvec
1867 real(kind=8) :: time00
1869 ! Compute friction and stochastic forces
1874 time_fric=time_fric+MPI_Wtime()-time00
1876 call stochastic_force(stochforcvec)
1877 time_stoch=time_stoch+MPI_Wtime()-time00
1880 ! Compute the acceleration due to friction forces (d_af_work) and stochastic
1881 ! forces (d_as_work)
1883 call ginv_mult(fric_work, d_af_work)
1884 call ginv_mult(stochforcvec, d_as_work)
1886 end subroutine sddir_precalc
1887 !-----------------------------------------------------------------------------
1888 subroutine sddir_verlet1
1889 ! Applying velocity Verlet algorithm - step 1 to velocities
1892 ! implicit real*8 (a-h,o-z)
1893 ! include 'DIMENSIONS'
1894 ! include 'COMMON.CONTROL'
1895 ! include 'COMMON.VAR'
1896 ! include 'COMMON.MD'
1898 ! include 'COMMON.LANGEVIN'
1900 ! include 'COMMON.LANGEVIN.lang0'
1902 ! include 'COMMON.CHAIN'
1903 ! include 'COMMON.DERIV'
1904 ! include 'COMMON.GEO'
1905 ! include 'COMMON.LOCAL'
1906 ! include 'COMMON.INTERACT'
1907 ! include 'COMMON.IOUNITS'
1908 ! include 'COMMON.NAMES'
1909 ! Revised 3/31/05 AL: correlation between random contributions to
1910 ! position and velocity increments included.
1911 real(kind=8) :: sqrt13 = 0.57735026918962576451d0 ! 1/sqrt(3)
1912 real(kind=8) :: adt,adt2
1913 integer :: i,j,ind,inres
1915 ! Add the contribution from BOTH friction and stochastic force to the
1916 ! coordinates, but ONLY the contribution from the friction forces to velocities
1919 adt=(d_a_old(j,0)+d_af_work(j))*d_time
1920 adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
1921 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1922 d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
1923 d_t(j,0)=d_t_old(j,0)+adt
1928 adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
1929 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1930 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1931 d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
1932 d_t(j,i)=d_t_old(j,i)+adt
1937 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1940 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
1941 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1942 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1943 d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
1944 d_t(j,inres)=d_t_old(j,inres)+adt
1950 end subroutine sddir_verlet1
1951 !-----------------------------------------------------------------------------
1952 subroutine sddir_verlet2
1953 ! Calculating the adjusted velocities for accelerations
1956 ! implicit real*8 (a-h,o-z)
1957 ! include 'DIMENSIONS'
1958 ! include 'COMMON.CONTROL'
1959 ! include 'COMMON.VAR'
1960 ! include 'COMMON.MD'
1962 ! include 'COMMON.LANGEVIN'
1964 ! include 'COMMON.LANGEVIN.lang0'
1966 ! include 'COMMON.CHAIN'
1967 ! include 'COMMON.DERIV'
1968 ! include 'COMMON.GEO'
1969 ! include 'COMMON.LOCAL'
1970 ! include 'COMMON.INTERACT'
1971 ! include 'COMMON.IOUNITS'
1972 ! include 'COMMON.NAMES'
1973 real(kind=8),dimension(6*nres) :: stochforcvec,d_as_work1 !(MAXRES6) maxres6=6*maxres
1974 real(kind=8) :: cos60 = 0.5d0, sin60 = 0.86602540378443864676d0
1975 integer :: i,j,ind,inres
1976 ! Revised 3/31/05 AL: correlation between random contributions to
1977 ! position and velocity increments included.
1978 ! The correlation coefficients are calculated at low-friction limit.
1979 ! Also, friction forces are now not calculated with new velocities.
1981 ! call friction_force
1982 call stochastic_force(stochforcvec)
1984 ! Compute the acceleration due to friction forces (d_af_work) and stochastic
1985 ! forces (d_as_work)
1987 call ginv_mult(stochforcvec, d_as_work1)
1993 d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j)) &
1994 +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
1999 d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j)) &
2000 +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
2005 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2008 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres) &
2009 +d_af_work(ind+j))+sin60*d_as_work(ind+j) &
2010 +cos60*d_as_work1(ind+j))*d_time
2016 end subroutine sddir_verlet2
2017 !-----------------------------------------------------------------------------
2018 subroutine max_accel
2020 ! Find the maximum difference in the accelerations of the the sites
2021 ! at the beginning and the end of the time step.
2024 ! implicit real*8 (a-h,o-z)
2025 ! include 'DIMENSIONS'
2026 ! include 'COMMON.CONTROL'
2027 ! include 'COMMON.VAR'
2028 ! include 'COMMON.MD'
2029 ! include 'COMMON.CHAIN'
2030 ! include 'COMMON.DERIV'
2031 ! include 'COMMON.GEO'
2032 ! include 'COMMON.LOCAL'
2033 ! include 'COMMON.INTERACT'
2034 ! include 'COMMON.IOUNITS'
2035 real(kind=8),dimension(3) :: aux,accel,accel_old
2036 real(kind=8) :: dacc
2040 ! aux(j)=d_a(j,0)-d_a_old(j,0)
2041 accel_old(j)=d_a_old(j,0)
2048 ! 7/3/08 changed to asymmetric difference
2050 ! accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
2051 accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
2052 accel(j)=accel(j)+0.5d0*d_a(j,i)
2053 ! if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
2054 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
2055 dacc=dabs(accel(j)-accel_old(j))
2056 ! write (iout,*) i,dacc
2057 if (dacc.gt.amax) amax=dacc
2065 accel_old(j)=d_a_old(j,0)
2070 accel_old(j)=accel_old(j)+d_a_old(j,1)
2071 accel(j)=accel(j)+d_a(j,1)
2075 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2077 ! accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
2078 accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
2079 accel(j)=accel(j)+d_a(j,i+nres)
2083 ! if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
2084 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
2085 dacc=dabs(accel(j)-accel_old(j))
2086 ! write (iout,*) "side-chain",i,dacc
2087 if (dacc.gt.amax) amax=dacc
2091 accel_old(j)=accel_old(j)+d_a_old(j,i)
2092 accel(j)=accel(j)+d_a(j,i)
2093 ! aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
2097 end subroutine max_accel
2098 !-----------------------------------------------------------------------------
2099 subroutine predict_edrift(epdrift)
2101 ! Predict the drift of the potential energy
2104 use control_data, only: lmuca
2105 ! implicit real*8 (a-h,o-z)
2106 ! include 'DIMENSIONS'
2107 ! include 'COMMON.CONTROL'
2108 ! include 'COMMON.VAR'
2109 ! include 'COMMON.MD'
2110 ! include 'COMMON.CHAIN'
2111 ! include 'COMMON.DERIV'
2112 ! include 'COMMON.GEO'
2113 ! include 'COMMON.LOCAL'
2114 ! include 'COMMON.INTERACT'
2115 ! include 'COMMON.IOUNITS'
2116 ! include 'COMMON.MUCA'
2117 real(kind=8) :: epdrift,epdriftij
2119 ! Drift of the potential energy
2125 epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
2126 if (lmuca) epdriftij=epdriftij*factor
2127 ! write (iout,*) "back",i,j,epdriftij
2128 if (epdriftij.gt.epdrift) epdrift=epdriftij
2132 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2135 dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
2136 if (lmuca) epdriftij=epdriftij*factor
2137 ! write (iout,*) "side",i,j,epdriftij
2138 if (epdriftij.gt.epdrift) epdrift=epdriftij
2142 epdrift=0.5d0*epdrift*d_time*d_time
2143 ! write (iout,*) "epdrift",epdrift
2145 end subroutine predict_edrift
2146 !-----------------------------------------------------------------------------
2147 subroutine verlet_bath
2149 ! Coupling to the thermostat by using the Berendsen algorithm
2152 ! implicit real*8 (a-h,o-z)
2153 ! include 'DIMENSIONS'
2154 ! include 'COMMON.CONTROL'
2155 ! include 'COMMON.VAR'
2156 ! include 'COMMON.MD'
2157 ! include 'COMMON.CHAIN'
2158 ! include 'COMMON.DERIV'
2159 ! include 'COMMON.GEO'
2160 ! include 'COMMON.LOCAL'
2161 ! include 'COMMON.INTERACT'
2162 ! include 'COMMON.IOUNITS'
2163 ! include 'COMMON.NAMES'
2164 real(kind=8) :: T_half,fact
2165 integer :: i,j,inres
2167 T_half=2.0d0/(dimen3*Rb)*EK
2168 fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
2169 ! write(iout,*) "T_half", T_half
2170 ! write(iout,*) "EK", EK
2171 ! write(iout,*) "fact", fact
2173 d_t(j,0)=fact*d_t(j,0)
2177 d_t(j,i)=fact*d_t(j,i)
2181 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2184 d_t(j,inres)=fact*d_t(j,inres)
2189 end subroutine verlet_bath
2190 !-----------------------------------------------------------------------------
2192 ! Set up the initial conditions of a MD simulation
2195 use control, only:tcpu
2196 !el use io_basic, only:ilen
2199 use minimm, only:minim_dc,minimize,sc_move
2200 use io_config, only:readrst
2201 use io, only:statout
2202 ! implicit real*8 (a-h,o-z)
2203 ! include 'DIMENSIONS'
2206 character(len=16) :: form
2207 integer :: IERROR,ERRCODE
2209 ! include 'COMMON.SETUP'
2210 ! include 'COMMON.CONTROL'
2211 ! include 'COMMON.VAR'
2212 ! include 'COMMON.MD'
2214 ! include 'COMMON.LANGEVIN'
2216 ! include 'COMMON.LANGEVIN.lang0'
2218 ! include 'COMMON.CHAIN'
2219 ! include 'COMMON.DERIV'
2220 ! include 'COMMON.GEO'
2221 ! include 'COMMON.LOCAL'
2222 ! include 'COMMON.INTERACT'
2223 ! include 'COMMON.IOUNITS'
2224 ! include 'COMMON.NAMES'
2225 ! include 'COMMON.REMD'
2226 real(kind=8),dimension(0:n_ene) :: energia_long,energia_short
2227 real(kind=8),dimension(3) :: vcm,incr,L
2228 real(kind=8) :: xv,sigv,lowb,highb
2229 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
2230 character(len=256) :: qstr
2233 character(len=50) :: tytul
2234 logical :: file_exist
2235 !el common /gucio/ cm
2236 integer :: i,j,ipos,iq,iw,nft_sc,iretcode,nfun,itime,ierr
2237 real(kind=8) :: etot,tt0
2241 ! write(iout,*) "d_time", d_time
2242 ! Compute the standard deviations of stochastic forces for Langevin dynamics
2243 ! if the friction coefficients do not depend on surface area
2244 if (lang.gt.0 .and. .not.surfarea) then
2246 stdforcp(i)=stdfp*dsqrt(gamp)
2249 stdforcsc(i)=stdfsc(iabs(itype(i))) &
2250 *dsqrt(gamsc(iabs(itype(i))))
2253 ! Open the pdb file for snapshotshots
2256 if (ilen(tmpdir).gt.0) &
2257 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2258 liczba(:ilen(liczba))//".pdb")
2260 file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2264 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) &
2265 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2266 liczba(:ilen(liczba))//".x")
2267 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2270 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) &
2271 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2272 liczba(:ilen(liczba))//".cx")
2273 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2279 if (ilen(tmpdir).gt.0) &
2280 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
2281 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
2283 if (ilen(tmpdir).gt.0) &
2284 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
2285 cartname=prefix(:ilen(prefix))//"_MD.cx"
2289 write (qstr,'(256(1h ))')
2292 iq = qinfrag(i,iset)*10
2293 iw = wfrag(i,iset)/100
2295 if(me.eq.king.or..not.out1file) &
2296 write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
2297 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
2302 iq = qinpair(i,iset)*10
2303 iw = wpair(i,iset)/100
2305 if(me.eq.king.or..not.out1file) &
2306 write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
2307 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
2311 ! pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
2313 ! cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
2315 ! cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
2317 ! statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
2321 if (restart1file) then
2323 inquire(file=mremd_rst_name,exist=file_exist)
2324 write (*,*) me," Before broadcast: file_exist",file_exist
2326 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,&
2329 write (*,*) me," After broadcast: file_exist",file_exist
2330 ! inquire(file=mremd_rst_name,exist=file_exist)
2331 if(me.eq.king.or..not.out1file) &
2332 write(iout,*) "Initial state read by master and distributed"
2334 if (ilen(tmpdir).gt.0) &
2335 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_' &
2336 //liczba(:ilen(liczba))//'.rst')
2337 inquire(file=rest2name,exist=file_exist)
2340 if(.not.restart1file) then
2341 if(me.eq.king.or..not.out1file) &
2342 write(iout,*) "Initial state will be read from file ",&
2343 rest2name(:ilen(rest2name))
2346 call rescale_weights(t_bath)
2348 if(me.eq.king.or..not.out1file)then
2349 if (restart1file) then
2350 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),&
2353 write(iout,*) "File ",rest2name(:ilen(rest2name)),&
2356 write(iout,*) "Initial velocities randomly generated"
2362 ! Generate initial velocities
2363 if(me.eq.king.or..not.out1file) &
2364 write(iout,*) "Initial velocities randomly generated"
2368 ! rest2name = prefix(:ilen(prefix))//'.rst'
2369 if(me.eq.king.or..not.out1file)then
2370 write (iout,*) "Initial velocities"
2372 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2373 (d_t(j,i+nres),j=1,3)
2375 ! Zeroing the total angular momentum of the system
2376 write(iout,*) "Calling the zero-angular momentum subroutine"
2379 ! Getting the potential energy and forces and velocities and accelerations
2381 ! write (iout,*) "velocity of the center of the mass:"
2382 ! write (iout,*) (vcm(j),j=1,3)
2384 d_t(j,0)=d_t(j,0)-vcm(j)
2386 ! Removing the velocity of the center of mass
2388 if(me.eq.king.or..not.out1file)then
2389 write (iout,*) "vcm right after adjustment:"
2390 write (iout,*) (vcm(j),j=1,3)
2394 if(iranconf.ne.0) then
2396 print *, 'Calling OVERLAP_SC'
2397 call overlap_sc(fail)
2400 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
2401 print *,'SC_move',nft_sc,etot
2402 if(me.eq.king.or..not.out1file) &
2403 write(iout,*) 'SC_move',nft_sc,etot
2407 print *, 'Calling MINIM_DC'
2408 call minim_dc(etot,iretcode,nfun)
2410 call geom_to_var(nvar,varia)
2411 print *,'Calling MINIMIZE.'
2412 call minimize(etot,varia,iretcode,nfun)
2413 call var_to_geom(nvar,varia)
2415 if(me.eq.king.or..not.out1file) &
2416 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2419 call chainbuild_cart
2424 kinetic_T=2.0d0/(dimen3*Rb)*EK
2425 if(me.eq.king.or..not.out1file)then
2435 call etotal(potEcomp)
2436 if (large) call enerprint(potEcomp)
2439 t_etotal=t_etotal+MPI_Wtime()-tt0
2441 t_etotal=t_etotal+tcpu()-tt0
2448 if (amax*d_time .gt. dvmax) then
2449 d_time=d_time*dvmax/amax
2450 if(me.eq.king.or..not.out1file) write (iout,*) &
2451 "Time step reduced to",d_time,&
2452 " because of too large initial acceleration."
2454 if(me.eq.king.or..not.out1file)then
2455 write(iout,*) "Potential energy and its components"
2456 call enerprint(potEcomp)
2457 ! write(iout,*) (potEcomp(i),i=0,n_ene)
2459 potE=potEcomp(0)-potEcomp(20)
2462 if (ntwe.ne.0) call statout(itime)
2463 if(me.eq.king.or..not.out1file) &
2464 write (iout,'(/a/3(a25,1pe14.5/))') "Initial:", &
2465 " Kinetic energy",EK," Potential energy",potE, &
2466 " Total energy",totE," Maximum acceleration ", &
2469 write (iout,*) "Initial coordinates"
2471 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),&
2474 write (iout,*) "Initial dC"
2476 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),&
2477 (dc(j,i+nres),j=1,3)
2479 write (iout,*) "Initial velocities"
2480 write (iout,"(13x,' backbone ',23x,' side chain')")
2482 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2483 (d_t(j,i+nres),j=1,3)
2485 write (iout,*) "Initial accelerations"
2487 ! write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2488 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),&
2489 (d_a(j,i+nres),j=1,3)
2495 d_t_old(j,i)=d_t(j,i)
2496 d_a_old(j,i)=d_a(j,i)
2498 ! write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
2507 call etotal_short(energia_short)
2508 if (large) call enerprint(potEcomp)
2511 t_eshort=t_eshort+MPI_Wtime()-tt0
2513 t_eshort=t_eshort+tcpu()-tt0
2518 if(.not.out1file .and. large) then
2519 write (iout,*) "energia_long",energia_long(0),&
2520 " energia_short",energia_short(0),&
2521 " total",energia_long(0)+energia_short(0)
2522 write (iout,*) "Initial fast-force accelerations"
2524 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
2525 (d_a(j,i+nres),j=1,3)
2528 ! 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
2531 d_a_short(j,i)=d_a(j,i)
2540 call etotal_long(energia_long)
2541 if (large) call enerprint(potEcomp)
2544 t_elong=t_elong+MPI_Wtime()-tt0
2546 t_elong=t_elong+tcpu()-tt0
2551 if(.not.out1file .and. large) then
2552 write (iout,*) "energia_long",energia_long(0)
2553 write (iout,*) "Initial slow-force accelerations"
2555 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
2556 (d_a(j,i+nres),j=1,3)
2560 t_enegrad=t_enegrad+MPI_Wtime()-tt0
2562 t_enegrad=t_enegrad+tcpu()-tt0
2566 end subroutine init_MD
2567 !-----------------------------------------------------------------------------
2568 subroutine random_vel
2570 ! implicit real*8 (a-h,o-z)
2572 use random, only:anorm_distr
2573 ! include 'DIMENSIONS'
2574 ! include 'COMMON.CONTROL'
2575 ! include 'COMMON.VAR'
2576 ! include 'COMMON.MD'
2578 ! include 'COMMON.LANGEVIN'
2580 ! include 'COMMON.LANGEVIN.lang0'
2582 ! include 'COMMON.CHAIN'
2583 ! include 'COMMON.DERIV'
2584 ! include 'COMMON.GEO'
2585 ! include 'COMMON.LOCAL'
2586 ! include 'COMMON.INTERACT'
2587 ! include 'COMMON.IOUNITS'
2588 ! include 'COMMON.NAMES'
2589 ! include 'COMMON.TIME1'
2590 real(kind=8) :: xv,sigv,lowb,highb ,Ek1
2591 integer :: i,j,ii,k,ind
2592 ! Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
2593 ! First generate velocities in the eigenspace of the G matrix
2594 ! write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
2601 sigv=dsqrt((Rb*t_bath)/geigen(i))
2604 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
2605 ! write (iout,*) "i",i," ii",ii," geigen",geigen(i),&
2606 ! " d_t_work_new",d_t_work_new(ii)
2615 ! Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
2618 ! write (iout,*) "Ek from eigenvectors",Ek1
2620 ! Transform velocities to UNRES coordinate space
2626 d_t_work(ind)=d_t_work(ind) &
2627 +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
2629 ! write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
2633 ! Transfer to the d_t vector
2635 d_t(j,0)=d_t_work(j)
2641 d_t(j,i)=d_t_work(ind)
2645 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2648 d_t(j,i+nres)=d_t_work(ind)
2653 ! write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",&
2654 ! 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
2657 end subroutine random_vel
2658 !-----------------------------------------------------------------------------
2660 subroutine sd_verlet_p_setup
2661 ! Sets up the parameters of stochastic Verlet algorithm
2662 ! implicit real*8 (a-h,o-z)
2663 ! include 'DIMENSIONS'
2664 use control, only: tcpu
2669 ! include 'COMMON.CONTROL'
2670 ! include 'COMMON.VAR'
2671 ! include 'COMMON.MD'
2673 ! include 'COMMON.LANGEVIN'
2675 ! include 'COMMON.LANGEVIN.lang0'
2677 ! include 'COMMON.CHAIN'
2678 ! include 'COMMON.DERIV'
2679 ! include 'COMMON.GEO'
2680 ! include 'COMMON.LOCAL'
2681 ! include 'COMMON.INTERACT'
2682 ! include 'COMMON.IOUNITS'
2683 ! include 'COMMON.NAMES'
2684 ! include 'COMMON.TIME1'
2685 real(kind=8),dimension(6*nres) :: emgdt !(MAXRES6) maxres6=6*maxres
2686 real(kind=8) :: pterm,vterm,rho,rhoc,vsig
2687 real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
2688 prand_vec,vrand_vec1,vrand_vec2 !(MAXRES6) maxres6=6*maxres
2689 logical :: lprn = .false.
2690 real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
2691 real(kind=8) :: ktm,gdt,egdt,gdt2,gdt3,gdt4,gdt5,gdt6,gdt7,gdt8,&
2693 integer :: i,maxres2
2700 ! AL 8/17/04 Code adapted from tinker
2702 ! Get the frictional and random terms for stochastic dynamics in the
2703 ! eigenspace of mass-scaled UNRES friction matrix
2707 gdt = fricgam(i) * d_time
2709 ! Stochastic dynamics reduces to simple MD for zero friction
2711 if (gdt .le. zero) then
2712 pfric_vec(i) = 1.0d0
2713 vfric_vec(i) = d_time
2714 afric_vec(i) = 0.5d0 * d_time * d_time
2715 prand_vec(i) = 0.0d0
2716 vrand_vec1(i) = 0.0d0
2717 vrand_vec2(i) = 0.0d0
2719 ! Analytical expressions when friction coefficient is large
2722 if (gdt .ge. gdt_radius) then
2725 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
2726 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
2727 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
2728 vterm = 1.0d0 - egdt**2
2729 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
2731 ! Use series expansions when friction coefficient is small
2742 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0 &
2743 - gdt5/120.0d0 + gdt6/720.0d0 &
2744 - gdt7/5040.0d0 + gdt8/40320.0d0 &
2745 - gdt9/362880.0d0) / fricgam(i)**2
2746 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
2747 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
2748 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0 &
2749 + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0 &
2750 + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0 &
2751 + 127.0d0*gdt9/90720.0d0
2752 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0 &
2753 - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0 &
2754 - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0 &
2755 - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
2756 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0 &
2757 - 17.0d0*gdt2/1280.0d0 &
2758 + 17.0d0*gdt3/6144.0d0 &
2759 + 40967.0d0*gdt4/34406400.0d0 &
2760 - 57203.0d0*gdt5/275251200.0d0 &
2761 - 1429487.0d0*gdt6/13212057600.0d0)
2764 ! Compute the scaling factors of random terms for the nonzero friction case
2766 ktm = 0.5d0*d_time/fricgam(i)
2767 psig = dsqrt(ktm*pterm) / fricgam(i)
2768 vsig = dsqrt(ktm*vterm)
2769 rhoc = dsqrt(1.0d0 - rho*rho)
2771 vrand_vec1(i) = vsig * rho
2772 vrand_vec2(i) = vsig * rhoc
2777 "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",&
2780 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),&
2781 afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
2785 ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
2788 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
2789 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
2790 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
2791 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
2792 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
2793 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
2796 t_sdsetup=t_sdsetup+MPI_Wtime()
2798 t_sdsetup=t_sdsetup+tcpu()-tt0
2801 end subroutine sd_verlet_p_setup
2802 !-----------------------------------------------------------------------------
2803 subroutine eigtransf1(n,ndim,ab,d,c)
2807 real(kind=8) :: ab(ndim,ndim,n),c(ndim,n),d(ndim)
2813 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
2818 end subroutine eigtransf1
2819 !-----------------------------------------------------------------------------
2820 subroutine eigtransf(n,ndim,a,b,d,c)
2824 real(kind=8) :: a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
2830 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
2835 end subroutine eigtransf
2836 !-----------------------------------------------------------------------------
2837 subroutine sd_verlet1
2839 ! Applying stochastic velocity Verlet algorithm - step 1 to velocities
2841 ! implicit real*8 (a-h,o-z)
2842 ! include 'DIMENSIONS'
2843 ! include 'COMMON.CONTROL'
2844 ! include 'COMMON.VAR'
2845 ! include 'COMMON.MD'
2847 ! include 'COMMON.LANGEVIN'
2849 ! include 'COMMON.LANGEVIN.lang0'
2851 ! include 'COMMON.CHAIN'
2852 ! include 'COMMON.DERIV'
2853 ! include 'COMMON.GEO'
2854 ! include 'COMMON.LOCAL'
2855 ! include 'COMMON.INTERACT'
2856 ! include 'COMMON.IOUNITS'
2857 ! include 'COMMON.NAMES'
2858 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
2859 !el common /stochcalc/ stochforcvec
2860 logical :: lprn = .false.
2861 real(kind=8) :: ddt1,ddt2
2862 integer :: i,j,ind,inres
2864 ! write (iout,*) "dc_old"
2866 ! write (iout,'(i5,3f10.5,5x,3f10.5)')
2867 ! & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
2870 dc_work(j)=dc_old(j,0)
2871 d_t_work(j)=d_t_old(j,0)
2872 d_a_work(j)=d_a_old(j,0)
2877 dc_work(ind+j)=dc_old(j,i)
2878 d_t_work(ind+j)=d_t_old(j,i)
2879 d_a_work(ind+j)=d_a_old(j,i)
2884 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2886 dc_work(ind+j)=dc_old(j,i+nres)
2887 d_t_work(ind+j)=d_t_old(j,i+nres)
2888 d_a_work(ind+j)=d_a_old(j,i+nres)
2896 "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",&
2900 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
2901 vfric_mat(i,j),afric_mat(i,j),&
2902 prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
2910 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) &
2911 +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
2912 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
2913 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
2915 d_t_work_new(i)=ddt1+0.5d0*ddt2
2916 d_t_work(i)=ddt1+ddt2
2921 d_t(j,0)=d_t_work(j)
2926 dc(j,i)=dc_work(ind+j)
2927 d_t(j,i)=d_t_work(ind+j)
2932 if (itype(i).ne.10) then
2935 dc(j,inres)=dc_work(ind+j)
2936 d_t(j,inres)=d_t_work(ind+j)
2942 end subroutine sd_verlet1
2943 !-----------------------------------------------------------------------------
2944 subroutine sd_verlet2
2946 ! Calculating the adjusted velocities for accelerations
2948 ! implicit real*8 (a-h,o-z)
2949 ! include 'DIMENSIONS'
2950 ! include 'COMMON.CONTROL'
2951 ! include 'COMMON.VAR'
2952 ! include 'COMMON.MD'
2954 ! include 'COMMON.LANGEVIN'
2956 ! include 'COMMON.LANGEVIN.lang0'
2958 ! include 'COMMON.CHAIN'
2959 ! include 'COMMON.DERIV'
2960 ! include 'COMMON.GEO'
2961 ! include 'COMMON.LOCAL'
2962 ! include 'COMMON.INTERACT'
2963 ! include 'COMMON.IOUNITS'
2964 ! include 'COMMON.NAMES'
2965 !el real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV !(MAXRES6) maxres6=6*maxres
2966 real(kind=8),dimension(6*nres) :: stochforcvecV !(MAXRES6) maxres6=6*maxres
2967 !el common /stochcalc/ stochforcvec
2969 real(kind=8) :: ddt1,ddt2
2970 integer :: i,j,ind,inres
2971 ! Compute the stochastic forces which contribute to velocity change
2973 call stochastic_force(stochforcvecV)
2980 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
2981 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+ &
2982 vrand_mat2(i,j)*stochforcvecV(j)
2984 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
2988 d_t(j,0)=d_t_work(j)
2993 d_t(j,i)=d_t_work(ind+j)
2998 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3001 d_t(j,inres)=d_t_work(ind+j)
3007 end subroutine sd_verlet2
3008 !-----------------------------------------------------------------------------
3009 subroutine sd_verlet_ciccotti_setup
3011 ! Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
3013 ! implicit real*8 (a-h,o-z)
3014 ! include 'DIMENSIONS'
3015 use control, only: tcpu
3020 ! include 'COMMON.CONTROL'
3021 ! include 'COMMON.VAR'
3022 ! include 'COMMON.MD'
3024 ! include 'COMMON.LANGEVIN'
3026 ! include 'COMMON.LANGEVIN.lang0'
3028 ! include 'COMMON.CHAIN'
3029 ! include 'COMMON.DERIV'
3030 ! include 'COMMON.GEO'
3031 ! include 'COMMON.LOCAL'
3032 ! include 'COMMON.INTERACT'
3033 ! include 'COMMON.IOUNITS'
3034 ! include 'COMMON.NAMES'
3035 ! include 'COMMON.TIME1'
3036 real(kind=8),dimension(6*nres) :: emgdt !(MAXRES6) maxres6=6*maxres
3037 real(kind=8) :: pterm,vterm,rho,rhoc,vsig
3038 real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
3039 prand_vec,vrand_vec1,vrand_vec2 !(MAXRES6) maxres6=6*maxres
3040 logical :: lprn = .false.
3041 real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
3042 real(kind=8) :: ktm,gdt,egdt,tt0
3043 integer :: i,maxres2
3050 ! AL 8/17/04 Code adapted from tinker
3052 ! Get the frictional and random terms for stochastic dynamics in the
3053 ! eigenspace of mass-scaled UNRES friction matrix
3057 write (iout,*) "i",i," fricgam",fricgam(i)
3058 gdt = fricgam(i) * d_time
3060 ! Stochastic dynamics reduces to simple MD for zero friction
3062 if (gdt .le. zero) then
3063 pfric_vec(i) = 1.0d0
3064 vfric_vec(i) = d_time
3065 afric_vec(i) = 0.5d0*d_time*d_time
3066 prand_vec(i) = afric_vec(i)
3067 vrand_vec2(i) = vfric_vec(i)
3069 ! Analytical expressions when friction coefficient is large
3074 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
3075 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
3076 prand_vec(i) = afric_vec(i)
3077 vrand_vec2(i) = vfric_vec(i)
3079 ! Compute the scaling factors of random terms for the nonzero friction case
3081 ! ktm = 0.5d0*d_time/fricgam(i)
3082 ! psig = dsqrt(ktm*pterm) / fricgam(i)
3083 ! vsig = dsqrt(ktm*vterm)
3084 ! prand_vec(i) = psig*afric_vec(i)
3085 ! vrand_vec2(i) = vsig*vfric_vec(i)
3090 "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",&
3093 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),&
3094 afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
3098 ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
3100 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
3101 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
3102 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
3103 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
3104 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
3106 t_sdsetup=t_sdsetup+MPI_Wtime()
3108 t_sdsetup=t_sdsetup+tcpu()-tt0
3111 end subroutine sd_verlet_ciccotti_setup
3112 !-----------------------------------------------------------------------------
3113 subroutine sd_verlet1_ciccotti
3115 ! Applying stochastic velocity Verlet algorithm - step 1 to velocities
3116 ! implicit real*8 (a-h,o-z)
3118 ! include 'DIMENSIONS'
3122 ! include 'COMMON.CONTROL'
3123 ! include 'COMMON.VAR'
3124 ! include 'COMMON.MD'
3126 ! include 'COMMON.LANGEVIN'
3128 ! include 'COMMON.LANGEVIN.lang0'
3130 ! include 'COMMON.CHAIN'
3131 ! include 'COMMON.DERIV'
3132 ! include 'COMMON.GEO'
3133 ! include 'COMMON.LOCAL'
3134 ! include 'COMMON.INTERACT'
3135 ! include 'COMMON.IOUNITS'
3136 ! include 'COMMON.NAMES'
3137 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
3138 !el common /stochcalc/ stochforcvec
3139 logical :: lprn = .false.
3140 real(kind=8) :: ddt1,ddt2
3141 integer :: i,j,ind,inres
3142 ! write (iout,*) "dc_old"
3144 ! write (iout,'(i5,3f10.5,5x,3f10.5)')
3145 ! & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3148 dc_work(j)=dc_old(j,0)
3149 d_t_work(j)=d_t_old(j,0)
3150 d_a_work(j)=d_a_old(j,0)
3155 dc_work(ind+j)=dc_old(j,i)
3156 d_t_work(ind+j)=d_t_old(j,i)
3157 d_a_work(ind+j)=d_a_old(j,i)
3162 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3164 dc_work(ind+j)=dc_old(j,i+nres)
3165 d_t_work(ind+j)=d_t_old(j,i+nres)
3166 d_a_work(ind+j)=d_a_old(j,i+nres)
3175 "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",&
3179 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
3180 vfric_mat(i,j),afric_mat(i,j),&
3181 prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3189 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) &
3190 +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3191 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3192 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3194 d_t_work_new(i)=ddt1+0.5d0*ddt2
3195 d_t_work(i)=ddt1+ddt2
3200 d_t(j,0)=d_t_work(j)
3205 dc(j,i)=dc_work(ind+j)
3206 d_t(j,i)=d_t_work(ind+j)
3211 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3214 dc(j,inres)=dc_work(ind+j)
3215 d_t(j,inres)=d_t_work(ind+j)
3221 end subroutine sd_verlet1_ciccotti
3222 !-----------------------------------------------------------------------------
3223 subroutine sd_verlet2_ciccotti
3225 ! Calculating the adjusted velocities for accelerations
3227 ! implicit real*8 (a-h,o-z)
3228 ! include 'DIMENSIONS'
3229 ! include 'COMMON.CONTROL'
3230 ! include 'COMMON.VAR'
3231 ! include 'COMMON.MD'
3233 ! include 'COMMON.LANGEVIN'
3235 ! include 'COMMON.LANGEVIN.lang0'
3237 ! include 'COMMON.CHAIN'
3238 ! include 'COMMON.DERIV'
3239 ! include 'COMMON.GEO'
3240 ! include 'COMMON.LOCAL'
3241 ! include 'COMMON.INTERACT'
3242 ! include 'COMMON.IOUNITS'
3243 ! include 'COMMON.NAMES'
3244 !el real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV !(MAXRES6) maxres6=6*maxres
3245 real(kind=8),dimension(6*nres) :: stochforcvecV !(MAXRES6) maxres6=6*maxres
3246 !el common /stochcalc/ stochforcvec
3247 real(kind=8) :: ddt1,ddt2
3248 integer :: i,j,ind,inres
3250 ! Compute the stochastic forces which contribute to velocity change
3252 call stochastic_force(stochforcvecV)
3259 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3260 ! ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
3261 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
3263 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
3267 d_t(j,0)=d_t_work(j)
3272 d_t(j,i)=d_t_work(ind+j)
3277 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3280 d_t(j,inres)=d_t_work(ind+j)
3286 end subroutine sd_verlet2_ciccotti
3288 !-----------------------------------------------------------------------------
3290 !-----------------------------------------------------------------------------
3291 subroutine inertia_tensor
3293 ! Calculating the intertia tensor for the entire protein in order to
3294 ! remove the perpendicular components of velocity matrix which cause
3295 ! the molecule to rotate.
3298 ! implicit real*8 (a-h,o-z)
3299 ! include 'DIMENSIONS'
3300 ! include 'COMMON.CONTROL'
3301 ! include 'COMMON.VAR'
3302 ! include 'COMMON.MD'
3303 ! include 'COMMON.CHAIN'
3304 ! include 'COMMON.DERIV'
3305 ! include 'COMMON.GEO'
3306 ! include 'COMMON.LOCAL'
3307 ! include 'COMMON.INTERACT'
3308 ! include 'COMMON.IOUNITS'
3309 ! include 'COMMON.NAMES'
3311 real(kind=8),dimension(3,3) :: Im,Imcp,eigvec,Id
3312 real(kind=8),dimension(3) :: pr,eigval,L,vp,vrot
3313 real(kind=8) :: M_SC,mag,mag2
3314 real(kind=8),dimension(3,0:nres) :: vpp !(3,0:MAXRES)
3315 real(kind=8),dimension(3) :: vs_p,pp,incr,v
3316 real(kind=8),dimension(3,3) :: pr1,pr2
3318 !el common /gucio/ cm
3319 integer :: iti,inres,i,j,k
3330 ! calculating the center of the mass of the protein
3333 cm(j)=cm(j)+c(j,i)+0.5d0*dc(j,i)
3342 M_SC=M_SC+msc(iabs(iti))
3345 cm(j)=cm(j)+msc(iabs(iti))*c(j,inres)
3349 cm(j)=cm(j)/(M_SC+(nct-nnt)*mp)
3354 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
3356 Im(1,1)=Im(1,1)+mp*(pr(2)*pr(2)+pr(3)*pr(3))
3357 Im(1,2)=Im(1,2)-mp*pr(1)*pr(2)
3358 Im(1,3)=Im(1,3)-mp*pr(1)*pr(3)
3359 Im(2,3)=Im(2,3)-mp*pr(2)*pr(3)
3360 Im(2,2)=Im(2,2)+mp*(pr(3)*pr(3)+pr(1)*pr(1))
3361 Im(3,3)=Im(3,3)+mp*(pr(1)*pr(1)+pr(2)*pr(2))
3368 pr(j)=c(j,inres)-cm(j)
3370 Im(1,1)=Im(1,1)+msc(iabs(iti))*(pr(2)*pr(2)+pr(3)*pr(3))
3371 Im(1,2)=Im(1,2)-msc(iabs(iti))*pr(1)*pr(2)
3372 Im(1,3)=Im(1,3)-msc(iabs(iti))*pr(1)*pr(3)
3373 Im(2,3)=Im(2,3)-msc(iabs(iti))*pr(2)*pr(3)
3374 Im(2,2)=Im(2,2)+msc(iabs(iti))*(pr(3)*pr(3)+pr(1)*pr(1))
3375 Im(3,3)=Im(3,3)+msc(iabs(iti))*(pr(1)*pr(1)+pr(2)*pr(2))
3379 Im(1,1)=Im(1,1)+Ip*(1-dc_norm(1,i)*dc_norm(1,i))* &
3380 vbld(i+1)*vbld(i+1)*0.25d0
3381 Im(1,2)=Im(1,2)+Ip*(-dc_norm(1,i)*dc_norm(2,i))* &
3382 vbld(i+1)*vbld(i+1)*0.25d0
3383 Im(1,3)=Im(1,3)+Ip*(-dc_norm(1,i)*dc_norm(3,i))* &
3384 vbld(i+1)*vbld(i+1)*0.25d0
3385 Im(2,3)=Im(2,3)+Ip*(-dc_norm(2,i)*dc_norm(3,i))* &
3386 vbld(i+1)*vbld(i+1)*0.25d0
3387 Im(2,2)=Im(2,2)+Ip*(1-dc_norm(2,i)*dc_norm(2,i))* &
3388 vbld(i+1)*vbld(i+1)*0.25d0
3389 Im(3,3)=Im(3,3)+Ip*(1-dc_norm(3,i)*dc_norm(3,i))* &
3390 vbld(i+1)*vbld(i+1)*0.25d0
3395 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3398 Im(1,1)=Im(1,1)+Isc(iti)*(1-dc_norm(1,inres)* &
3399 dc_norm(1,inres))*vbld(inres)*vbld(inres)
3400 Im(1,2)=Im(1,2)-Isc(iti)*(dc_norm(1,inres)* &
3401 dc_norm(2,inres))*vbld(inres)*vbld(inres)
3402 Im(1,3)=Im(1,3)-Isc(iti)*(dc_norm(1,inres)* &
3403 dc_norm(3,inres))*vbld(inres)*vbld(inres)
3404 Im(2,3)=Im(2,3)-Isc(iti)*(dc_norm(2,inres)* &
3405 dc_norm(3,inres))*vbld(inres)*vbld(inres)
3406 Im(2,2)=Im(2,2)+Isc(iti)*(1-dc_norm(2,inres)* &
3407 dc_norm(2,inres))*vbld(inres)*vbld(inres)
3408 Im(3,3)=Im(3,3)+Isc(iti)*(1-dc_norm(3,inres)* &
3409 dc_norm(3,inres))*vbld(inres)*vbld(inres)
3414 ! write(iout,*) "The angular momentum before adjustment:"
3415 ! write(iout,*) (L(j),j=1,3)
3421 ! Copying the Im matrix for the djacob subroutine
3429 ! Finding the eigenvectors and eignvalues of the inertia tensor
3430 call djacob(3,3,10000,1.0d-10,Imcp,eigvec,eigval)
3431 ! write (iout,*) "Eigenvalues & Eigenvectors"
3432 ! write (iout,'(5x,3f10.5)') (eigval(i),i=1,3)
3435 ! write (iout,'(i5,3f10.5)') i,(eigvec(i,j),j=1,3)
3437 ! Constructing the diagonalized matrix
3439 if (dabs(eigval(i)).gt.1.0d-15) then
3440 Id(i,i)=1.0d0/eigval(i)
3447 Imcp(i,j)=eigvec(j,i)
3453 pr1(i,j)=pr1(i,j)+Id(i,k)*Imcp(k,j)
3460 pr2(i,j)=pr2(i,j)+eigvec(i,k)*pr1(k,j)
3464 ! Calculating the total rotational velocity of the molecule
3467 vrot(i)=vrot(i)+pr2(i,j)*L(j)
3470 ! Resetting the velocities
3472 call vecpr(vrot(1),dc(1,i),vp)
3474 d_t(j,i)=d_t(j,i)-vp(j)
3478 if(itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3480 call vecpr(vrot(1),dc(1,inres),vp)
3482 d_t(j,inres)=d_t(j,inres)-vp(j)
3487 ! write(iout,*) "The angular momentum after adjustment:"
3488 ! write(iout,*) (L(j),j=1,3)
3491 end subroutine inertia_tensor
3492 !-----------------------------------------------------------------------------
3493 subroutine angmom(cm,L)
3496 ! implicit real*8 (a-h,o-z)
3497 ! include 'DIMENSIONS'
3498 ! include 'COMMON.CONTROL'
3499 ! include 'COMMON.VAR'
3500 ! include 'COMMON.MD'
3501 ! include 'COMMON.CHAIN'
3502 ! include 'COMMON.DERIV'
3503 ! include 'COMMON.GEO'
3504 ! include 'COMMON.LOCAL'
3505 ! include 'COMMON.INTERACT'
3506 ! include 'COMMON.IOUNITS'
3507 ! include 'COMMON.NAMES'
3509 real(kind=8),dimension(3) :: L,cm,pr,vp,vrot,incr,v,pp
3510 integer :: iti,inres,i,j
3511 ! Calculate the angular momentum
3520 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
3523 v(j)=incr(j)+0.5d0*d_t(j,i)
3526 incr(j)=incr(j)+d_t(j,i)
3528 call vecpr(pr(1),v(1),vp)
3534 pp(j)=0.5d0*d_t(j,i)
3536 call vecpr(pr(1),pp(1),vp)
3548 pr(j)=c(j,inres)-cm(j)
3550 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3552 v(j)=incr(j)+d_t(j,inres)
3559 call vecpr(pr(1),v(1),vp)
3560 ! write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),
3561 ! & " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
3563 L(j)=L(j)+msc(iabs(iti))*vp(j)
3565 ! write (iout,*) "L",(l(j),j=1,3)
3566 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3568 v(j)=incr(j)+d_t(j,inres)
3570 call vecpr(dc(1,inres),d_t(1,inres),vp)
3572 L(j)=L(j)+Isc(iti)*vp(j)
3576 incr(j)=incr(j)+d_t(j,i)
3580 end subroutine angmom
3581 !-----------------------------------------------------------------------------
3582 subroutine vcm_vel(vcm)
3585 ! implicit real*8 (a-h,o-z)
3586 ! include 'DIMENSIONS'
3587 ! include 'COMMON.VAR'
3588 ! include 'COMMON.MD'
3589 ! include 'COMMON.CHAIN'
3590 ! include 'COMMON.DERIV'
3591 ! include 'COMMON.GEO'
3592 ! include 'COMMON.LOCAL'
3593 ! include 'COMMON.INTERACT'
3594 ! include 'COMMON.IOUNITS'
3595 real(kind=8),dimension(3) :: vcm,vv
3596 real(kind=8) :: summas,amas
3608 vcm(j)=vcm(j)+mp*(vv(j)+0.5d0*d_t(j,i))
3611 amas=msc(iabs(itype(i)))
3613 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3615 vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres))
3619 vcm(j)=vcm(j)+amas*vv(j)
3623 vv(j)=vv(j)+d_t(j,i)
3626 ! write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas
3628 vcm(j)=vcm(j)/summas
3631 end subroutine vcm_vel
3632 !-----------------------------------------------------------------------------
3634 !-----------------------------------------------------------------------------
3636 ! RATTLE algorithm for velocity Verlet - step 1, UNRES
3640 ! implicit real*8 (a-h,o-z)
3641 ! include 'DIMENSIONS'
3643 ! include 'COMMON.CONTROL'
3644 ! include 'COMMON.VAR'
3645 ! include 'COMMON.MD'
3647 ! include 'COMMON.LANGEVIN'
3649 ! include 'COMMON.LANGEVIN.lang0'
3651 ! include 'COMMON.CHAIN'
3652 ! include 'COMMON.DERIV'
3653 ! include 'COMMON.GEO'
3654 ! include 'COMMON.LOCAL'
3655 ! include 'COMMON.INTERACT'
3656 ! include 'COMMON.IOUNITS'
3657 ! include 'COMMON.NAMES'
3658 ! include 'COMMON.TIME1'
3659 !el real(kind=8) :: gginv(2*nres,2*nres),&
3660 !el gdc(3,2*nres,2*nres)
3661 real(kind=8) :: dC_uncor(3,2*nres) !,&
3662 !el real(kind=8) :: Cmat(2*nres,2*nres)
3663 real(kind=8) :: x(2*nres),xcorr(3,2*nres) !maxres2=2*maxres
3664 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
3665 !el common /przechowalnia/ nbond
3666 integer :: max_rattle = 5
3667 logical :: lprn = .false., lprn1 = .false., not_done
3668 real(kind=8) :: tol_rattle = 1.0d-5
3670 integer :: ii,i,j,jj,l,ind,ind1,nres2
3673 !el /common/ przechowalnia
3675 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
3676 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
3677 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
3679 if (lprn) write (iout,*) "RATTLE1"
3682 if (itype(i).ne.10) nbond=nbond+1
3684 ! Make a folded form of the Ginv-matrix
3697 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
3701 if (itype(k).ne.10) then
3705 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
3712 if (itype(i).ne.10) then
3722 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
3726 if (itype(k).ne.10) then
3730 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
3738 write (iout,*) "Matrix GGinv"
3739 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
3745 if (iter.gt.max_rattle) then
3746 write (iout,*) "Error - too many iterations in RATTLE."
3749 ! Calculate the matrix C = GG**(-1) dC_old o dC
3754 dC_uncor(j,ind1)=dC(j,i)
3758 if (itype(i).ne.10) then
3761 dC_uncor(j,ind1)=dC(j,i+nres)
3770 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
3774 if (itype(k).ne.10) then
3777 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
3782 ! Calculate deviations from standard virtual-bond lengths
3786 x(ind)=vbld(i+1)**2-vbl**2
3789 if (itype(i).ne.10) then
3791 x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i))**2
3795 write (iout,*) "Coordinates and violations"
3797 write(iout,'(i5,3f10.5,5x,e15.5)') &
3798 i,(dC_uncor(j,i),j=1,3),x(i)
3800 write (iout,*) "Velocities and violations"
3804 write (iout,'(2i5,3f10.5,5x,e15.5)') &
3805 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
3808 if (itype(i).ne.10) then
3810 write (iout,'(2i5,3f10.5,5x,e15.5)') &
3811 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
3812 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
3815 ! write (iout,*) "gdc"
3817 ! write (iout,*) "i",i
3819 ! write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
3825 if (dabs(x(i)).gt.xmax) then
3829 if (xmax.lt.tol_rattle) then
3833 ! Calculate the matrix of the system of equations
3838 Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
3843 write (iout,*) "Matrix Cmat"
3844 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
3846 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
3847 ! Add constraint term to positions
3854 xx = xx+x(ii)*gdc(j,ind,ii)
3858 d_t_new(j,i)=d_t_new(j,i)-xx/d_time
3862 if (itype(i).ne.10) then
3867 xx = xx+x(ii)*gdc(j,ind,ii)
3870 dC(j,i+nres)=dC(j,i+nres)-xx
3871 d_t_new(j,i+nres)=d_t_new(j,i+nres)-xx/d_time
3875 ! Rebuild the chain using the new coordinates
3876 call chainbuild_cart
3878 write (iout,*) "New coordinates, Lagrange multipliers,",&
3879 " and differences between actual and standard bond lengths"
3883 xx=vbld(i+1)**2-vbl**2
3884 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
3885 i,(dC(j,i),j=1,3),x(ind),xx
3888 if (itype(i).ne.10) then
3890 xx=vbld(i+nres)**2-vbldsc0(1,itype(i))**2
3891 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
3892 i,(dC(j,i+nres),j=1,3),x(ind),xx
3895 write (iout,*) "Velocities and violations"
3899 write (iout,'(2i5,3f10.5,5x,e15.5)') &
3900 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
3903 if (itype(i).ne.10) then
3905 write (iout,'(2i5,3f10.5,5x,e15.5)') &
3906 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
3907 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
3914 10 write (iout,*) "Error - singularity in solving the system",&
3915 " of equations for Lagrange multipliers."
3919 "RATTLE inactive; use -DRATTLE switch at compile time."
3922 end subroutine rattle1
3923 !-----------------------------------------------------------------------------
3925 ! RATTLE algorithm for velocity Verlet - step 2, UNRES
3929 ! implicit real*8 (a-h,o-z)
3930 ! include 'DIMENSIONS'
3932 ! include 'COMMON.CONTROL'
3933 ! include 'COMMON.VAR'
3934 ! include 'COMMON.MD'
3936 ! include 'COMMON.LANGEVIN'
3938 ! include 'COMMON.LANGEVIN.lang0'
3940 ! include 'COMMON.CHAIN'
3941 ! include 'COMMON.DERIV'
3942 ! include 'COMMON.GEO'
3943 ! include 'COMMON.LOCAL'
3944 ! include 'COMMON.INTERACT'
3945 ! include 'COMMON.IOUNITS'
3946 ! include 'COMMON.NAMES'
3947 ! include 'COMMON.TIME1'
3948 !el real(kind=8) :: gginv(2*nres,2*nres),&
3949 !el gdc(3,2*nres,2*nres)
3950 real(kind=8) :: dC_uncor(3,2*nres) !,&
3951 !el Cmat(2*nres,2*nres)
3952 real(kind=8) :: x(2*nres) !maxres2=2*maxres
3953 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
3954 !el common /przechowalnia/ nbond
3955 integer :: max_rattle = 5
3956 logical :: lprn = .false., lprn1 = .false., not_done
3957 real(kind=8) :: tol_rattle = 1.0d-5
3961 !el /common/ przechowalnia
3962 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
3963 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
3964 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
3966 if (lprn) write (iout,*) "RATTLE2"
3967 if (lprn) write (iout,*) "Velocity correction"
3968 ! Calculate the matrix G dC
3974 gdc(j,i,ind)=GGinv(i,ind)*dC(j,k)
3978 if (itype(k).ne.10) then
3981 gdc(j,i,ind)=GGinv(i,ind)*dC(j,k+nres)
3987 ! write (iout,*) "gdc"
3989 ! write (iout,*) "i",i
3991 ! write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
3995 ! Calculate the matrix of the system of equations
4002 Cmat(ind,j)=Cmat(ind,j)+dC(k,i)*gdc(k,ind,j)
4007 if (itype(i).ne.10) then
4012 Cmat(ind,j)=Cmat(ind,j)+dC(k,i+nres)*gdc(k,ind,j)
4017 ! Calculate the scalar product dC o d_t_new
4021 x(ind)=scalar(d_t(1,i),dC(1,i))
4024 if (itype(i).ne.10) then
4026 x(ind)=scalar(d_t(1,i+nres),dC(1,i+nres))
4030 write (iout,*) "Velocities and violations"
4034 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4035 i,ind,(d_t(j,i),j=1,3),x(ind)
4038 if (itype(i).ne.10) then
4040 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4041 i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind)
4047 if (dabs(x(i)).gt.xmax) then
4051 if (xmax.lt.tol_rattle) then
4056 write (iout,*) "Matrix Cmat"
4057 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4059 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4060 ! Add constraint term to velocities
4067 xx = xx+x(ii)*gdc(j,ind,ii)
4069 d_t(j,i)=d_t(j,i)-xx
4073 if (itype(i).ne.10) then
4078 xx = xx+x(ii)*gdc(j,ind,ii)
4080 d_t(j,i+nres)=d_t(j,i+nres)-xx
4086 "New velocities, Lagrange multipliers violations"
4090 if (lprn) write (iout,'(2i5,3f10.5,5x,2e15.5)') &
4091 i,ind,(d_t(j,i),j=1,3),x(ind),scalar(d_t(1,i),dC(1,i))
4094 if (itype(i).ne.10) then
4096 write (iout,'(2i5,3f10.5,5x,2e15.5)') &
4097 i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind),&
4098 scalar(d_t(1,i+nres),dC(1,i+nres))
4104 10 write (iout,*) "Error - singularity in solving the system",&
4105 " of equations for Lagrange multipliers."
4109 "RATTLE inactive; use -DRATTLE option at compile time."
4112 end subroutine rattle2
4113 !-----------------------------------------------------------------------------
4114 subroutine rattle_brown
4115 ! RATTLE/LINCS algorithm for Brownian dynamics, UNRES
4119 ! implicit real*8 (a-h,o-z)
4120 ! include 'DIMENSIONS'
4122 ! include 'COMMON.CONTROL'
4123 ! include 'COMMON.VAR'
4124 ! include 'COMMON.MD'
4126 ! include 'COMMON.LANGEVIN'
4128 ! include 'COMMON.LANGEVIN.lang0'
4130 ! include 'COMMON.CHAIN'
4131 ! include 'COMMON.DERIV'
4132 ! include 'COMMON.GEO'
4133 ! include 'COMMON.LOCAL'
4134 ! include 'COMMON.INTERACT'
4135 ! include 'COMMON.IOUNITS'
4136 ! include 'COMMON.NAMES'
4137 ! include 'COMMON.TIME1'
4138 !el real(kind=8) :: gginv(2*nres,2*nres),&
4139 !el gdc(3,2*nres,2*nres)
4140 real(kind=8) :: dC_uncor(3,2*nres) !,&
4141 !el real(kind=8) :: Cmat(2*nres,2*nres)
4142 real(kind=8) :: x(2*nres) !maxres2=2*maxres
4143 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
4144 !el common /przechowalnia/ nbond
4145 integer :: max_rattle = 5
4146 logical :: lprn = .true., lprn1 = .true., not_done
4147 real(kind=8) :: tol_rattle = 1.0d-5
4151 !el /common/ przechowalnia
4152 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4153 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4154 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4157 if (lprn) write (iout,*) "RATTLE_BROWN"
4160 if (itype(i).ne.10) nbond=nbond+1
4162 ! Make a folded form of the Ginv-matrix
4175 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4179 if (itype(k).ne.10) then
4183 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4190 if (itype(i).ne.10) then
4200 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4204 if (itype(k).ne.10) then
4208 if (j.eq.1 .and. l.eq.1)GGinv(ii,jj)=fricmat(ind,ind1)
4216 write (iout,*) "Matrix GGinv"
4217 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
4223 if (iter.gt.max_rattle) then
4224 write (iout,*) "Error - too many iterations in RATTLE."
4227 ! Calculate the matrix C = GG**(-1) dC_old o dC
4232 dC_uncor(j,ind1)=dC(j,i)
4236 if (itype(i).ne.10) then
4239 dC_uncor(j,ind1)=dC(j,i+nres)
4248 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
4252 if (itype(k).ne.10) then
4255 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
4260 ! Calculate deviations from standard virtual-bond lengths
4264 x(ind)=vbld(i+1)**2-vbl**2
4267 if (itype(i).ne.10) then
4269 x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i))**2
4273 write (iout,*) "Coordinates and violations"
4275 write(iout,'(i5,3f10.5,5x,e15.5)') &
4276 i,(dC_uncor(j,i),j=1,3),x(i)
4278 write (iout,*) "Velocities and violations"
4282 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4283 i,ind,(d_t(j,i),j=1,3),scalar(d_t(1,i),dC_old(1,i))
4286 if (itype(i).ne.10) then
4288 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4289 i+nres,ind,(d_t(j,i+nres),j=1,3),&
4290 scalar(d_t(1,i+nres),dC_old(1,i+nres))
4293 write (iout,*) "gdc"
4295 write (iout,*) "i",i
4297 write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4303 if (dabs(x(i)).gt.xmax) then
4307 if (xmax.lt.tol_rattle) then
4311 ! Calculate the matrix of the system of equations
4316 Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
4321 write (iout,*) "Matrix Cmat"
4322 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4324 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4325 ! Add constraint term to positions
4332 xx = xx+x(ii)*gdc(j,ind,ii)
4335 d_t(j,i)=d_t(j,i)+xx/d_time
4340 if (itype(i).ne.10) then
4345 xx = xx+x(ii)*gdc(j,ind,ii)
4348 d_t(j,i+nres)=d_t(j,i+nres)+xx/d_time
4349 dC(j,i+nres)=dC(j,i+nres)+xx
4353 ! Rebuild the chain using the new coordinates
4354 call chainbuild_cart
4356 write (iout,*) "New coordinates, Lagrange multipliers,",&
4357 " and differences between actual and standard bond lengths"
4361 xx=vbld(i+1)**2-vbl**2
4362 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4363 i,(dC(j,i),j=1,3),x(ind),xx
4366 if (itype(i).ne.10) then
4368 xx=vbld(i+nres)**2-vbldsc0(1,itype(i))**2
4369 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4370 i,(dC(j,i+nres),j=1,3),x(ind),xx
4373 write (iout,*) "Velocities and violations"
4377 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4378 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4381 if (itype(i).ne.10) then
4383 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4384 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4385 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4392 10 write (iout,*) "Error - singularity in solving the system",&
4393 " of equations for Lagrange multipliers."
4397 "RATTLE inactive; use -DRATTLE option at compile time"
4400 end subroutine rattle_brown
4401 !-----------------------------------------------------------------------------
4403 !-----------------------------------------------------------------------------
4404 subroutine friction_force
4409 ! implicit real*8 (a-h,o-z)
4410 ! include 'DIMENSIONS'
4411 ! include 'COMMON.VAR'
4412 ! include 'COMMON.CHAIN'
4413 ! include 'COMMON.DERIV'
4414 ! include 'COMMON.GEO'
4415 ! include 'COMMON.LOCAL'
4416 ! include 'COMMON.INTERACT'
4417 ! include 'COMMON.MD'
4419 ! include 'COMMON.LANGEVIN'
4421 ! include 'COMMON.LANGEVIN.lang0'
4423 ! include 'COMMON.IOUNITS'
4424 !el real(kind=8),dimension(6*nres) :: gamvec !(MAXRES6) maxres6=6*maxres
4425 !el common /syfek/ gamvec
4426 real(kind=8) :: vv(3),vvtot(3,nres),v_work(6*nres) !,&
4427 !el ginvfric(2*nres,2*nres) !maxres2=2*maxres
4428 !el common /przechowalnia/ ginvfric
4430 logical :: lprn = .false., checkmode = .false.
4431 integer :: i,j,ind,k,nres2,nres6
4435 if(.not.allocated(gamvec)) allocate(gamvec(nres6)) !(MAXRES6)
4436 if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
4444 d_t_work(j)=d_t(j,0)
4449 d_t_work(ind+j)=d_t(j,i)
4454 if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
4456 d_t_work(ind+j)=d_t(j,i+nres)
4462 call fricmat_mult(d_t_work,fric_work)
4464 if (.not.checkmode) return
4467 write (iout,*) "d_t_work and fric_work"
4469 write (iout,'(i3,2e15.5)') i,d_t_work(i),fric_work(i)
4473 friction(j,0)=fric_work(j)
4478 friction(j,i)=fric_work(ind+j)
4483 if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
4485 friction(j,i+nres)=fric_work(ind+j)
4491 write(iout,*) "Friction backbone"
4493 write(iout,'(i5,3e15.5,5x,3e15.5)') &
4494 i,(friction(j,i),j=1,3),(d_t(j,i),j=1,3)
4496 write(iout,*) "Friction side chain"
4498 write(iout,'(i5,3e15.5,5x,3e15.5)') &
4499 i,(friction(j,i+nres),j=1,3),(d_t(j,i+nres),j=1,3)
4508 vvtot(j,i)=vv(j)+0.5d0*d_t(j,i)
4509 vvtot(j,i+nres)=vv(j)+d_t(j,i+nres)
4510 vv(j)=vv(j)+d_t(j,i)
4513 write (iout,*) "vvtot backbone and sidechain"
4515 write (iout,'(i5,3e15.5,5x,3e15.5)') i,(vvtot(j,i),j=1,3),&
4516 (vvtot(j,i+nres),j=1,3)
4521 v_work(ind+j)=vvtot(j,i)
4527 v_work(ind+j)=vvtot(j,i+nres)
4531 write (iout,*) "v_work gamvec and site-based friction forces"
4533 write (iout,'(i5,3e15.5)') i,v_work(i),gamvec(i),&
4537 ! fric_work1(i)=0.0d0
4539 ! fric_work1(i)=fric_work1(i)-A(j,i)*gamvec(j)*v_work(j)
4542 ! write (iout,*) "fric_work and fric_work1"
4544 ! write (iout,'(i5,2e15.5)') i,fric_work(i),fric_work1(i)
4550 ginvfric(i,j)=ginvfric(i,j)+ginv(i,k)*fricmat(k,j)
4554 write (iout,*) "ginvfric"
4556 write (iout,'(i5,100f8.3)') i,(ginvfric(i,j),j=1,dimen)
4558 write (iout,*) "symmetry check"
4561 write (iout,*) i,j,ginvfric(i,j)-ginvfric(j,i)
4566 end subroutine friction_force
4567 !-----------------------------------------------------------------------------
4568 subroutine setup_fricmat
4572 use control_data, only:time_Bcast
4573 use control, only:tcpu
4575 ! implicit real*8 (a-h,o-z)
4579 real(kind=8) :: time00
4581 ! include 'DIMENSIONS'
4582 ! include 'COMMON.VAR'
4583 ! include 'COMMON.CHAIN'
4584 ! include 'COMMON.DERIV'
4585 ! include 'COMMON.GEO'
4586 ! include 'COMMON.LOCAL'
4587 ! include 'COMMON.INTERACT'
4588 ! include 'COMMON.MD'
4589 ! include 'COMMON.SETUP'
4590 ! include 'COMMON.TIME1'
4591 ! integer licznik /0/
4594 ! include 'COMMON.LANGEVIN'
4596 ! include 'COMMON.LANGEVIN.lang0'
4598 ! include 'COMMON.IOUNITS'
4600 integer :: i,j,ind,ind1,m
4601 logical :: lprn = .false.
4602 real(kind=8) :: dtdi !el ,gamvec(2*nres)
4603 !el real(kind=8),dimension(2*nres,2*nres) :: ginvfric,fcopy
4604 real(kind=8),dimension(2*nres,2*nres) :: fcopy
4605 !el real(kind=8),dimension(2*nres*(2*nres+1)/2) :: Ghalf !(mmaxres2) (mmaxres2=(maxres2*(maxres2+1)/2))
4606 !el common /syfek/ gamvec
4607 real(kind=8) :: work(8*2*nres)
4608 integer :: iwork(2*nres)
4609 !el common /przechowalnia/ ginvfric,Ghalf,fcopy
4610 integer :: ii,iti,k,l,nzero,nres2,nres6,ierr
4612 if (fg_rank.ne.king) goto 10
4617 if(.not.allocated(gamvec)) allocate(gamvec(nres2)) !(MAXRES2)
4618 if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
4619 !el if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
4620 !el allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
4621 if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !maxres2=2*maxres
4623 !el if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
4624 ! Zeroing out fricmat
4630 ! Load the friction coefficients corresponding to peptide groups
4636 ! Load the friction coefficients corresponding to side chains
4644 gamvec(ii)=gamsc(iabs(iti))
4646 if (surfarea) call sdarea(gamvec)
4648 ! write (iout,*) "Matrix A and vector gamma"
4650 ! write (iout,'(i2,$)') i
4652 ! write (iout,'(f4.1,$)') A(i,j)
4654 ! write (iout,'(f8.3)') gamvec(i)
4658 write (iout,*) "Vector gamvec"
4660 write (iout,'(i5,f10.5)') i, gamvec(i)
4664 ! The friction matrix
4669 dtdi=dtdi+A(j,k)*A(j,i)*gamvec(j)
4676 write (iout,'(//a)') "Matrix fricmat"
4677 call matout2(dimen,dimen,nres2,nres2,fricmat)
4679 if (lang.eq.2 .or. lang.eq.3) then
4680 ! Mass-scale the friction matrix if non-direct integration will be performed
4686 Ginvfric(i,j)=Ginvfric(i,j)+ &
4687 Gsqrm(i,k)*Gsqrm(l,j)*fricmat(k,l)
4692 ! Diagonalize the friction matrix
4697 Ghalf(ind)=Ginvfric(i,j)
4700 call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
4703 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
4704 " mass-scaled friction matrix"
4705 call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
4707 ! Precompute matrices for tinker stochastic integrator
4714 mt1(i,j)=mt1(i,j)+fricvec(k,i)*gsqrm(k,j)
4715 mt2(i,j)=mt2(i,j)+fricvec(k,i)*gsqrp(k,j)
4721 else if (lang.eq.4) then
4722 ! Diagonalize the friction matrix
4727 Ghalf(ind)=fricmat(i,j)
4730 call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
4733 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
4735 call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
4737 ! Determine the number of zero eigenvalues of the friction matrix
4738 nzero=max0(dimen-dimen1,0)
4739 ! do while (fricgam(nzero+1).le.1.0d-5 .and. nzero.lt.dimen)
4742 write (iout,*) "Number of zero eigenvalues:",nzero
4747 fricmat(i,j)=fricmat(i,j) &
4748 +fricvec(i,k)*fricvec(j,k)/fricgam(k)
4753 write (iout,'(//a)') "Generalized inverse of fricmat"
4754 call matout(dimen,dimen,nres6,nres6,fricmat)
4759 if (nfgtasks.gt.1) then
4760 if (fg_rank.eq.0) then
4761 ! The matching BROADCAST for fg processors is called in ERGASTULUM
4767 call MPI_Bcast(10,1,MPI_INTEGER,king,FG_COMM,IERROR)
4769 time_Bcast=time_Bcast+MPI_Wtime()-time00
4771 time_Bcast=time_Bcast+tcpu()-time00
4773 ! print *,"Processor",myrank,
4774 ! & " BROADCAST iorder in SETUP_FRICMAT"
4777 write (iout,*) "setup_fricmat licznik"!,licznik !sp
4783 ! Scatter the friction matrix
4784 call MPI_Scatterv(fricmat(1,1),nginv_counts(0),&
4785 nginv_start(0),MPI_DOUBLE_PRECISION,fcopy(1,1),&
4786 myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
4789 time_scatter=time_scatter+MPI_Wtime()-time00
4790 time_scatter_fmat=time_scatter_fmat+MPI_Wtime()-time00
4792 time_scatter=time_scatter+tcpu()-time00
4793 time_scatter_fmat=time_scatter_fmat+tcpu()-time00
4797 do j=1,2*my_ng_count
4798 fricmat(j,i)=fcopy(i,j)
4801 ! write (iout,*) "My chunk of fricmat"
4802 ! call MATOUT2(my_ng_count,dimen,maxres2,maxres2,fcopy)
4806 end subroutine setup_fricmat
4807 !-----------------------------------------------------------------------------
4808 subroutine stochastic_force(stochforcvec)
4811 use random, only:anorm_distr
4812 ! implicit real*8 (a-h,o-z)
4813 ! include 'DIMENSIONS'
4814 use control, only: tcpu
4819 ! include 'COMMON.VAR'
4820 ! include 'COMMON.CHAIN'
4821 ! include 'COMMON.DERIV'
4822 ! include 'COMMON.GEO'
4823 ! include 'COMMON.LOCAL'
4824 ! include 'COMMON.INTERACT'
4825 ! include 'COMMON.MD'
4826 ! include 'COMMON.TIME1'
4828 ! include 'COMMON.LANGEVIN'
4830 ! include 'COMMON.LANGEVIN.lang0'
4832 ! include 'COMMON.IOUNITS'
4834 real(kind=8) :: x,sig,lowb,highb
4835 real(kind=8) :: ff(3),force(3,0:2*nres),zeta2,lowb2
4836 real(kind=8) :: highb2,sig2,forcvec(6*nres),stochforcvec(6*nres)
4837 real(kind=8) :: time00
4838 logical :: lprn = .false.
4843 stochforc(j,i)=0.0d0
4853 ! Compute the stochastic forces acting on bodies. Store in force.
4859 force(j,i)=anorm_distr(x,sig,lowb,highb)
4867 force(j,i+nres)=anorm_distr(x,sig2,lowb2,highb2)
4871 time_fsample=time_fsample+MPI_Wtime()-time00
4873 time_fsample=time_fsample+tcpu()-time00
4875 ! Compute the stochastic forces acting on virtual-bond vectors.
4881 stochforc(j,i)=ff(j)+0.5d0*force(j,i)
4884 ff(j)=ff(j)+force(j,i)
4886 if (itype(i+1).ne.ntyp1) then
4888 stochforc(j,i)=stochforc(j,i)+force(j,i+nres+1)
4889 ff(j)=ff(j)+force(j,i+nres+1)
4894 stochforc(j,0)=ff(j)+force(j,nnt+nres)
4897 if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
4899 stochforc(j,i+nres)=force(j,i+nres)
4905 stochforcvec(j)=stochforc(j,0)
4910 stochforcvec(ind+j)=stochforc(j,i)
4915 if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
4917 stochforcvec(ind+j)=stochforc(j,i+nres)
4923 write (iout,*) "stochforcvec"
4925 write(iout,'(i5,e15.5)') i,stochforcvec(i)
4927 write(iout,*) "Stochastic forces backbone"
4929 write(iout,'(i5,3e15.5)') i,(stochforc(j,i),j=1,3)
4931 write(iout,*) "Stochastic forces side chain"
4933 write(iout,'(i5,3e15.5)') &
4934 i,(stochforc(j,i+nres),j=1,3)
4942 write (iout,*) i,ind
4944 forcvec(ind+j)=force(j,i)
4949 write (iout,*) i,ind
4951 forcvec(j+ind)=force(j,i+nres)
4956 write (iout,*) "forcvec"
4960 write (iout,'(2i3,2f10.5)') i,j,force(j,i),&
4967 write (iout,'(2i3,2f10.5)') i,j,force(j,i+nres),&
4976 end subroutine stochastic_force
4977 !-----------------------------------------------------------------------------
4978 subroutine sdarea(gamvec)
4980 ! Scale the friction coefficients according to solvent accessible surface areas
4981 ! Code adapted from TINKER
4985 ! implicit real*8 (a-h,o-z)
4986 ! include 'DIMENSIONS'
4987 ! include 'COMMON.CONTROL'
4988 ! include 'COMMON.VAR'
4989 ! include 'COMMON.MD'
4991 ! include 'COMMON.LANGEVIN'
4993 ! include 'COMMON.LANGEVIN.lang0'
4995 ! include 'COMMON.CHAIN'
4996 ! include 'COMMON.DERIV'
4997 ! include 'COMMON.GEO'
4998 ! include 'COMMON.LOCAL'
4999 ! include 'COMMON.INTERACT'
5000 ! include 'COMMON.IOUNITS'
5001 ! include 'COMMON.NAMES'
5002 real(kind=8),dimension(2*nres) :: radius,gamvec !(maxres2)
5003 real(kind=8),parameter :: twosix = 1.122462048309372981d0
5004 logical :: lprn = .false.
5005 real(kind=8) :: probe,area,ratio
5006 integer :: i,j,ind,iti
5008 ! determine new friction coefficients every few SD steps
5010 ! set the atomic radii to estimates of sigma values
5012 ! print *,"Entered sdarea"
5018 ! Load peptide group radii
5022 ! Load side chain radii
5025 radius(i+nres)=restok(iti)
5028 ! write (iout,*) "i",i," radius",radius(i)
5031 radius(i) = radius(i) / twosix
5032 if (radius(i) .ne. 0.0d0) radius(i) = radius(i) + probe
5035 ! scale atomic friction coefficients by accessible area
5037 if (lprn) write (iout,*) &
5038 "Original gammas, surface areas, scaling factors, new gammas, ",&
5039 "std's of stochastic forces"
5042 if (radius(i).gt.0.0d0) then
5043 call surfatom (i,area,radius)
5044 ratio = dmax1(area/(4.0d0*pi*radius(i)**2),1.0d-1)
5045 if (lprn) write (iout,'(i5,3f10.5,$)') &
5046 i,gamvec(ind+1),area,ratio
5049 gamvec(ind) = ratio * gamvec(ind)
5051 stdforcp(i)=stdfp*dsqrt(gamvec(ind))
5052 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcp(i)
5056 if (radius(i+nres).gt.0.0d0) then
5057 call surfatom (i+nres,area,radius)
5058 ratio = dmax1(area/(4.0d0*pi*radius(i+nres)**2),1.0d-1)
5059 if (lprn) write (iout,'(i5,3f10.5,$)') &
5060 i,gamvec(ind+1),area,ratio
5063 gamvec(ind) = ratio * gamvec(ind)
5065 stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamvec(ind))
5066 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcsc(i)
5071 end subroutine sdarea
5072 !-----------------------------------------------------------------------------
5074 !-----------------------------------------------------------------------------
5077 ! ###################################################
5078 ! ## COPYRIGHT (C) 1996 by Jay William Ponder ##
5079 ! ## All Rights Reserved ##
5080 ! ###################################################
5082 ! ################################################################
5084 ! ## subroutine surfatom -- exposed surface area of an atom ##
5086 ! ################################################################
5089 ! "surfatom" performs an analytical computation of the surface
5090 ! area of a specified atom; a simplified version of "surface"
5092 ! literature references:
5094 ! T. J. Richmond, "Solvent Accessible Surface Area and
5095 ! Excluded Volume in Proteins", Journal of Molecular Biology,
5098 ! L. Wesson and D. Eisenberg, "Atomic Solvation Parameters
5099 ! Applied to Molecular Dynamics of Proteins in Solution",
5100 ! Protein Science, 1, 227-235 (1992)
5102 ! variables and parameters:
5104 ! ir number of atom for which area is desired
5105 ! area accessible surface area of the atom
5106 ! radius radii of each of the individual atoms
5109 subroutine surfatom(ir,area,radius)
5111 ! implicit real*8 (a-h,o-z)
5112 ! include 'DIMENSIONS'
5114 ! include 'COMMON.GEO'
5115 ! include 'COMMON.IOUNITS'
5117 integer :: nsup,nstart_sup
5118 ! double precision c,dc,dc_old,d_c_work,xloc,xrot,dc_norm
5119 ! common /chain/ c(3,maxres2+2),dc(3,0:maxres2),dc_old(3,0:maxres2),
5120 ! & xloc(3,maxres),xrot(3,maxres),dc_norm(3,0:maxres2),
5121 ! & dc_work(MAXRES6),nres,nres0
5122 integer,parameter :: maxarc=300
5126 integer :: mi,ni,narc
5127 integer :: key(maxarc)
5128 integer :: intag(maxarc)
5129 integer :: intag1(maxarc)
5130 real(kind=8) :: area,arcsum
5131 real(kind=8) :: arclen,exang
5132 real(kind=8) :: delta,delta2
5133 real(kind=8) :: eps,rmove
5134 real(kind=8) :: xr,yr,zr
5135 real(kind=8) :: rr,rrsq
5136 real(kind=8) :: rplus,rminus
5137 real(kind=8) :: axx,axy,axz
5138 real(kind=8) :: ayx,ayy
5139 real(kind=8) :: azx,azy,azz
5140 real(kind=8) :: uxj,uyj,uzj
5141 real(kind=8) :: tx,ty,tz
5142 real(kind=8) :: txb,tyb,td
5143 real(kind=8) :: tr2,tr,txr,tyr
5144 real(kind=8) :: tk1,tk2
5145 real(kind=8) :: thec,the,t,tb
5146 real(kind=8) :: txk,tyk,tzk
5147 real(kind=8) :: t1,ti,tf,tt
5148 real(kind=8) :: txj,tyj,tzj
5149 real(kind=8) :: ccsq,cc,xysq
5150 real(kind=8) :: bsqk,bk,cosine
5151 real(kind=8) :: dsqj,gi,pix2
5152 real(kind=8) :: therk,dk,gk
5153 real(kind=8) :: risqk,rik
5154 real(kind=8) :: radius(2*nres) !(maxatm) (maxatm=maxres2)
5155 real(kind=8) :: ri(maxarc),risq(maxarc)
5156 real(kind=8) :: ux(maxarc),uy(maxarc),uz(maxarc)
5157 real(kind=8) :: xc(maxarc),yc(maxarc),zc(maxarc)
5158 real(kind=8) :: xc1(maxarc),yc1(maxarc),zc1(maxarc)
5159 real(kind=8) :: dsq(maxarc),bsq(maxarc)
5160 real(kind=8) :: dsq1(maxarc),bsq1(maxarc)
5161 real(kind=8) :: arci(maxarc),arcf(maxarc)
5162 real(kind=8) :: ex(maxarc),lt(maxarc),gr(maxarc)
5163 real(kind=8) :: b(maxarc),b1(maxarc),bg(maxarc)
5164 real(kind=8) :: kent(maxarc),kout(maxarc)
5165 real(kind=8) :: ther(maxarc)
5166 logical :: moved,top
5167 logical :: omit(maxarc)
5170 maxatm = 2*nres !maxres2 maxres2=2*maxres
5176 ! zero out the surface area for the sphere of interest
5179 ! write (2,*) "ir",ir," radius",radius(ir)
5180 if (radius(ir) .eq. 0.0d0) return
5182 ! set the overlap significance and connectivity shift
5186 delta2 = delta * delta
5191 ! store coordinates and radius of the sphere of interest
5199 ! initialize values of some counters and summations
5208 ! test each sphere to see if it overlaps the sphere of interest
5211 if (i.eq.ir .or. radius(i).eq.0.0d0) goto 30
5212 rplus = rr + radius(i)
5214 if (abs(tx) .ge. rplus) goto 30
5216 if (abs(ty) .ge. rplus) goto 30
5218 if (abs(tz) .ge. rplus) goto 30
5220 ! check for sphere overlap by testing distance against radii
5222 xysq = tx*tx + ty*ty
5223 if (xysq .lt. delta2) then
5230 if (rplus-cc .le. delta) goto 30
5231 rminus = rr - radius(i)
5233 ! check to see if sphere of interest is completely buried
5235 if (cc-abs(rminus) .le. delta) then
5236 if (rminus .le. 0.0d0) goto 170
5240 ! check for too many overlaps with sphere of interest
5242 if (io .ge. maxarc) then
5244 20 format (/,' SURFATOM -- Increase the Value of MAXARC')
5248 ! get overlap between current sphere and sphere of interest
5257 gr(io) = (ccsq+rplus*rminus) / (2.0d0*rr*b1(io))
5263 ! case where no other spheres overlap the sphere of interest
5266 area = 4.0d0 * pi * rrsq
5270 ! case where only one sphere overlaps the sphere of interest
5273 area = pix2 * (1.0d0 + gr(1))
5274 area = mod(area,4.0d0*pi) * rrsq
5278 ! case where many spheres intersect the sphere of interest;
5279 ! sort the intersecting spheres by their degree of overlap
5281 call sort2 (io,gr,key)
5284 intag(i) = intag1(k)
5293 ! get radius of each overlap circle on surface of the sphere
5298 risq(i) = rrsq - gi*gi
5299 ri(i) = sqrt(risq(i))
5300 ther(i) = 0.5d0*pi - asin(min(1.0d0,max(-1.0d0,gr(i))))
5303 ! find boundary of inaccessible area on sphere of interest
5306 if (.not. omit(k)) then
5313 ! check to see if J circle is intersecting K circle;
5314 ! get distance between circle centers and sum of radii
5317 if (omit(j)) goto 60
5318 cc = (txk*xc(j)+tyk*yc(j)+tzk*zc(j))/(bk*b(j))
5319 cc = acos(min(1.0d0,max(-1.0d0,cc)))
5320 td = therk + ther(j)
5322 ! check to see if circles enclose separate regions
5324 if (cc .ge. td) goto 60
5326 ! check for circle J completely inside circle K
5328 if (cc+ther(j) .lt. therk) goto 40
5330 ! check for circles that are essentially parallel
5332 if (cc .gt. delta) goto 50
5337 ! check to see if sphere of interest is completely buried
5340 if (pix2-cc .le. td) goto 170
5346 ! find T value of circle intersections
5349 if (omit(k)) goto 110
5364 ! rotation matrix elements
5376 if (.not. omit(j)) then
5381 ! rotate spheres so K vector colinear with z-axis
5383 uxj = txj*axx + tyj*axy - tzj*axz
5384 uyj = tyj*ayy - txj*ayx
5385 uzj = txj*azx + tyj*azy + tzj*azz
5386 cosine = min(1.0d0,max(-1.0d0,uzj/b(j)))
5387 if (acos(cosine) .lt. therk+ther(j)) then
5388 dsqj = uxj*uxj + uyj*uyj
5393 tr2 = risqk*dsqj - tb*tb
5399 ! get T values of intersection for K circle
5402 tb = min(1.0d0,max(-1.0d0,tb))
5404 if (tyb-txr .lt. 0.0d0) tk1 = pix2 - tk1
5406 tb = min(1.0d0,max(-1.0d0,tb))
5408 if (tyb+txr .lt. 0.0d0) tk2 = pix2 - tk2
5409 thec = (rrsq*uzj-gk*bg(j)) / (rik*ri(j)*b(j))
5410 if (abs(thec) .lt. 1.0d0) then
5412 else if (thec .ge. 1.0d0) then
5414 else if (thec .le. -1.0d0) then
5418 ! see if "tk1" is entry or exit point; check t=0 point;
5419 ! "ti" is exit point, "tf" is entry point
5421 cosine = min(1.0d0,max(-1.0d0, &
5422 (uzj*gk-uxj*rik)/(b(j)*rr)))
5423 if ((acos(cosine)-ther(j))*(tk2-tk1) .le. 0.0d0) then
5431 if (narc .ge. maxarc) then
5433 70 format (/,' SURFATOM -- Increase the Value',&
5437 if (tf .le. ti) then
5458 ! special case; K circle without intersections
5460 if (narc .le. 0) goto 90
5462 ! general case; sum up arclength and set connectivity code
5464 call sort2 (narc,arci,key)
5469 if (narc .gt. 1) then
5472 if (t .lt. arci(j)) then
5473 arcsum = arcsum + arci(j) - t
5474 exang = exang + ex(ni)
5476 if (jb .ge. maxarc) then
5478 80 format (/,' SURFATOM -- Increase the Value',&
5483 kent(jb) = maxarc*i + k
5485 kout(jb) = maxarc*k + i
5494 arcsum = arcsum + pix2 - t
5496 exang = exang + ex(ni)
5499 kent(jb) = maxarc*i + k
5501 kout(jb) = maxarc*k + i
5508 arclen = arclen + gr(k)*arcsum
5511 if (arclen .eq. 0.0d0) goto 170
5512 if (jb .eq. 0) goto 150
5514 ! find number of independent boundaries and check connectivity
5518 if (kout(k) .ne. 0) then
5525 if (m .eq. kent(ii)) then
5528 if (j .eq. jb) goto 150
5540 ! attempt to fix connectivity error by moving atom slightly
5544 140 format (/,' SURFATOM -- Connectivity Error at Atom',i6)
5553 ! compute the exposed surface area for the sphere of interest
5556 area = ib*pix2 + exang + arclen
5557 area = mod(area,4.0d0*pi) * rrsq
5559 ! attempt to fix negative area by moving atom slightly
5561 if (area .lt. 0.0d0) then
5564 160 format (/,' SURFATOM -- Negative Area at Atom',i6)
5575 end subroutine surfatom
5576 !----------------------------------------------------------------
5577 !----------------------------------------------------------------
5578 subroutine alloc_MD_arrays
5579 !EL Allocation of arrays used by MD module
5581 integer :: nres2,nres6
5584 !----------------------
5588 allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
5589 allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
5590 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
5591 allocate(fricvec(nres2,nres2))
5592 allocate(pfric_mat(nres2,nres2),vfric_mat(nres2,nres2))
5593 allocate(afric_mat(nres2,nres2),prand_mat(nres2,nres2))
5594 allocate(vrand_mat1(nres2,nres2),vrand_mat2(nres2,nres2)) !(MAXRES2,MAXRES2)
5595 allocate(pfric0_mat(nres2,nres2,0:maxflag_stoch))
5596 allocate(afric0_mat(nres2,nres2,0:maxflag_stoch))
5597 allocate(vfric0_mat(nres2,nres2,0:maxflag_stoch))
5598 allocate(prand0_mat(nres2,nres2,0:maxflag_stoch))
5599 allocate(vrand0_mat1(nres2,nres2,0:maxflag_stoch))
5600 allocate(vrand0_mat2(nres2,nres2,0:maxflag_stoch)) !(MAXRES2,MAXRES2,0:maxflag_stoch)
5601 allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
5603 allocate(mt1(nres2,nres2),mt2(nres2,nres2),mt3(nres2,nres2)) !(maxres2,maxres2)
5604 !----------------------
5606 ! commom.langevin.lang0
5608 allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
5609 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
5610 allocate(fricvec(nres2,nres2)) !(MAXRES2,MAXRES2)
5611 allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
5612 allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
5615 !el if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2))
5616 !----------------------
5617 ! commom.hairpin in CSA module
5618 !----------------------
5619 ! common.mce in MCM_MD module
5620 !----------------------
5622 ! common /mdgrad/ in module.energy
5623 ! common /back_constr/ in module.energy
5624 ! common /qmeas/ in module.energy
5627 allocate(potEcomp(0:n_ene+4)) !(0:n_ene+4)
5629 allocate(d_t(3,0:nres2),d_a(3,0:nres2),d_t_old(3,0:nres2)) !(3,0:MAXRES2)
5630 allocate(d_a_work(nres6)) !(6*MAXRES)
5631 allocate(Gmat(nres2,nres2),A(nres2,nres2))
5632 if(.not.allocated(Ginv)) allocate(Ginv(nres2,nres2)) !in control: ergastulum
5633 allocate(Gsqrp(nres2,nres2),Gsqrm(nres2,nres2),Gvec(nres2,nres2)) !(maxres2,maxres2)
5634 allocate(Geigen(nres2)) !(maxres2)
5635 if(.not.allocated(vtot)) allocate(vtot(nres2)) !(maxres2)
5636 ! common /inertia/ in io_conf: parmread
5637 ! real(kind=8),dimension(:),allocatable :: ISC,msc !(ntyp+1)
5638 ! common /langevin/in io read_MDpar
5639 ! real(kind=8),dimension(:),allocatable :: gamsc !(ntyp1)
5640 ! real(kind=8),dimension(:),allocatable :: stdfsc !(ntyp)
5641 ! in io_conf: parmread
5642 ! real(kind=8),dimension(:),allocatable :: restok !(ntyp+1)
5643 ! common /mdpmpi/ in control: ergastulum
5644 if(.not.allocated(ng_start)) allocate(ng_start(0:nfgtasks-1))
5645 if(.not.allocated(ng_counts)) allocate(ng_counts(0:nfgtasks-1))
5646 if(.not.allocated(nginv_counts)) allocate(nginv_counts(0:nfgtasks-1)) !(0:MaxProcs-1)
5647 if(.not.allocated(nginv_start)) allocate(nginv_start(0:nfgtasks)) !(0:MaxProcs)
5648 !----------------------
5649 ! common.muca in read_muca
5650 ! common /double_muca/
5651 ! real(kind=8) :: elow,ehigh,factor,hbin,factor_min
5652 ! real(kind=8),dimension(:),allocatable :: emuca,nemuca,&
5653 ! nemuca2,hist !(4*maxres)
5654 ! common /integer_muca/
5655 ! integer :: nmuca,imtime,muca_smooth
5657 ! real(kind=8),dimension(:),allocatable :: elowi,ehighi !(maxprocs)
5658 !----------------------
5660 ! common /mdgrad/ in module.energy
5661 ! common /back_constr/ in module.energy
5662 ! common /qmeas/ in module.energy
5666 allocate(d_t_work(nres6),d_t_work_new(nres6),d_af_work(nres6))
5667 allocate(d_as_work(nres6),kinetic_force(nres6)) !(MAXRES6)
5668 allocate(d_t_new(3,0:nres2),d_a_old(3,0:nres2),d_a_short(3,0:nres2)) !,d_a !(3,0:MAXRES2)
5669 allocate(stdforcp(nres),stdforcsc(nres)) !(MAXRES)
5670 !----------------------
5672 allocate(D_ban(nres6)) !(MAXRES6) maxres6=6*maxres
5673 ! common /stochcalc/ stochforcvec
5674 allocate(stochforcvec(nres6)) !(MAXRES6) maxres6=6*maxres
5675 !----------------------
5677 end subroutine alloc_MD_arrays
5678 !-----------------------------------------------------------------------------
5679 !-----------------------------------------------------------------------------