2 !-----------------------------------------------------------------------------
15 !-----------------------------------------------------------------------------
17 ! common /mdgrad/ in module.energy
18 ! common /back_constr/ in module.energy
19 ! common /qmeas/ in module.energy
23 real(kind=8),dimension(:),allocatable :: d_t_work,&
24 d_t_work_new,d_af_work,d_as_work,kinetic_force !(MAXRES6)
25 real(kind=8),dimension(:,:),allocatable :: d_t_new,&
26 d_a_old,d_a_short!,d_a !(3,0:MAXRES2)
27 ! real(kind=8),dimension(:),allocatable :: d_a_work !(6*MAXRES)
28 ! real(kind=8),dimension(:,:),allocatable :: Gmat,Ginv,A,&
29 ! Gsqrp,Gsqrm,Gvec !(maxres2,maxres2)
30 ! real(kind=8),dimension(:),allocatable :: Geigen !(maxres2)
31 ! integer :: dimen,dimen1,dimen3
32 ! integer :: lang,count_reset_moment,count_reset_vel
33 ! logical :: reset_moment,reset_vel,rattle,RESPA
36 ! real(kind=8) :: rwat,etawat,stdfp,cPoise
37 ! real(kind=8),dimension(:),allocatable :: gamsc !(ntyp1)
38 ! real(kind=8),dimension(:),allocatable :: stdfsc !(ntyp)
39 real(kind=8),dimension(:),allocatable :: stdforcp,stdforcsc !(MAXRES)
40 !-----------------------------------------------------------------------------
44 ! ###################################################
45 ! ## COPYRIGHT (C) 1992 by Jay William Ponder ##
46 ! ## All Rights Reserved ##
47 ! ###################################################
49 ! #############################################################
51 ! ## sizes.i -- parameter values to set array dimensions ##
53 ! #############################################################
56 ! "sizes.i" sets values for critical array dimensions used
57 ! throughout the software; these parameters will fix the size
58 ! of the largest systems that can be handled; values too large
59 ! for the computer's memory and/or swap space to accomodate
60 ! will result in poor performance or outright failure
62 ! parameter: maximum allowed number of:
64 ! maxatm atoms in the molecular system
65 ! maxval atoms directly bonded to an atom
66 ! maxgrp ! user-defined groups of atoms
67 ! maxtyp force field atom type definitions
68 ! maxclass force field atom class definitions
69 ! maxkey lines in the keyword file
70 ! maxrot bonds for torsional rotation
71 ! maxvar optimization variables (vector storage)
72 ! maxopt optimization variables (matrix storage)
73 ! maxhess off-diagonal Hessian elements
74 ! maxlight sites for method of lights neighbors
75 ! maxvib vibrational frequencies
76 ! maxgeo distance geometry points
77 ! maxcell unit cells in replicated crystal
78 ! maxring 3-, 4-, or 5-membered rings
79 ! maxfix geometric restraints
80 ! maxbio biopolymer atom definitions
81 ! maxres residues in the macromolecule
82 ! maxamino amino acid residue types
83 ! maxnuc nucleic acid residue types
84 ! maxbnd covalent bonds in molecular system
85 ! maxang bond angles in molecular system
86 ! maxtors torsional angles in molecular system
87 ! maxpi atoms in conjugated pisystem
88 ! maxpib covalent bonds involving pisystem
89 ! maxpit torsional angles involving pisystem
92 !el integer maxatm,maxval,maxgrp
93 !el integer maxtyp,maxclass,maxkey
94 !el integer maxrot,maxopt
95 !el integer maxhess,maxlight,maxvib
96 !el integer maxgeo,maxcell,maxring
97 !el integer maxfix,maxbio
98 !el integer maxamino,maxnuc,maxbnd
99 !el integer maxang,maxtors,maxpi
100 !el integer maxpib,maxpit
101 ! integer :: maxatm !=2*nres !maxres2 maxres2=2*maxres
102 ! integer,parameter :: maxval=8
103 ! integer,parameter :: maxgrp=1000
104 ! integer,parameter :: maxtyp=3000
105 ! integer,parameter :: maxclass=500
106 ! integer,parameter :: maxkey=10000
107 ! integer,parameter :: maxrot=1000
108 ! integer,parameter :: maxopt=1000
109 ! integer,parameter :: maxhess=1000000
110 ! integer :: maxlight !=8*maxatm
111 ! integer,parameter :: maxvib=1000
112 ! integer,parameter :: maxgeo=1000
113 ! integer,parameter :: maxcell=10000
114 ! integer,parameter :: maxring=10000
115 ! integer,parameter :: maxfix=10000
116 ! integer,parameter :: maxbio=10000
117 ! integer,parameter :: maxamino=31
118 ! integer,parameter :: maxnuc=12
119 ! integer :: maxbnd !=2*maxatm
120 ! integer :: maxang !=3*maxatm
121 ! integer :: maxtors !=4*maxatm
122 ! integer,parameter :: maxpi=100
123 ! integer,parameter :: maxpib=2*maxpi
124 ! integer,parameter :: maxpit=4*maxpi
125 !-----------------------------------------------------------------------------
126 ! Maximum number of seed
127 ! integer,parameter :: max_seed=1
128 !-----------------------------------------------------------------------------
129 real(kind=8),dimension(:),allocatable :: stochforcvec !(MAXRES6) maxres6=6*maxres
130 ! common /stochcalc/ stochforcvec
131 !-----------------------------------------------------------------------------
132 ! common /przechowalnia/ subroutines: rattle1,rattle2,rattle_brown
133 real(kind=8),dimension(:,:),allocatable :: GGinv !(2*nres,2*nres) maxres2=2*maxres
134 real(kind=8),dimension(:,:,:),allocatable :: gdc !(3,2*nres,2*nres) maxres2=2*maxres
135 real(kind=8),dimension(:,:),allocatable :: Cmat !(2*nres,2*nres) maxres2=2*maxres
136 !-----------------------------------------------------------------------------
137 ! common /syfek/ subroutines: friction_force,setup_fricmat
138 !el real(kind=8),dimension(:),allocatable :: gamvec !(MAXRES6) or (MAXRES2)
139 !-----------------------------------------------------------------------------
140 ! common /przechowalnia/ subroutines: friction_force,setup_fricmat
141 real(kind=8),dimension(:,:),allocatable :: ginvfric !(2*nres,2*nres) !maxres2=2*maxres
142 !-----------------------------------------------------------------------------
143 ! common /przechowalnia/ subroutine: setup_fricmat
144 real(kind=8),dimension(:,:),allocatable :: fcopy !(2*nres,2*nres)
145 !-----------------------------------------------------------------------------
148 !-----------------------------------------------------------------------------
150 !-----------------------------------------------------------------------------
152 !-----------------------------------------------------------------------------
153 subroutine brown_step(itime)
154 !------------------------------------------------
155 ! Perform a single Euler integration step of Brownian dynamics
156 !------------------------------------------------
157 ! implicit real*8 (a-h,o-z)
159 use control, only: tcpu
162 ! use io_conf, only:cartprint
163 ! include 'DIMENSIONS'
167 ! include 'COMMON.CONTROL'
168 ! include 'COMMON.VAR'
169 ! include 'COMMON.MD'
171 ! include 'COMMON.LANGEVIN'
173 ! include 'COMMON.LANGEVIN.lang0'
175 ! include 'COMMON.CHAIN'
176 ! include 'COMMON.DERIV'
177 ! include 'COMMON.GEO'
178 ! include 'COMMON.LOCAL'
179 ! include 'COMMON.INTERACT'
180 ! include 'COMMON.IOUNITS'
181 ! include 'COMMON.NAMES'
182 ! include 'COMMON.TIME1'
183 real(kind=8),dimension(6*nres) :: zapas !(MAXRES6) maxres6=6*maxres
184 integer :: rstcount !ilen,
186 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
187 ! real(kind=8),dimension(6*nres,2*nres) :: Bmat,GBmat,Tmat !(MAXRES6,MAXRES2) (maxres2=2*maxres,maxres6=6*maxres)
188 ! real(kind=8),dimension(2*nres,2*nres) :: Cmat_,Cinv !(maxres2,maxres2) maxres2=2*maxres
189 ! real(kind=8),dimension(6*nres,6*nres) :: Pmat !(maxres6,maxres6) maxres6=6*maxres
190 real(kind=8),dimension(:,:),allocatable :: Bmat,GBmat,Tmat !(MAXRES6,MAXRES2) (maxres2=2*maxres,maxres6=6*maxres)
191 real(kind=8),dimension(:,:),allocatable :: Cmat_,Cinv !(maxres2,maxres2) maxres2=2*maxres
192 real(kind=8),dimension(:,:),allocatable :: Pmat !(maxres6,maxres6) maxres6=6*maxres
193 real(kind=8),dimension(6*nres) :: Td !(maxres6) maxres6=6*maxres
194 real(kind=8),dimension(2*nres) :: ppvec !(maxres2) maxres2=2*maxres
195 !el common /stochcalc/ stochforcvec
196 !el real(kind=8),dimension(3) :: cm !el
197 !el common /gucio/ cm
199 logical :: lprn = .false.,lprn1 = .false.
200 integer :: maxiter = 5
201 real(kind=8) :: difftol = 1.0d-5
202 real(kind=8) :: xx,diffmax,blen2,diffbond,tt0
203 integer :: i,j,nbond,k,ind,ind1,iter
204 integer :: nres2,nres6
208 allocate(Bmat(6*nres,2*nres),GBmat(6*nres,2*nres),Tmat(6*nres,2*nres))
209 allocate(Cmat_(2*nres,2*nres),Cinv(2*nres,2*nres))
210 allocate(Pmat(6*nres,6*nres))
212 if (.not.allocated(stochforcvec)) allocate(stochforcvec(nres6)) !(MAXRES6) maxres6=6*maxres
216 if (itype(i,1).ne.10) nbond=nbond+1
220 write (iout,*) "Generalized inverse of fricmat"
221 call matout(dimen,dimen,nres6,nres6,fricmat)
233 Bmat(ind+j,ind1)=dC_norm(j,i)
238 if (itype(i,1).ne.10) then
241 Bmat(ind+j,ind1)=dC_norm(j,i+nres)
247 write (iout,*) "Matrix Bmat"
248 call MATOUT(nbond,dimen,nres6,nres6,Bmat)
254 GBmat(i,j)=GBmat(i,j)+fricmat(i,k)*Bmat(k,j)
259 write (iout,*) "Matrix GBmat"
260 call MATOUT(nbond,dimen,nres6,nres2,Gbmat)
266 Cmat_(i,j)=Cmat_(i,j)+Bmat(k,i)*GBmat(k,j)
271 write (iout,*) "Matrix Cmat"
272 call MATOUT(nbond,nbond,nres2,nres2,Cmat_)
274 call matinvert(nbond,nres2,Cmat_,Cinv,osob)
276 write (iout,*) "Matrix Cinv"
277 call MATOUT(nbond,nbond,nres2,nres2,Cinv)
283 Tmat(i,j)=Tmat(i,j)+GBmat(i,k)*Cinv(k,j)
288 write (iout,*) "Matrix Tmat"
289 call MATOUT(nbond,dimen,nres6,nres2,Tmat)
299 Pmat(i,j)=Pmat(i,j)-Tmat(i,k)*Bmat(j,k)
304 write (iout,*) "Matrix Pmat"
305 call MATOUT(dimen,dimen,nres6,nres6,Pmat)
312 Td(i)=Td(i)+vbl*Tmat(i,ind)
315 if (itype(k,1).ne.10) then
317 Td(i)=Td(i)+vbldsc0(1,itype(k,1))*Tmat(i,ind)
322 write (iout,*) "Vector Td"
324 write (iout,'(i5,f10.5)') i,Td(i)
327 call stochastic_force(stochforcvec)
329 write (iout,*) "stochforcvec"
331 write (iout,*) i,stochforcvec(i)
335 zapas(j)=-gcart(j,0)+stochforcvec(j)
337 dC_work(j)=dC_old(j,0)
343 zapas(ind)=-gcart(j,i)+stochforcvec(ind)
344 dC_work(ind)=dC_old(j,i)
348 if (itype(i,1).ne.10) then
351 zapas(ind)=-gxcart(j,i)+stochforcvec(ind)
352 dC_work(ind)=dC_old(j,i+nres)
358 write (iout,*) "Initial d_t_work"
360 write (iout,*) i,d_t_work(i)
367 d_t_work(i)=d_t_work(i)+fricmat(i,j)*zapas(j)
374 zapas(i)=zapas(i)+Pmat(i,j)*(dC_work(j)+d_t_work(j)*d_time)
378 write (iout,*) "Final d_t_work and zapas"
380 write (iout,*) i,d_t_work(i),zapas(i)
394 dc_work(ind+j)=dc(j,i)
400 d_t(j,i+nres)=d_t_work(ind+j)
401 dc(j,i+nres)=zapas(ind+j)
402 dc_work(ind+j)=dc(j,i+nres)
408 write (iout,*) "Before correction for rotational lengthening"
409 write (iout,*) "New coordinates",&
410 " and differences between actual and standard bond lengths"
415 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
419 if (itype(i,1).ne.10) then
421 xx=vbld(i+nres)-vbldsc0(1,itype(i,1))
422 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
423 i,(dC(j,i+nres),j=1,3),xx
427 ! Second correction (rotational lengthening)
433 blen2 = scalar(dc(1,i),dc(1,i))
434 ppvec(ind)=2*vbl**2-blen2
435 diffbond=dabs(vbl-dsqrt(blen2))
436 if (diffbond.gt.diffmax) diffmax=diffbond
437 if (ppvec(ind).gt.0.0d0) then
438 ppvec(ind)=dsqrt(ppvec(ind))
443 write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind)
447 if (itype(i,1).ne.10) then
449 blen2 = scalar(dc(1,i+nres),dc(1,i+nres))
450 ppvec(ind)=2*vbldsc0(1,itype(i,1))**2-blen2
451 diffbond=dabs(vbldsc0(1,itype(i,1))-dsqrt(blen2))
452 if (diffbond.gt.diffmax) diffmax=diffbond
453 if (ppvec(ind).gt.0.0d0) then
454 ppvec(ind)=dsqrt(ppvec(ind))
459 write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind)
463 if (lprn) write (iout,*) "iter",iter," diffmax",diffmax
464 if (diffmax.lt.difftol) goto 10
468 Td(i)=Td(i)+ppvec(j)*Tmat(i,j)
474 zapas(i)=zapas(i)+Pmat(i,j)*dc_work(j)
485 dc_work(ind+j)=zapas(ind+j)
490 if (itype(i,1).ne.10) then
492 dc(j,i+nres)=zapas(ind+j)
493 dc_work(ind+j)=zapas(ind+j)
498 ! Building the chain from the newly calculated coordinates
501 if (large.and. mod(itime,ntwe).eq.0) then
502 write (iout,*) "Cartesian and internal coordinates: step 1"
505 write (iout,'(a)') "Potential forces"
507 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(-gcart(j,i),j=1,3),&
510 write (iout,'(a)') "Stochastic forces"
512 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(stochforc(j,i),j=1,3),&
513 (stochforc(j,i+nres),j=1,3)
515 write (iout,'(a)') "Velocities"
517 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
518 (d_t(j,i+nres),j=1,3)
523 write (iout,*) "After correction for rotational lengthening"
524 write (iout,*) "New coordinates",&
525 " and differences between actual and standard bond lengths"
530 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
534 if (itype(i,1).ne.10) then
536 xx=vbld(i+nres)-vbldsc0(1,itype(i,1))
537 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
538 i,(dC(j,i+nres),j=1,3),xx
543 ! write (iout,*) "Too many attempts at correcting the bonds"
551 ! Calculate energy and forces
553 call etotal(potEcomp)
554 potE=potEcomp(0)-potEcomp(20)
558 ! Calculate the kinetic and total energy and the kinetic temperature
561 t_enegrad=t_enegrad+MPI_Wtime()-tt0
563 t_enegrad=t_enegrad+tcpu()-tt0
566 kinetic_T=2.0d0/(dimen*Rb)*EK
568 end subroutine brown_step
569 !-----------------------------------------------------------------------------
571 !-----------------------------------------------------------------------------
572 subroutine gauss(RO,AP,MT,M,N,*)
574 ! CALCULATES (RO**(-1))*AP BY GAUSS ELIMINATION
575 ! RO IS A SQUARE MATRIX
576 ! THE CALCULATED PRODUCT IS STORED IN AP
577 ! ABNORMAL EXIT IF RO IS SINGULAR
579 integer :: MT, M, N, M1,I,J,IM,&
581 real(kind=8) :: RO(MT,M),AP(MT,N),X,RM,PR,Y
587 if(dabs(X).le.1.0D-13) return 1
599 if(DABS(RO(J,I)).LE.RM) goto 2
613 if(dabs(X).le.1.0E-13) return 1
622 8 AP(J,K)=AP(J,K)-Y*AP(I,K)
624 9 RO(J,K)=RO(J,K)-Y*RO(I,K)
628 if(dabs(X).le.1.0E-13) return 1
638 15 X=X-AP(K,J)*RO(MI,K)
643 !-----------------------------------------------------------------------------
645 !-----------------------------------------------------------------------------
646 subroutine kinetic(KE_total)
647 !----------------------------------------------------------------
648 ! This subroutine calculates the total kinetic energy of the chain
649 !-----------------------------------------------------------------
651 ! implicit real*8 (a-h,o-z)
652 ! include 'DIMENSIONS'
653 ! include 'COMMON.VAR'
654 ! include 'COMMON.CHAIN'
655 ! include 'COMMON.DERIV'
656 ! include 'COMMON.GEO'
657 ! include 'COMMON.LOCAL'
658 ! include 'COMMON.INTERACT'
659 ! include 'COMMON.MD'
660 ! include 'COMMON.IOUNITS'
661 real(kind=8) :: KE_total,mscab
663 integer :: i,j,k,iti,mnum
664 real(kind=8) :: KEt_p,KEt_sc,KEr_p,KEr_sc,incr(3),&
667 write (iout,*) "Velocities, kietic"
669 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
670 (d_t(j,i+nres),j=1,3)
675 ! write (iout,*) "ISC",(isc(itype(i,1)),i=1,nres)
676 ! The translational part for peptide virtual bonds
682 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
683 ! write (iout,*) "Kinetic trp:",i,(incr(j),j=1,3)
685 v(j)=incr(j)+0.5d0*d_t(j,i)
687 vtot(i)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
688 KEt_p=KEt_p+mp(mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
690 incr(j)=incr(j)+d_t(j,i)
693 ! write(iout,*) 'KEt_p', KEt_p
694 ! The translational part for the side chain virtual bond
695 ! Only now we can initialize incr with zeros. It must be equal
696 ! to the velocities of the first Calpha.
702 iti=iabs(itype(i,mnum))
708 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
709 if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
710 .or.(mnum.eq.5)) then
716 v(j)=incr(j)+d_t(j,nres+i)
719 ! write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3)
720 ! write (iout,*) "i",i," msc",msc(iti)," v",(v(j),j=1,3)
721 KEt_sc=KEt_sc+mscab*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
722 vtot(i+nres)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
724 incr(j)=incr(j)+d_t(j,i)
728 ! write(iout,*) 'KEt_sc', KEt_sc
729 ! The part due to stretching and rotation of the peptide groups
733 ! write (iout,*) "i",i
734 ! write (iout,*) "i",i," mag1",mag1," mag2",mag2
738 ! write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3)
739 KEr_p=KEr_p+Ip(mnum)*(incr(1)*incr(1)+incr(2)*incr(2) &
743 ! write(iout,*) 'KEr_p', KEr_p
744 ! The rotational part of the side chain virtual bond
748 iti=iabs(itype(i,mnum))
749 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
750 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
751 .and.(mnum.ne.5)) then
753 incr(j)=d_t(j,nres+i)
755 ! write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3)
756 KEr_sc=KEr_sc+Isc(iti,mnum)*(incr(1)*incr(1)+incr(2)*incr(2)+ &
760 ! The total kinetic energy
762 ! write(iout,*) 'KEr_sc', KEr_sc
763 KE_total=0.5d0*(KEt_p+KEt_sc+0.25d0*KEr_p+KEr_sc)
764 ! write (iout,*) "KE_total",KE_total
766 end subroutine kinetic
767 !-----------------------------------------------------------------------------
769 !-----------------------------------------------------------------------------
771 !------------------------------------------------
772 ! The driver for molecular dynamics subroutines
773 !------------------------------------------------
776 use control, only:tcpu,ovrtim
777 ! use io_comm, only:ilen
779 use compare, only:secondary2,hairpin
780 use io, only:cartout,statout
781 ! implicit real*8 (a-h,o-z)
782 ! include 'DIMENSIONS'
785 integer :: IERROR,ERRCODE
787 ! include 'COMMON.SETUP'
788 ! include 'COMMON.CONTROL'
789 ! include 'COMMON.VAR'
790 ! include 'COMMON.MD'
792 ! include 'COMMON.LANGEVIN'
794 ! include 'COMMON.LANGEVIN.lang0'
796 ! include 'COMMON.CHAIN'
797 ! include 'COMMON.DERIV'
798 ! include 'COMMON.GEO'
799 ! include 'COMMON.LOCAL'
800 ! include 'COMMON.INTERACT'
801 ! include 'COMMON.IOUNITS'
802 ! include 'COMMON.NAMES'
803 ! include 'COMMON.TIME1'
804 ! include 'COMMON.HAIRPIN'
805 real(kind=8),dimension(3) :: L,vcm
807 real(kind=8),dimension(6*nres) :: v_work,v_transf !(maxres6) maxres6=6*maxres
809 integer :: rstcount !ilen,
811 character(len=50) :: tytul
812 !el common /gucio/ cm
813 integer :: itime,i,j,nharp
814 integer,dimension(4,nres/3) :: iharp !(4,nres/3)(4,maxres/3)
816 real(kind=8) :: tt0,scalfac
822 print *,"MY tmpdir",tmpdir,ilen(tmpdir)
823 if (ilen(tmpdir).gt.0) &
824 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_" &
825 //liczba(:ilen(liczba))//'.rst')
827 if (ilen(tmpdir).gt.0) &
828 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"//'.rst')
835 write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
841 print *,"just befor setup matix",nres
842 ! Determine the inverse of the inertia matrix.
843 call setup_MD_matrices
845 print *,"AFTER SETUP MATRICES"
847 print *,"AFTER INIT MD"
850 t_MDsetup = MPI_Wtime()-tt0
852 t_MDsetup = tcpu()-tt0
855 ! Entering the MD loop
861 if (lang.eq.2 .or. lang.eq.3) then
865 call sd_verlet_p_setup
867 call sd_verlet_ciccotti_setup
871 pfric0_mat(i,j,0)=pfric_mat(i,j)
872 afric0_mat(i,j,0)=afric_mat(i,j)
873 vfric0_mat(i,j,0)=vfric_mat(i,j)
874 prand0_mat(i,j,0)=prand_mat(i,j)
875 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
876 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
881 flag_stoch(i)=.false.
885 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
887 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
891 else if (lang.eq.1 .or. lang.eq.4) then
892 print *,"before setup_fricmat"
894 print *,"after setup_fricmat"
897 t_langsetup=MPI_Wtime()-tt0
900 t_langsetup=tcpu()-tt0
903 do itime=1,n_timestep
905 if (large.and. mod(itime,ntwe).eq.0) &
906 write (iout,*) "itime",itime
908 if (lang.gt.0 .and. surfarea .and. &
909 mod(itime,reset_fricmat).eq.0) then
910 if (lang.eq.2 .or. lang.eq.3) then
914 call sd_verlet_p_setup
916 call sd_verlet_ciccotti_setup
920 pfric0_mat(i,j,0)=pfric_mat(i,j)
921 afric0_mat(i,j,0)=afric_mat(i,j)
922 vfric0_mat(i,j,0)=vfric_mat(i,j)
923 prand0_mat(i,j,0)=prand_mat(i,j)
924 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
925 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
930 flag_stoch(i)=.false.
933 else if (lang.eq.1 .or. lang.eq.4) then
936 write (iout,'(a,i10)') &
937 "Friction matrix reset based on surface area, itime",itime
939 if (reset_vel .and. tbf .and. lang.eq.0 &
940 .and. mod(itime,count_reset_vel).eq.0) then
942 write(iout,'(a,f20.2)') &
943 "Velocities reset to random values, time",totT
946 d_t_old(j,i)=d_t(j,i)
950 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
954 d_t(j,0)=d_t(j,0)-vcm(j)
957 kinetic_T=2.0d0/(dimen3*Rb)*EK
958 scalfac=dsqrt(T_bath/kinetic_T)
959 write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
962 d_t_old(j,i)=scalfac*d_t(j,i)
968 ! Time-reversible RESPA algorithm
969 ! (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
970 call RESPA_step(itime)
972 ! Variable time step algorithm.
973 call velverlet_step(itime)
977 call brown_step(itime)
979 print *,"Brown dynamics not here!"
981 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
987 if (mod(itime,ntwe).eq.0) call statout(itime)
1001 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1.and.mnum.ne.5) then
1004 v_work(ind)=d_t(j,i+nres)
1009 write (66,'(80f10.5)') &
1010 ((d_t(j,i),j=1,3),i=0,nres-1),((d_t(j,i+nres),j=1,3),i=1,nres)
1014 v_transf(i)=v_transf(i)+gvec(j,i)*v_work(j)
1016 v_transf(i)= v_transf(i)*dsqrt(geigen(i))
1018 write (67,'(80f10.5)') (v_transf(i),i=1,ind)
1021 if (mod(itime,ntwx).eq.0) then
1022 write (tytul,'("time",f8.2)') totT
1024 call hairpin(.true.,nharp,iharp)
1025 call secondary2(.true.)
1026 call pdbout(potE,tytul,ipdb)
1031 if (rstcount.eq.1000.or.itime.eq.n_timestep) then
1032 open(irest2,file=rest2name,status='unknown')
1033 write(irest2,*) totT,EK,potE,totE,t_bath
1035 ! AL 4/17/17: Now writing d_t(0,:) too
1037 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
1039 ! AL 4/17/17: Now writing d_c(0,:) too
1041 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
1049 t_MD=MPI_Wtime()-tt0
1053 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))') &
1055 'MD calculations setup:',t_MDsetup,&
1056 'Energy & gradient evaluation:',t_enegrad,&
1057 'Stochastic MD setup:',t_langsetup,&
1058 'Stochastic MD step setup:',t_sdsetup,&
1060 write (iout,'(/28(1h=),a25,27(1h=))') &
1061 ' End of MD calculation '
1063 write (iout,*) "time for etotal",t_etotal," elong",t_elong,&
1065 write (iout,*) "time_fric",time_fric," time_stoch",time_stoch,&
1066 " time_fricmatmult",time_fricmatmult," time_fsample ",&
1071 !-----------------------------------------------------------------------------
1072 subroutine velverlet_step(itime)
1073 !-------------------------------------------------------------------------------
1074 ! Perform a single velocity Verlet step; the time step can be rescaled if
1075 ! increments in accelerations exceed the threshold
1076 !-------------------------------------------------------------------------------
1077 ! implicit real*8 (a-h,o-z)
1078 ! include 'DIMENSIONS'
1080 use control, only:tcpu
1084 integer :: ierror,ierrcode
1085 real(kind=8) :: errcode
1087 ! include 'COMMON.SETUP'
1088 ! include 'COMMON.VAR'
1089 ! include 'COMMON.MD'
1091 ! include 'COMMON.LANGEVIN'
1093 ! include 'COMMON.LANGEVIN.lang0'
1095 ! include 'COMMON.CHAIN'
1096 ! include 'COMMON.DERIV'
1097 ! include 'COMMON.GEO'
1098 ! include 'COMMON.LOCAL'
1099 ! include 'COMMON.INTERACT'
1100 ! include 'COMMON.IOUNITS'
1101 ! include 'COMMON.NAMES'
1102 ! include 'COMMON.TIME1'
1103 ! include 'COMMON.MUCA'
1104 real(kind=8),dimension(3) :: vcm,incr
1105 real(kind=8),dimension(3) :: L
1106 integer :: count,rstcount !ilen,
1108 character(len=50) :: tytul
1109 integer :: maxcount_scale = 20
1110 !el common /gucio/ cm
1111 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1112 !el common /stochcalc/ stochforcvec
1113 integer :: itime,icount_scale,itime_scal,i,j,ifac_time
1115 real(kind=8) :: epdrift,tt0,fac_time
1117 if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres)) !(MAXRES6) maxres6=6*maxres
1123 else if (lang.eq.2 .or. lang.eq.3) then
1125 call stochastic_force(stochforcvec)
1128 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
1130 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1137 icount_scale=icount_scale+1
1138 if (icount_scale.gt.maxcount_scale) then
1140 "ERROR: too many attempts at scaling down the time step. ",&
1141 "amax=",amax,"epdrift=",epdrift,&
1142 "damax=",damax,"edriftmax=",edriftmax,&
1146 call MPI_Abort(MPI_COMM_WORLD,IERROR,IERRCODE)
1150 ! First step of the velocity Verlet algorithm
1155 else if (lang.eq.3) then
1157 call sd_verlet1_ciccotti
1159 else if (lang.eq.1) then
1164 ! Build the chain from the newly calculated coordinates
1165 call chainbuild_cart
1166 if (rattle) call rattle1
1168 if (large.and. mod(itime,ntwe).eq.0) then
1169 write (iout,*) "Cartesian and internal coordinates: step 1"
1174 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),&
1175 (dc(j,i+nres),j=1,3)
1177 write (iout,*) "Accelerations"
1179 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1180 (d_a(j,i+nres),j=1,3)
1182 write (iout,*) "Velocities, step 1"
1184 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1185 (d_t(j,i+nres),j=1,3)
1194 ! Calculate energy and forces
1196 call etotal(potEcomp)
1197 ! AL 4/17/17: Reduce the steps if NaNs occurred.
1198 if (potEcomp(0).gt.0.99e20 .or. isnan(potEcomp(0)).gt.0) then
1203 if (large.and. mod(itime,ntwe).eq.0) &
1204 call enerprint(potEcomp)
1207 t_etotal=t_etotal+MPI_Wtime()-tt0
1209 t_etotal=t_etotal+tcpu()-tt0
1212 potE=potEcomp(0)-potEcomp(20)
1214 ! Get the new accelerations
1217 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1219 t_enegrad=t_enegrad+tcpu()-tt0
1221 ! Determine maximum acceleration and scale down the timestep if needed
1223 amax=amax/(itime_scal+1)**2
1224 call predict_edrift(epdrift)
1225 if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then
1226 ! Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
1228 ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax)) &
1230 itime_scal=itime_scal+ifac_time
1231 ! fac_time=dmin1(damax/amax,0.5d0)
1232 fac_time=0.5d0**ifac_time
1233 d_time=d_time*fac_time
1234 if (lang.eq.2 .or. lang.eq.3) then
1236 ! write (iout,*) "Calling sd_verlet_setup: 1"
1237 ! Rescale the stochastic forces and recalculate or restore
1238 ! the matrices of tinker integrator
1239 if (itime_scal.gt.maxflag_stoch) then
1240 if (large) write (iout,'(a,i5,a)') &
1241 "Calculate matrices for stochastic step;",&
1242 " itime_scal ",itime_scal
1244 call sd_verlet_p_setup
1246 call sd_verlet_ciccotti_setup
1248 write (iout,'(2a,i3,a,i3,1h.)') &
1249 "Warning: cannot store matrices for stochastic",&
1250 " integration because the index",itime_scal,&
1251 " is greater than",maxflag_stoch
1252 write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",&
1253 " integration Langevin algorithm for better efficiency."
1254 else if (flag_stoch(itime_scal)) then
1255 if (large) write (iout,'(a,i5,a,l1)') &
1256 "Restore matrices for stochastic step; itime_scal ",&
1257 itime_scal," flag ",flag_stoch(itime_scal)
1260 pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
1261 afric_mat(i,j)=afric0_mat(i,j,itime_scal)
1262 vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
1263 prand_mat(i,j)=prand0_mat(i,j,itime_scal)
1264 vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
1265 vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
1269 if (large) write (iout,'(2a,i5,a,l1)') &
1270 "Calculate & store matrices for stochastic step;",&
1271 " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
1273 call sd_verlet_p_setup
1275 call sd_verlet_ciccotti_setup
1277 flag_stoch(ifac_time)=.true.
1280 pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
1281 afric0_mat(i,j,itime_scal)=afric_mat(i,j)
1282 vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
1283 prand0_mat(i,j,itime_scal)=prand_mat(i,j)
1284 vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
1285 vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
1289 fac_time=1.0d0/dsqrt(fac_time)
1291 stochforcvec(i)=fac_time*stochforcvec(i)
1294 else if (lang.eq.1) then
1295 ! Rescale the accelerations due to stochastic forces
1296 fac_time=1.0d0/dsqrt(fac_time)
1298 d_as_work(i)=d_as_work(i)*fac_time
1301 if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)') &
1302 "itime",itime," Timestep scaled down to ",&
1303 d_time," ifac_time",ifac_time," itime_scal",itime_scal
1305 ! Second step of the velocity Verlet algorithm
1310 else if (lang.eq.3) then
1312 call sd_verlet2_ciccotti
1314 else if (lang.eq.1) then
1319 if (rattle) call rattle2
1322 if (d_time.ne.d_time0) then
1325 if (lang.eq.2 .or. lang.eq.3) then
1326 if (large) write (iout,'(a)') &
1327 "Restore original matrices for stochastic step"
1328 ! write (iout,*) "Calling sd_verlet_setup: 2"
1329 ! Restore the matrices of tinker integrator if the time step has been restored
1332 pfric_mat(i,j)=pfric0_mat(i,j,0)
1333 afric_mat(i,j)=afric0_mat(i,j,0)
1334 vfric_mat(i,j)=vfric0_mat(i,j,0)
1335 prand_mat(i,j)=prand0_mat(i,j,0)
1336 vrand_mat1(i,j)=vrand0_mat1(i,j,0)
1337 vrand_mat2(i,j)=vrand0_mat2(i,j,0)
1346 ! Calculate the kinetic and the total energy and the kinetic temperature
1350 ! call kinetic1(EK1)
1351 ! write (iout,*) "step",itime," EK",EK," EK1",EK1
1353 ! Couple the system to Berendsen bath if needed
1354 if (tbf .and. lang.eq.0) then
1357 kinetic_T=2.0d0/(dimen3*Rb)*EK
1358 ! Backup the coordinates, velocities, and accelerations
1362 d_t_old(j,i)=d_t(j,i)
1363 d_a_old(j,i)=d_a(j,i)
1367 if (mod(itime,ntwe).eq.0 .and. large) then
1368 write (iout,*) "Velocities, step 2"
1370 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1371 (d_t(j,i+nres),j=1,3)
1376 end subroutine velverlet_step
1377 !-----------------------------------------------------------------------------
1378 subroutine RESPA_step(itime)
1379 !-------------------------------------------------------------------------------
1380 ! Perform a single RESPA step.
1381 !-------------------------------------------------------------------------------
1382 ! implicit real*8 (a-h,o-z)
1383 ! include 'DIMENSIONS'
1387 use control, only:tcpu
1389 ! use io_conf, only:cartprint
1392 integer :: IERROR,ERRCODE
1394 ! include 'COMMON.SETUP'
1395 ! include 'COMMON.CONTROL'
1396 ! include 'COMMON.VAR'
1397 ! include 'COMMON.MD'
1399 ! include 'COMMON.LANGEVIN'
1401 ! include 'COMMON.LANGEVIN.lang0'
1403 ! include 'COMMON.CHAIN'
1404 ! include 'COMMON.DERIV'
1405 ! include 'COMMON.GEO'
1406 ! include 'COMMON.LOCAL'
1407 ! include 'COMMON.INTERACT'
1408 ! include 'COMMON.IOUNITS'
1409 ! include 'COMMON.NAMES'
1410 ! include 'COMMON.TIME1'
1411 real(kind=8),dimension(0:n_ene) :: energia_short,energia_long
1412 real(kind=8),dimension(3) :: L,vcm,incr
1413 real(kind=8),dimension(3,0:2*nres) :: dc_old0,d_t_old0,d_a_old0 !(3,0:maxres2) maxres2=2*maxres
1414 logical :: PRINT_AMTS_MSG = .false.
1415 integer :: count,rstcount !ilen,
1417 character(len=50) :: tytul
1418 integer :: maxcount_scale = 10
1419 !el common /gucio/ cm
1420 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1421 !el common /stochcalc/ stochforcvec
1422 integer :: itime,itt,i,j,itsplit
1424 !el common /cipiszcze/ itt
1426 real(kind=8) :: epdrift,tt0,epdriftmax
1429 if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres)) !(MAXRES6) maxres6=6*maxres
1433 if (large.and. mod(itime,ntwe).eq.0) then
1434 write (iout,*) "***************** RESPA itime",itime
1435 write (iout,*) "Cartesian and internal coordinates: step 0"
1437 call pdbout(0.0d0,"cipiszcze",iout)
1439 write (iout,*) "Accelerations from long-range forces"
1441 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1442 (d_a(j,i+nres),j=1,3)
1444 write (iout,*) "Velocities, step 0"
1446 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1447 (d_t(j,i+nres),j=1,3)
1452 ! Perform the initial RESPA step (increment velocities)
1453 ! write (iout,*) "*********************** RESPA ini"
1456 if (mod(itime,ntwe).eq.0 .and. large) then
1457 write (iout,*) "Velocities, end"
1459 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1460 (d_t(j,i+nres),j=1,3)
1464 ! Compute the short-range forces
1470 ! 7/2/2009 commented out
1472 ! call etotal_short(energia_short)
1475 ! 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
1478 d_a(j,i)=d_a_short(j,i)
1482 if (large.and. mod(itime,ntwe).eq.0) then
1483 write (iout,*) "energia_short",energia_short(0)
1484 write (iout,*) "Accelerations from short-range forces"
1486 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1487 (d_a(j,i+nres),j=1,3)
1492 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1494 t_enegrad=t_enegrad+tcpu()-tt0
1499 d_t_old(j,i)=d_t(j,i)
1500 d_a_old(j,i)=d_a(j,i)
1503 ! 6/30/08 A-MTS: attempt at increasing the split number
1506 dc_old0(j,i)=dc_old(j,i)
1507 d_t_old0(j,i)=d_t_old(j,i)
1508 d_a_old0(j,i)=d_a_old(j,i)
1511 if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
1512 if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
1519 ! write (iout,*) "itime",itime," ntime_split",ntime_split
1520 ! Split the time step
1521 d_time=d_time0/ntime_split
1522 ! Perform the short-range RESPA steps (velocity Verlet increments of
1523 ! positions and velocities using short-range forces)
1524 ! write (iout,*) "*********************** RESPA split"
1525 do itsplit=1,ntime_split
1528 else if (lang.eq.2 .or. lang.eq.3) then
1530 call stochastic_force(stochforcvec)
1533 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
1535 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1540 ! First step of the velocity Verlet algorithm
1545 else if (lang.eq.3) then
1547 call sd_verlet1_ciccotti
1549 else if (lang.eq.1) then
1554 ! Build the chain from the newly calculated coordinates
1555 call chainbuild_cart
1556 if (rattle) call rattle1
1558 if (large.and. mod(itime,ntwe).eq.0) then
1559 write (iout,*) "***** ITSPLIT",itsplit
1560 write (iout,*) "Cartesian and internal coordinates: step 1"
1561 call pdbout(0.0d0,"cipiszcze",iout)
1564 write (iout,*) "Velocities, step 1"
1566 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1567 (d_t(j,i+nres),j=1,3)
1576 ! Calculate energy and forces
1578 call etotal_short(energia_short)
1579 ! AL 4/17/17: Exit itime_split loop when energy goes infinite
1580 if (energia_short(0).gt.0.99e20 .or. isnan(energia_short(0)) ) then
1581 if (PRINT_AMTS_MSG) &
1582 write (iout,*) "Infinities/NaNs in energia_short",energia_short(0),"; increasing ntime_split to",ntime_split
1583 ntime_split=ntime_split*2
1584 if (ntime_split.gt.maxtime_split) then
1587 "Cannot rescue the run; aborting job. Retry with a smaller time step"
1589 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1592 "Cannot rescue the run; terminating. Retry with a smaller time step"
1598 if (large.and. mod(itime,ntwe).eq.0) &
1599 call enerprint(energia_short)
1602 t_eshort=t_eshort+MPI_Wtime()-tt0
1604 t_eshort=t_eshort+tcpu()-tt0
1608 ! Get the new accelerations
1610 ! 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
1613 d_a_short(j,i)=d_a(j,i)
1617 if (large.and. mod(itime,ntwe).eq.0) then
1618 write (iout,*)"energia_short",energia_short(0)
1619 write (iout,*) "Accelerations from short-range forces"
1621 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1622 (d_a(j,i+nres),j=1,3)
1627 ! Determine maximum acceleration and scale down the timestep if needed
1629 amax=amax/ntime_split**2
1630 call predict_edrift(epdrift)
1631 if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0) &
1632 write (iout,*) "amax",amax," damax",damax,&
1633 " epdrift",epdrift," epdriftmax",epdriftmax
1634 ! Exit loop and try with increased split number if the change of
1635 ! acceleration is too big
1636 if (amax.gt.damax .or. epdrift.gt.edriftmax) then
1637 if (ntime_split.lt.maxtime_split) then
1639 ntime_split=ntime_split*2
1640 ! AL 4/17/17: We should exit the itime_split loop when acceleration change is too big
1644 dc_old(j,i)=dc_old0(j,i)
1645 d_t_old(j,i)=d_t_old0(j,i)
1646 d_a_old(j,i)=d_a_old0(j,i)
1649 if (PRINT_AMTS_MSG) then
1650 write (iout,*) "acceleration/energy drift too large",amax,&
1651 epdrift," split increased to ",ntime_split," itime",itime,&
1657 "Uh-hu. Bumpy landscape. Maximum splitting number",&
1659 " already reached!!! Trying to carry on!"
1663 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1665 t_enegrad=t_enegrad+tcpu()-tt0
1667 ! Second step of the velocity Verlet algorithm
1672 else if (lang.eq.3) then
1674 call sd_verlet2_ciccotti
1676 else if (lang.eq.1) then
1681 if (rattle) call rattle2
1682 ! Backup the coordinates, velocities, and accelerations
1686 d_t_old(j,i)=d_t(j,i)
1687 d_a_old(j,i)=d_a(j,i)
1694 ! Restore the time step
1696 ! Compute long-range forces
1703 call etotal_long(energia_long)
1704 if (energia_long(0).gt.0.99e20 .or. isnan(energia_long(0))) then
1707 "Infinitied/NaNs in energia_long, Aborting MPI job."
1709 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1711 write (iout,*) "Infinitied/NaNs in energia_long, terminating."
1715 if (large.and. mod(itime,ntwe).eq.0) &
1716 call enerprint(energia_long)
1719 t_elong=t_elong+MPI_Wtime()-tt0
1721 t_elong=t_elong+tcpu()-tt0
1727 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1729 t_enegrad=t_enegrad+tcpu()-tt0
1731 ! Compute accelerations from long-range forces
1733 if (large.and. mod(itime,ntwe).eq.0) then
1734 write (iout,*) "energia_long",energia_long(0)
1735 write (iout,*) "Cartesian and internal coordinates: step 2"
1737 call pdbout(0.0d0,"cipiszcze",iout)
1739 write (iout,*) "Accelerations from long-range forces"
1741 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1742 (d_a(j,i+nres),j=1,3)
1744 write (iout,*) "Velocities, step 2"
1746 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1747 (d_t(j,i+nres),j=1,3)
1751 ! Compute the final RESPA step (increment velocities)
1752 ! write (iout,*) "*********************** RESPA fin"
1754 ! Compute the complete potential energy
1756 potEcomp(i)=energia_short(i)+energia_long(i)
1758 potE=potEcomp(0)-potEcomp(20)
1759 ! potE=energia_short(0)+energia_long(0)
1762 ! Calculate the kinetic and the total energy and the kinetic temperature
1765 ! Couple the system to Berendsen bath if needed
1766 if (tbf .and. lang.eq.0) then
1769 kinetic_T=2.0d0/(dimen3*Rb)*EK
1770 ! Backup the coordinates, velocities, and accelerations
1772 if (mod(itime,ntwe).eq.0 .and. large) then
1773 write (iout,*) "Velocities, end"
1775 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1776 (d_t(j,i+nres),j=1,3)
1781 end subroutine RESPA_step
1782 !-----------------------------------------------------------------------------
1783 subroutine RESPA_vel
1784 ! First and last RESPA step (incrementing velocities using long-range
1787 ! implicit real*8 (a-h,o-z)
1788 ! include 'DIMENSIONS'
1789 ! include 'COMMON.CONTROL'
1790 ! include 'COMMON.VAR'
1791 ! include 'COMMON.MD'
1792 ! include 'COMMON.CHAIN'
1793 ! include 'COMMON.DERIV'
1794 ! include 'COMMON.GEO'
1795 ! include 'COMMON.LOCAL'
1796 ! include 'COMMON.INTERACT'
1797 ! include 'COMMON.IOUNITS'
1798 ! include 'COMMON.NAMES'
1799 integer :: i,j,inres,mnum
1802 d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
1806 d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
1811 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
1812 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
1813 .and.(mnum.ne.5)) then
1816 d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
1821 end subroutine RESPA_vel
1822 !-----------------------------------------------------------------------------
1824 ! Applying velocity Verlet algorithm - step 1 to coordinates
1826 ! implicit real*8 (a-h,o-z)
1827 ! include 'DIMENSIONS'
1828 ! include 'COMMON.CONTROL'
1829 ! include 'COMMON.VAR'
1830 ! include 'COMMON.MD'
1831 ! include 'COMMON.CHAIN'
1832 ! include 'COMMON.DERIV'
1833 ! include 'COMMON.GEO'
1834 ! include 'COMMON.LOCAL'
1835 ! include 'COMMON.INTERACT'
1836 ! include 'COMMON.IOUNITS'
1837 ! include 'COMMON.NAMES'
1838 real(kind=8) :: adt,adt2
1839 integer :: i,j,inres,mnum
1842 write (iout,*) "VELVERLET1 START: DC"
1844 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),&
1845 (dc(j,i+nres),j=1,3)
1849 adt=d_a_old(j,0)*d_time
1851 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1852 d_t_new(j,0)=d_t_old(j,0)+adt2
1853 d_t(j,0)=d_t_old(j,0)+adt
1857 adt=d_a_old(j,i)*d_time
1859 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1860 d_t_new(j,i)=d_t_old(j,i)+adt2
1861 d_t(j,i)=d_t_old(j,i)+adt
1866 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
1867 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
1868 .and.(mnum.ne.5)) then
1871 adt=d_a_old(j,inres)*d_time
1873 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1874 d_t_new(j,inres)=d_t_old(j,inres)+adt2
1875 d_t(j,inres)=d_t_old(j,inres)+adt
1880 write (iout,*) "VELVERLET1 END: DC"
1882 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),&
1883 (dc(j,i+nres),j=1,3)
1887 end subroutine verlet1
1888 !-----------------------------------------------------------------------------
1890 ! Step 2 of the velocity Verlet algorithm: update velocities
1892 ! implicit real*8 (a-h,o-z)
1893 ! include 'DIMENSIONS'
1894 ! include 'COMMON.CONTROL'
1895 ! include 'COMMON.VAR'
1896 ! include 'COMMON.MD'
1897 ! include 'COMMON.CHAIN'
1898 ! include 'COMMON.DERIV'
1899 ! include 'COMMON.GEO'
1900 ! include 'COMMON.LOCAL'
1901 ! include 'COMMON.INTERACT'
1902 ! include 'COMMON.IOUNITS'
1903 ! include 'COMMON.NAMES'
1904 integer :: i,j,inres,mnum
1907 d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
1911 d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
1916 ! iti=iabs(itype(i,mnum))
1917 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
1918 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
1919 .and.(mnum.ne.5)) then
1922 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
1927 end subroutine verlet2
1928 !-----------------------------------------------------------------------------
1929 subroutine sddir_precalc
1930 ! Applying velocity Verlet algorithm - step 1 to coordinates
1931 ! implicit real*8 (a-h,o-z)
1932 ! include 'DIMENSIONS'
1938 ! include 'COMMON.CONTROL'
1939 ! include 'COMMON.VAR'
1940 ! include 'COMMON.MD'
1942 ! include 'COMMON.LANGEVIN'
1944 ! include 'COMMON.LANGEVIN.lang0'
1946 ! include 'COMMON.CHAIN'
1947 ! include 'COMMON.DERIV'
1948 ! include 'COMMON.GEO'
1949 ! include 'COMMON.LOCAL'
1950 ! include 'COMMON.INTERACT'
1951 ! include 'COMMON.IOUNITS'
1952 ! include 'COMMON.NAMES'
1953 ! include 'COMMON.TIME1'
1954 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1955 !el common /stochcalc/ stochforcvec
1956 real(kind=8) :: time00
1958 ! Compute friction and stochastic forces
1963 time_fric=time_fric+MPI_Wtime()-time00
1965 call stochastic_force(stochforcvec)
1966 time_stoch=time_stoch+MPI_Wtime()-time00
1969 ! Compute the acceleration due to friction forces (d_af_work) and stochastic
1970 ! forces (d_as_work)
1972 call ginv_mult(fric_work, d_af_work)
1973 call ginv_mult(stochforcvec, d_as_work)
1975 end subroutine sddir_precalc
1976 !-----------------------------------------------------------------------------
1977 subroutine sddir_verlet1
1978 ! Applying velocity Verlet algorithm - step 1 to velocities
1981 ! implicit real*8 (a-h,o-z)
1982 ! include 'DIMENSIONS'
1983 ! include 'COMMON.CONTROL'
1984 ! include 'COMMON.VAR'
1985 ! include 'COMMON.MD'
1987 ! include 'COMMON.LANGEVIN'
1989 ! include 'COMMON.LANGEVIN.lang0'
1991 ! include 'COMMON.CHAIN'
1992 ! include 'COMMON.DERIV'
1993 ! include 'COMMON.GEO'
1994 ! include 'COMMON.LOCAL'
1995 ! include 'COMMON.INTERACT'
1996 ! include 'COMMON.IOUNITS'
1997 ! include 'COMMON.NAMES'
1998 ! Revised 3/31/05 AL: correlation between random contributions to
1999 ! position and velocity increments included.
2000 real(kind=8) :: sqrt13 = 0.57735026918962576451d0 ! 1/sqrt(3)
2001 real(kind=8) :: adt,adt2
2002 integer :: i,j,ind,inres,mnum
2004 ! Add the contribution from BOTH friction and stochastic force to the
2005 ! coordinates, but ONLY the contribution from the friction forces to velocities
2008 adt=(d_a_old(j,0)+d_af_work(j))*d_time
2009 adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
2010 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
2011 d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
2012 d_t(j,0)=d_t_old(j,0)+adt
2017 adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
2018 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
2019 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
2020 d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
2021 d_t(j,i)=d_t_old(j,i)+adt
2027 ! iti=iabs(itype(i,mnum))
2028 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2029 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2030 .and.(mnum.ne.5)) then
2033 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
2034 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
2035 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
2036 d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
2037 d_t(j,inres)=d_t_old(j,inres)+adt
2043 end subroutine sddir_verlet1
2044 !-----------------------------------------------------------------------------
2045 subroutine sddir_verlet2
2046 ! Calculating the adjusted velocities for accelerations
2049 ! implicit real*8 (a-h,o-z)
2050 ! include 'DIMENSIONS'
2051 ! include 'COMMON.CONTROL'
2052 ! include 'COMMON.VAR'
2053 ! include 'COMMON.MD'
2055 ! include 'COMMON.LANGEVIN'
2057 ! include 'COMMON.LANGEVIN.lang0'
2059 ! include 'COMMON.CHAIN'
2060 ! include 'COMMON.DERIV'
2061 ! include 'COMMON.GEO'
2062 ! include 'COMMON.LOCAL'
2063 ! include 'COMMON.INTERACT'
2064 ! include 'COMMON.IOUNITS'
2065 ! include 'COMMON.NAMES'
2066 real(kind=8),dimension(6*nres) :: stochforcvec,d_as_work1 !(MAXRES6) maxres6=6*maxres
2067 real(kind=8) :: cos60 = 0.5d0, sin60 = 0.86602540378443864676d0
2068 integer :: i,j,ind,inres,mnum
2069 ! Revised 3/31/05 AL: correlation between random contributions to
2070 ! position and velocity increments included.
2071 ! The correlation coefficients are calculated at low-friction limit.
2072 ! Also, friction forces are now not calculated with new velocities.
2074 ! call friction_force
2075 call stochastic_force(stochforcvec)
2077 ! Compute the acceleration due to friction forces (d_af_work) and stochastic
2078 ! forces (d_as_work)
2080 call ginv_mult(stochforcvec, d_as_work1)
2086 d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j)) &
2087 +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
2092 d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j)) &
2093 +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
2099 ! iti=iabs(itype(i,mnum))
2100 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2101 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2102 .and.(mnum.ne.5)) then
2105 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres) &
2106 +d_af_work(ind+j))+sin60*d_as_work(ind+j) &
2107 +cos60*d_as_work1(ind+j))*d_time
2113 end subroutine sddir_verlet2
2114 !-----------------------------------------------------------------------------
2115 subroutine max_accel
2117 ! Find the maximum difference in the accelerations of the the sites
2118 ! at the beginning and the end of the time step.
2121 ! implicit real*8 (a-h,o-z)
2122 ! include 'DIMENSIONS'
2123 ! include 'COMMON.CONTROL'
2124 ! include 'COMMON.VAR'
2125 ! include 'COMMON.MD'
2126 ! include 'COMMON.CHAIN'
2127 ! include 'COMMON.DERIV'
2128 ! include 'COMMON.GEO'
2129 ! include 'COMMON.LOCAL'
2130 ! include 'COMMON.INTERACT'
2131 ! include 'COMMON.IOUNITS'
2132 real(kind=8),dimension(3) :: aux,accel,accel_old
2133 real(kind=8) :: dacc
2137 ! aux(j)=d_a(j,0)-d_a_old(j,0)
2138 accel_old(j)=d_a_old(j,0)
2145 ! 7/3/08 changed to asymmetric difference
2147 ! accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
2148 accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
2149 accel(j)=accel(j)+0.5d0*d_a(j,i)
2150 ! if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
2151 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
2152 dacc=dabs(accel(j)-accel_old(j))
2153 ! write (iout,*) i,dacc
2154 if (dacc.gt.amax) amax=dacc
2162 accel_old(j)=d_a_old(j,0)
2167 accel_old(j)=accel_old(j)+d_a_old(j,1)
2168 accel(j)=accel(j)+d_a(j,1)
2173 ! iti=iabs(itype(i,mnum))
2174 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2175 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2176 .and.(mnum.ne.5)) then
2178 ! accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
2179 accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
2180 accel(j)=accel(j)+d_a(j,i+nres)
2184 ! if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
2185 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
2186 dacc=dabs(accel(j)-accel_old(j))
2187 ! write (iout,*) "side-chain",i,dacc
2188 if (dacc.gt.amax) amax=dacc
2192 accel_old(j)=accel_old(j)+d_a_old(j,i)
2193 accel(j)=accel(j)+d_a(j,i)
2194 ! aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
2198 end subroutine max_accel
2199 !-----------------------------------------------------------------------------
2200 subroutine predict_edrift(epdrift)
2202 ! Predict the drift of the potential energy
2205 use control_data, only: lmuca
2206 ! implicit real*8 (a-h,o-z)
2207 ! include 'DIMENSIONS'
2208 ! include 'COMMON.CONTROL'
2209 ! include 'COMMON.VAR'
2210 ! include 'COMMON.MD'
2211 ! include 'COMMON.CHAIN'
2212 ! include 'COMMON.DERIV'
2213 ! include 'COMMON.GEO'
2214 ! include 'COMMON.LOCAL'
2215 ! include 'COMMON.INTERACT'
2216 ! include 'COMMON.IOUNITS'
2217 ! include 'COMMON.MUCA'
2218 real(kind=8) :: epdrift,epdriftij
2220 ! Drift of the potential energy
2226 epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
2227 if (lmuca) epdriftij=epdriftij*factor
2228 ! write (iout,*) "back",i,j,epdriftij
2229 if (epdriftij.gt.epdrift) epdrift=epdriftij
2233 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1.and.&
2234 molnum(i).ne.5) then
2237 dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
2238 if (lmuca) epdriftij=epdriftij*factor
2239 ! write (iout,*) "side",i,j,epdriftij
2240 if (epdriftij.gt.epdrift) epdrift=epdriftij
2244 epdrift=0.5d0*epdrift*d_time*d_time
2245 ! write (iout,*) "epdrift",epdrift
2247 end subroutine predict_edrift
2248 !-----------------------------------------------------------------------------
2249 subroutine verlet_bath
2251 ! Coupling to the thermostat by using the Berendsen algorithm
2254 ! implicit real*8 (a-h,o-z)
2255 ! include 'DIMENSIONS'
2256 ! include 'COMMON.CONTROL'
2257 ! include 'COMMON.VAR'
2258 ! include 'COMMON.MD'
2259 ! include 'COMMON.CHAIN'
2260 ! include 'COMMON.DERIV'
2261 ! include 'COMMON.GEO'
2262 ! include 'COMMON.LOCAL'
2263 ! include 'COMMON.INTERACT'
2264 ! include 'COMMON.IOUNITS'
2265 ! include 'COMMON.NAMES'
2266 real(kind=8) :: T_half,fact
2267 integer :: i,j,inres,mnum
2269 T_half=2.0d0/(dimen3*Rb)*EK
2270 fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
2271 ! write(iout,*) "T_half", T_half
2272 ! write(iout,*) "EK", EK
2273 ! write(iout,*) "fact", fact
2275 d_t(j,0)=fact*d_t(j,0)
2279 d_t(j,i)=fact*d_t(j,i)
2284 ! iti=iabs(itype(i,mnum))
2285 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2286 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2287 .and.(mnum.ne.5)) then
2290 d_t(j,inres)=fact*d_t(j,inres)
2295 end subroutine verlet_bath
2296 !-----------------------------------------------------------------------------
2298 ! Set up the initial conditions of a MD simulation
2301 use control, only:tcpu
2302 !el use io_basic, only:ilen
2305 use minimm, only:minim_dc,minimize,sc_move
2306 use io_config, only:readrst
2307 use io, only:statout
2308 ! implicit real*8 (a-h,o-z)
2309 ! include 'DIMENSIONS'
2312 character(len=16) :: form
2313 integer :: IERROR,ERRCODE
2315 ! include 'COMMON.SETUP'
2316 ! include 'COMMON.CONTROL'
2317 ! include 'COMMON.VAR'
2318 ! include 'COMMON.MD'
2320 ! include 'COMMON.LANGEVIN'
2322 ! include 'COMMON.LANGEVIN.lang0'
2324 ! include 'COMMON.CHAIN'
2325 ! include 'COMMON.DERIV'
2326 ! include 'COMMON.GEO'
2327 ! include 'COMMON.LOCAL'
2328 ! include 'COMMON.INTERACT'
2329 ! include 'COMMON.IOUNITS'
2330 ! include 'COMMON.NAMES'
2331 ! include 'COMMON.REMD'
2332 real(kind=8),dimension(0:n_ene) :: energia_long,energia_short
2333 real(kind=8),dimension(3) :: vcm,incr,L
2334 real(kind=8) :: xv,sigv,lowb,highb
2335 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
2336 character(len=256) :: qstr
2339 character(len=50) :: tytul
2340 logical :: file_exist
2341 !el common /gucio/ cm
2342 integer :: i,j,ipos,iq,iw,nft_sc,iretcode,nfun,itime,ierr,mnum
2343 real(kind=8) :: etot,tt0
2347 ! write(iout,*) "d_time", d_time
2348 ! Compute the standard deviations of stochastic forces for Langevin dynamics
2349 ! if the friction coefficients do not depend on surface area
2350 if (lang.gt.0 .and. .not.surfarea) then
2353 stdforcp(i)=stdfp(mnum)*dsqrt(gamp(mnum))
2357 stdforcsc(i)=stdfsc(iabs(itype(i,mnum)),mnum) &
2358 *dsqrt(gamsc(iabs(itype(i,mnum)),mnum))
2361 ! Open the pdb file for snapshotshots
2364 if (ilen(tmpdir).gt.0) &
2365 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2366 liczba(:ilen(liczba))//".pdb")
2368 file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2372 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) &
2373 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2374 liczba(:ilen(liczba))//".x")
2375 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2378 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) &
2379 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2380 liczba(:ilen(liczba))//".cx")
2381 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2387 if (ilen(tmpdir).gt.0) &
2388 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
2389 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
2391 if (ilen(tmpdir).gt.0) &
2392 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
2393 cartname=prefix(:ilen(prefix))//"_MD.cx"
2397 write (qstr,'(256(1h ))')
2400 iq = qinfrag(i,iset)*10
2401 iw = wfrag(i,iset)/100
2403 if(me.eq.king.or..not.out1file) &
2404 write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
2405 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
2410 iq = qinpair(i,iset)*10
2411 iw = wpair(i,iset)/100
2413 if(me.eq.king.or..not.out1file) &
2414 write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
2415 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
2419 ! pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
2421 ! cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
2423 ! cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
2425 ! statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
2429 if (restart1file) then
2431 inquire(file=mremd_rst_name,exist=file_exist)
2432 write (*,*) me," Before broadcast: file_exist",file_exist
2434 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,&
2437 write (*,*) me," After broadcast: file_exist",file_exist
2438 ! inquire(file=mremd_rst_name,exist=file_exist)
2439 if(me.eq.king.or..not.out1file) &
2440 write(iout,*) "Initial state read by master and distributed"
2442 if (ilen(tmpdir).gt.0) &
2443 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_' &
2444 //liczba(:ilen(liczba))//'.rst')
2445 inquire(file=rest2name,exist=file_exist)
2448 if(.not.restart1file) then
2449 if(me.eq.king.or..not.out1file) &
2450 write(iout,*) "Initial state will be read from file ",&
2451 rest2name(:ilen(rest2name))
2454 call rescale_weights(t_bath)
2456 if(me.eq.king.or..not.out1file)then
2457 if (restart1file) then
2458 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),&
2461 write(iout,*) "File ",rest2name(:ilen(rest2name)),&
2464 write(iout,*) "Initial velocities randomly generated"
2471 ! Generate initial velocities
2472 if(me.eq.king.or..not.out1file) &
2473 write(iout,*) "Initial velocities randomly generated"
2478 ! rest2name = prefix(:ilen(prefix))//'.rst'
2479 if(me.eq.king.or..not.out1file)then
2480 write (iout,*) "Initial velocities"
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 ! Zeroing the total angular momentum of the system
2486 write(iout,*) "Calling the zero-angular momentum subroutine"
2489 ! Getting the potential energy and forces and velocities and accelerations
2491 ! write (iout,*) "velocity of the center of the mass:"
2492 ! write (iout,*) (vcm(j),j=1,3)
2494 d_t(j,0)=d_t(j,0)-vcm(j)
2496 ! Removing the velocity of the center of mass
2498 if(me.eq.king.or..not.out1file)then
2499 write (iout,*) "vcm right after adjustment:"
2500 write (iout,*) (vcm(j),j=1,3)
2502 if ((.not.rest).and.(indpdb.eq.0)) then
2504 if(iranconf.ne.0) then
2506 print *, 'Calling OVERLAP_SC'
2507 call overlap_sc(fail)
2510 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
2511 print *,'SC_move',nft_sc,etot
2512 if(me.eq.king.or..not.out1file) &
2513 write(iout,*) 'SC_move',nft_sc,etot
2517 print *, 'Calling MINIM_DC'
2518 call minim_dc(etot,iretcode,nfun)
2520 call geom_to_var(nvar,varia)
2521 print *,'Calling MINIMIZE.'
2522 call minimize(etot,varia,iretcode,nfun)
2523 call var_to_geom(nvar,varia)
2525 if(me.eq.king.or..not.out1file) &
2526 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2529 call chainbuild_cart
2534 kinetic_T=2.0d0/(dimen3*Rb)*EK
2535 if(me.eq.king.or..not.out1file)then
2545 call etotal(potEcomp)
2546 if (large) call enerprint(potEcomp)
2549 t_etotal=t_etotal+MPI_Wtime()-tt0
2551 t_etotal=t_etotal+tcpu()-tt0
2558 if (amax*d_time .gt. dvmax) then
2559 d_time=d_time*dvmax/amax
2560 if(me.eq.king.or..not.out1file) write (iout,*) &
2561 "Time step reduced to",d_time,&
2562 " because of too large initial acceleration."
2564 if(me.eq.king.or..not.out1file)then
2565 write(iout,*) "Potential energy and its components"
2566 call enerprint(potEcomp)
2567 ! write(iout,*) (potEcomp(i),i=0,n_ene)
2569 potE=potEcomp(0)-potEcomp(20)
2572 if (ntwe.ne.0) call statout(itime)
2573 if(me.eq.king.or..not.out1file) &
2574 write (iout,'(/a/3(a25,1pe14.5/))') "Initial:", &
2575 " Kinetic energy",EK," Potential energy",potE, &
2576 " Total energy",totE," Maximum acceleration ", &
2579 write (iout,*) "Initial coordinates"
2581 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),&
2584 write (iout,*) "Initial dC"
2586 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),&
2587 (dc(j,i+nres),j=1,3)
2589 write (iout,*) "Initial velocities"
2590 write (iout,"(13x,' backbone ',23x,' side chain')")
2592 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2593 (d_t(j,i+nres),j=1,3)
2595 write (iout,*) "Initial accelerations"
2597 ! write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2598 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),&
2599 (d_a(j,i+nres),j=1,3)
2605 d_t_old(j,i)=d_t(j,i)
2606 d_a_old(j,i)=d_a(j,i)
2608 ! write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
2617 call etotal_short(energia_short)
2618 if (large) call enerprint(potEcomp)
2621 t_eshort=t_eshort+MPI_Wtime()-tt0
2623 t_eshort=t_eshort+tcpu()-tt0
2628 if(.not.out1file .and. large) then
2629 write (iout,*) "energia_long",energia_long(0),&
2630 " energia_short",energia_short(0),&
2631 " total",energia_long(0)+energia_short(0)
2632 write (iout,*) "Initial fast-force accelerations"
2634 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
2635 (d_a(j,i+nres),j=1,3)
2638 ! 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
2641 d_a_short(j,i)=d_a(j,i)
2650 call etotal_long(energia_long)
2651 if (large) call enerprint(potEcomp)
2654 t_elong=t_elong+MPI_Wtime()-tt0
2656 t_elong=t_elong+tcpu()-tt0
2661 if(.not.out1file .and. large) then
2662 write (iout,*) "energia_long",energia_long(0)
2663 write (iout,*) "Initial slow-force accelerations"
2665 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
2666 (d_a(j,i+nres),j=1,3)
2670 t_enegrad=t_enegrad+MPI_Wtime()-tt0
2672 t_enegrad=t_enegrad+tcpu()-tt0
2676 end subroutine init_MD
2677 !-----------------------------------------------------------------------------
2678 subroutine random_vel
2680 ! implicit real*8 (a-h,o-z)
2682 use random, only:anorm_distr
2684 ! include 'DIMENSIONS'
2685 ! include 'COMMON.CONTROL'
2686 ! include 'COMMON.VAR'
2687 ! include 'COMMON.MD'
2689 ! include 'COMMON.LANGEVIN'
2691 ! include 'COMMON.LANGEVIN.lang0'
2693 ! include 'COMMON.CHAIN'
2694 ! include 'COMMON.DERIV'
2695 ! include 'COMMON.GEO'
2696 ! include 'COMMON.LOCAL'
2697 ! include 'COMMON.INTERACT'
2698 ! include 'COMMON.IOUNITS'
2699 ! include 'COMMON.NAMES'
2700 ! include 'COMMON.TIME1'
2701 real(kind=8) :: xv,sigv,lowb,highb ,Ek1
2703 real(kind=8) ,allocatable, dimension(:) :: DDU1,DDU2,DL2,DL1,xsolv,DML,rs
2704 real(kind=8) :: sumx
2706 real(kind=8) ,allocatable, dimension(:) :: rsold
2707 real (kind=8),allocatable,dimension(:,:) :: matold
2710 integer :: i,j,ii,k,ind,mark,imark,mnum
2711 ! Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
2712 ! First generate velocities in the eigenspace of the G matrix
2713 ! write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
2720 sigv=dsqrt((Rb*t_bath)/geigen(i))
2723 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
2725 write (iout,*) "i",i," ii",ii," geigen",geigen(i),&
2726 " d_t_work_new",d_t_work_new(ii)
2737 Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
2740 write (iout,*) "Ek from eigenvectors",Ek1
2741 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
2745 ! Transform velocities to UNRES coordinate space
2746 allocate (DL1(2*nres))
2747 allocate (DDU1(2*nres))
2748 allocate (DL2(2*nres))
2749 allocate (DDU2(2*nres))
2750 allocate (xsolv(2*nres))
2751 allocate (DML(2*nres))
2752 allocate (rs(2*nres))
2754 allocate (rsold(2*nres))
2755 allocate (matold(2*nres,2*nres))
2757 matold(1,1)=DMorig(1)
2758 matold(1,2)=DU1orig(1)
2759 matold(1,3)=DU2orig(1)
2760 write (*,*) DMorig(1),DU1orig(1),DU2orig(1)
2765 matold(i,j)=DMorig(i)
2766 matold(i,j-1)=DU1orig(i-1)
2768 matold(i,j-2)=DU2orig(i-2)
2776 matold(i,j+1)=DU1orig(i)
2782 matold(i,j+2)=DU2orig(i)
2786 matold(dimen,dimen)=DMorig(dimen)
2787 matold(dimen,dimen-1)=DU1orig(dimen-1)
2788 matold(dimen,dimen-2)=DU2orig(dimen-2)
2789 write(iout,*) "old gmatrix"
2790 call matout(dimen,dimen,2*nres,2*nres,matold)
2794 ! Find the ith eigenvector of the pentadiagonal inertiq matrix
2798 DML(j)=DMorig(j)-geigen(i)
2801 DML(j-1)=DMorig(j)-geigen(i)
2806 DDU1(imark-1)=DU2orig(imark-1)
2807 do j=imark+1,dimen-1
2808 DDU1(j-1)=DU1orig(j)
2816 DDU2(j)=DU2orig(j+1)
2825 write (iout,*) "DL2,DL1,DML,DDU1,DDU2"
2826 write(iout,'(10f10.5)') (DL2(k),k=3,dimen-1)
2827 write(iout,'(10f10.5)') (DL1(k),k=2,dimen-1)
2828 write(iout,'(10f10.5)') (DML(k),k=1,dimen-1)
2829 write(iout,'(10f10.5)') (DDU1(k),k=1,dimen-2)
2830 write(iout,'(10f10.5)') (DDU2(k),k=1,dimen-3)
2833 if (imark.gt.2) rs(imark-2)=-DU2orig(imark-2)
2834 if (imark.gt.1) rs(imark-1)=-DU1orig(imark-1)
2835 if (imark.lt.dimen) rs(imark)=-DU1orig(imark)
2836 if (imark.lt.dimen-1) rs(imark+1)=-DU2orig(imark)
2840 ! write (iout,*) "Vector rs"
2842 ! write (iout,*) j,rs(j)
2845 call FDIAG(dimen-1,DL2,DL1,DML,DDU1,DDU2,rs,xsolv,mark)
2852 sumx=-geigen(i)*xsolv(j)
2854 sumx=sumx+matold(j,k)*xsolv(k)
2857 sumx=sumx+matold(j,k)*xsolv(k-1)
2859 write(iout,'(i5,3f10.5)') j,sumx,rsold(j),sumx-rsold(j)
2862 sumx=-geigen(i)*xsolv(j-1)
2864 sumx=sumx+matold(j,k)*xsolv(k)
2867 sumx=sumx+matold(j,k)*xsolv(k-1)
2869 write(iout,'(i5,3f10.5)') j-1,sumx,rsold(j-1),sumx-rsold(j-1)
2873 "Solution of equations system",i," for eigenvalue",geigen(i)
2875 write(iout,'(i5,f10.5)') j,xsolv(j)
2878 do j=dimen-1,imark,-1
2883 write (iout,*) "Un-normalized eigenvector",i," for eigenvalue",geigen(i)
2885 write(iout,'(i5,f10.5)') j,xsolv(j)
2888 ! Normalize ith eigenvector
2891 sumx=sumx+xsolv(j)**2
2895 xsolv(j)=xsolv(j)/sumx
2898 write (iout,*) "Eigenvector",i," for eigenvalue",geigen(i)
2900 write(iout,'(i5,3f10.5)') j,xsolv(j)
2903 ! All done at this point for eigenvector i, exit loop
2911 write (iout,*) "Unable to find eigenvector",i
2914 ! write (iout,*) "i=",i
2916 ! write (iout,*) "k=",k
2919 ! write(iout,*) "ind",ind," ind1",3*(i-1)+k
2920 d_t_work(ind)=d_t_work(ind) &
2921 +xsolv(j)*d_t_work_new(3*(i-1)+k)
2924 enddo ! i (loop over the eigenvectors)
2927 write (iout,*) "d_t_work"
2929 write (iout,"(i5,f10.5)") i,d_t_work(i)
2934 ! if (itype(i,1).eq.10) then
2936 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
2937 iti=iabs(itype(i,mnum))
2938 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2939 if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
2940 .or.(mnum.eq.5)) then
2947 ! write (iout,*) "k",k," ii+k",ii+k," ii+j+k",ii+j+k,"EK1",Ek1
2948 Ek1=Ek1+0.5d0*mp(mnum)*((d_t_work(ii+j+k)+d_t_work_new(ii+k))/2)**2+&
2949 0.5d0*0.25d0*IP(mnum)*(d_t_work(ii+j+k)-d_t_work_new(ii+k))**2
2952 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2953 .and.(mnum.ne.5)) ii=ii+3
2954 write (iout,*) "i",i," itype",itype(i,mnum)," mass",msc(itype(i,mnum),mnum)
2955 write (iout,*) "ii",ii
2958 write (iout,*) "k",k," ii",ii,"EK1",EK1
2959 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2961 Ek1=Ek1+0.5d0*Isc(iabs(itype(i,mnum),mnum))*(d_t_work(ii)-d_t_work(ii-3))**2
2962 Ek1=Ek1+0.5d0*msc(iabs(itype(i,mnum)),mnum)*d_t_work(ii)**2
2964 write (iout,*) "i",i," ii",ii
2966 write (iout,*) "Ek from d_t_work",Ek1
2967 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
2983 d_t(k,j)=d_t_work(ind)
2987 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2988 .and.(mnum.ne.5)) then
2990 d_t(k,j+nres)=d_t_work(ind)
2996 write (iout,*) "Random velocities in the Calpha,SC space"
2998 write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3)
3001 write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3)
3008 ! if (itype(i,1).eq.10) then
3010 if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
3011 .or.(mnum.eq.5)) then
3013 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
3017 d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
3018 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
3023 write (iout,*)"Random velocities in the virtual-bond-vector space"
3025 write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3)
3028 write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3)
3031 write (iout,*) "Ek from d_t_work",Ek1
3032 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
3040 d_t_work(ind)=d_t_work(ind) &
3041 +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
3043 ! write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
3047 ! Transfer to the d_t vector
3049 d_t(j,0)=d_t_work(j)
3055 d_t(j,i)=d_t_work(ind)
3060 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3061 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3062 .and.(mnum.ne.5)) then
3065 d_t(j,i+nres)=d_t_work(ind)
3071 ! write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",&
3072 ! 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
3076 end subroutine random_vel
3077 !-----------------------------------------------------------------------------
3079 subroutine sd_verlet_p_setup
3080 ! Sets up the parameters of stochastic Verlet algorithm
3081 ! implicit real*8 (a-h,o-z)
3082 ! include 'DIMENSIONS'
3083 use control, only: tcpu
3088 ! include 'COMMON.CONTROL'
3089 ! include 'COMMON.VAR'
3090 ! include 'COMMON.MD'
3092 ! include 'COMMON.LANGEVIN'
3094 ! include 'COMMON.LANGEVIN.lang0'
3096 ! include 'COMMON.CHAIN'
3097 ! include 'COMMON.DERIV'
3098 ! include 'COMMON.GEO'
3099 ! include 'COMMON.LOCAL'
3100 ! include 'COMMON.INTERACT'
3101 ! include 'COMMON.IOUNITS'
3102 ! include 'COMMON.NAMES'
3103 ! include 'COMMON.TIME1'
3104 real(kind=8),dimension(6*nres) :: emgdt !(MAXRES6) maxres6=6*maxres
3105 real(kind=8) :: pterm,vterm,rho,rhoc,vsig
3106 real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
3107 prand_vec,vrand_vec1,vrand_vec2 !(MAXRES6) maxres6=6*maxres
3108 logical :: lprn = .false.
3109 real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
3110 real(kind=8) :: ktm,gdt,egdt,gdt2,gdt3,gdt4,gdt5,gdt6,gdt7,gdt8,&
3112 integer :: i,maxres2
3119 ! AL 8/17/04 Code adapted from tinker
3121 ! Get the frictional and random terms for stochastic dynamics in the
3122 ! eigenspace of mass-scaled UNRES friction matrix
3126 gdt = fricgam(i) * d_time
3128 ! Stochastic dynamics reduces to simple MD for zero friction
3130 if (gdt .le. zero) then
3131 pfric_vec(i) = 1.0d0
3132 vfric_vec(i) = d_time
3133 afric_vec(i) = 0.5d0 * d_time * d_time
3134 prand_vec(i) = 0.0d0
3135 vrand_vec1(i) = 0.0d0
3136 vrand_vec2(i) = 0.0d0
3138 ! Analytical expressions when friction coefficient is large
3141 if (gdt .ge. gdt_radius) then
3144 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
3145 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
3146 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
3147 vterm = 1.0d0 - egdt**2
3148 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
3150 ! Use series expansions when friction coefficient is small
3161 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0 &
3162 - gdt5/120.0d0 + gdt6/720.0d0 &
3163 - gdt7/5040.0d0 + gdt8/40320.0d0 &
3164 - gdt9/362880.0d0) / fricgam(i)**2
3165 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
3166 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
3167 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0 &
3168 + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0 &
3169 + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0 &
3170 + 127.0d0*gdt9/90720.0d0
3171 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0 &
3172 - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0 &
3173 - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0 &
3174 - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
3175 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0 &
3176 - 17.0d0*gdt2/1280.0d0 &
3177 + 17.0d0*gdt3/6144.0d0 &
3178 + 40967.0d0*gdt4/34406400.0d0 &
3179 - 57203.0d0*gdt5/275251200.0d0 &
3180 - 1429487.0d0*gdt6/13212057600.0d0)
3183 ! Compute the scaling factors of random terms for the nonzero friction case
3185 ktm = 0.5d0*d_time/fricgam(i)
3186 psig = dsqrt(ktm*pterm) / fricgam(i)
3187 vsig = dsqrt(ktm*vterm)
3188 rhoc = dsqrt(1.0d0 - rho*rho)
3190 vrand_vec1(i) = vsig * rho
3191 vrand_vec2(i) = vsig * rhoc
3196 "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",&
3199 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),&
3200 afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
3204 ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
3207 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
3208 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
3209 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
3210 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
3211 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
3212 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
3215 t_sdsetup=t_sdsetup+MPI_Wtime()
3217 t_sdsetup=t_sdsetup+tcpu()-tt0
3220 end subroutine sd_verlet_p_setup
3221 !-----------------------------------------------------------------------------
3222 subroutine eigtransf1(n,ndim,ab,d,c)
3226 real(kind=8) :: ab(ndim,ndim,n),c(ndim,n),d(ndim)
3232 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
3237 end subroutine eigtransf1
3238 !-----------------------------------------------------------------------------
3239 subroutine eigtransf(n,ndim,a,b,d,c)
3243 real(kind=8) :: a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
3249 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
3254 end subroutine eigtransf
3255 !-----------------------------------------------------------------------------
3256 subroutine sd_verlet1
3258 ! Applying stochastic velocity Verlet algorithm - step 1 to velocities
3260 ! implicit real*8 (a-h,o-z)
3261 ! include 'DIMENSIONS'
3262 ! include 'COMMON.CONTROL'
3263 ! include 'COMMON.VAR'
3264 ! include 'COMMON.MD'
3266 ! include 'COMMON.LANGEVIN'
3268 ! include 'COMMON.LANGEVIN.lang0'
3270 ! include 'COMMON.CHAIN'
3271 ! include 'COMMON.DERIV'
3272 ! include 'COMMON.GEO'
3273 ! include 'COMMON.LOCAL'
3274 ! include 'COMMON.INTERACT'
3275 ! include 'COMMON.IOUNITS'
3276 ! include 'COMMON.NAMES'
3277 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
3278 !el common /stochcalc/ stochforcvec
3279 logical :: lprn = .false.
3280 real(kind=8) :: ddt1,ddt2
3281 integer :: i,j,ind,inres
3283 ! write (iout,*) "dc_old"
3285 ! write (iout,'(i5,3f10.5,5x,3f10.5)')
3286 ! & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3289 dc_work(j)=dc_old(j,0)
3290 d_t_work(j)=d_t_old(j,0)
3291 d_a_work(j)=d_a_old(j,0)
3296 dc_work(ind+j)=dc_old(j,i)
3297 d_t_work(ind+j)=d_t_old(j,i)
3298 d_a_work(ind+j)=d_a_old(j,i)
3304 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3305 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3306 .and.(mnum.ne.5)) then
3308 dc_work(ind+j)=dc_old(j,i+nres)
3309 d_t_work(ind+j)=d_t_old(j,i+nres)
3310 d_a_work(ind+j)=d_a_old(j,i+nres)
3318 "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",&
3322 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
3323 vfric_mat(i,j),afric_mat(i,j),&
3324 prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3332 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) &
3333 +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3334 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3335 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3337 d_t_work_new(i)=ddt1+0.5d0*ddt2
3338 d_t_work(i)=ddt1+ddt2
3343 d_t(j,0)=d_t_work(j)
3348 dc(j,i)=dc_work(ind+j)
3349 d_t(j,i)=d_t_work(ind+j)
3355 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3356 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3357 .and.(mnum.ne.5)) then
3360 dc(j,inres)=dc_work(ind+j)
3361 d_t(j,inres)=d_t_work(ind+j)
3367 end subroutine sd_verlet1
3368 !-----------------------------------------------------------------------------
3369 subroutine sd_verlet2
3371 ! Calculating the adjusted velocities for accelerations
3373 ! implicit real*8 (a-h,o-z)
3374 ! include 'DIMENSIONS'
3375 ! include 'COMMON.CONTROL'
3376 ! include 'COMMON.VAR'
3377 ! include 'COMMON.MD'
3379 ! include 'COMMON.LANGEVIN'
3381 ! include 'COMMON.LANGEVIN.lang0'
3383 ! include 'COMMON.CHAIN'
3384 ! include 'COMMON.DERIV'
3385 ! include 'COMMON.GEO'
3386 ! include 'COMMON.LOCAL'
3387 ! include 'COMMON.INTERACT'
3388 ! include 'COMMON.IOUNITS'
3389 ! include 'COMMON.NAMES'
3390 !el real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV !(MAXRES6) maxres6=6*maxres
3391 real(kind=8),dimension(6*nres) :: stochforcvecV !(MAXRES6) maxres6=6*maxres
3392 !el common /stochcalc/ stochforcvec
3394 real(kind=8) :: ddt1,ddt2
3395 integer :: i,j,ind,inres
3396 ! Compute the stochastic forces which contribute to velocity change
3398 call stochastic_force(stochforcvecV)
3405 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3406 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+ &
3407 vrand_mat2(i,j)*stochforcvecV(j)
3409 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
3413 d_t(j,0)=d_t_work(j)
3418 d_t(j,i)=d_t_work(ind+j)
3424 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3425 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3426 .and.(mnum.ne.5)) then
3429 d_t(j,inres)=d_t_work(ind+j)
3435 end subroutine sd_verlet2
3436 !-----------------------------------------------------------------------------
3437 subroutine sd_verlet_ciccotti_setup
3439 ! Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
3441 ! implicit real*8 (a-h,o-z)
3442 ! include 'DIMENSIONS'
3443 use control, only: tcpu
3448 ! include 'COMMON.CONTROL'
3449 ! include 'COMMON.VAR'
3450 ! include 'COMMON.MD'
3452 ! include 'COMMON.LANGEVIN'
3454 ! include 'COMMON.LANGEVIN.lang0'
3456 ! include 'COMMON.CHAIN'
3457 ! include 'COMMON.DERIV'
3458 ! include 'COMMON.GEO'
3459 ! include 'COMMON.LOCAL'
3460 ! include 'COMMON.INTERACT'
3461 ! include 'COMMON.IOUNITS'
3462 ! include 'COMMON.NAMES'
3463 ! include 'COMMON.TIME1'
3464 real(kind=8),dimension(6*nres) :: emgdt !(MAXRES6) maxres6=6*maxres
3465 real(kind=8) :: pterm,vterm,rho,rhoc,vsig
3466 real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
3467 prand_vec,vrand_vec1,vrand_vec2 !(MAXRES6) maxres6=6*maxres
3468 logical :: lprn = .false.
3469 real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
3470 real(kind=8) :: ktm,gdt,egdt,tt0
3471 integer :: i,maxres2
3478 ! AL 8/17/04 Code adapted from tinker
3480 ! Get the frictional and random terms for stochastic dynamics in the
3481 ! eigenspace of mass-scaled UNRES friction matrix
3485 write (iout,*) "i",i," fricgam",fricgam(i)
3486 gdt = fricgam(i) * d_time
3488 ! Stochastic dynamics reduces to simple MD for zero friction
3490 if (gdt .le. zero) then
3491 pfric_vec(i) = 1.0d0
3492 vfric_vec(i) = d_time
3493 afric_vec(i) = 0.5d0*d_time*d_time
3494 prand_vec(i) = afric_vec(i)
3495 vrand_vec2(i) = vfric_vec(i)
3497 ! Analytical expressions when friction coefficient is large
3502 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
3503 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
3504 prand_vec(i) = afric_vec(i)
3505 vrand_vec2(i) = vfric_vec(i)
3507 ! Compute the scaling factors of random terms for the nonzero friction case
3509 ! ktm = 0.5d0*d_time/fricgam(i)
3510 ! psig = dsqrt(ktm*pterm) / fricgam(i)
3511 ! vsig = dsqrt(ktm*vterm)
3512 ! prand_vec(i) = psig*afric_vec(i)
3513 ! vrand_vec2(i) = vsig*vfric_vec(i)
3518 "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",&
3521 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),&
3522 afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
3526 ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
3528 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
3529 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
3530 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
3531 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
3532 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
3534 t_sdsetup=t_sdsetup+MPI_Wtime()
3536 t_sdsetup=t_sdsetup+tcpu()-tt0
3539 end subroutine sd_verlet_ciccotti_setup
3540 !-----------------------------------------------------------------------------
3541 subroutine sd_verlet1_ciccotti
3543 ! Applying stochastic velocity Verlet algorithm - step 1 to velocities
3544 ! implicit real*8 (a-h,o-z)
3546 ! include 'DIMENSIONS'
3550 ! include 'COMMON.CONTROL'
3551 ! include 'COMMON.VAR'
3552 ! include 'COMMON.MD'
3554 ! include 'COMMON.LANGEVIN'
3556 ! include 'COMMON.LANGEVIN.lang0'
3558 ! include 'COMMON.CHAIN'
3559 ! include 'COMMON.DERIV'
3560 ! include 'COMMON.GEO'
3561 ! include 'COMMON.LOCAL'
3562 ! include 'COMMON.INTERACT'
3563 ! include 'COMMON.IOUNITS'
3564 ! include 'COMMON.NAMES'
3565 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
3566 !el common /stochcalc/ stochforcvec
3567 logical :: lprn = .false.
3568 real(kind=8) :: ddt1,ddt2
3569 integer :: i,j,ind,inres
3570 ! write (iout,*) "dc_old"
3572 ! write (iout,'(i5,3f10.5,5x,3f10.5)')
3573 ! & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3576 dc_work(j)=dc_old(j,0)
3577 d_t_work(j)=d_t_old(j,0)
3578 d_a_work(j)=d_a_old(j,0)
3583 dc_work(ind+j)=dc_old(j,i)
3584 d_t_work(ind+j)=d_t_old(j,i)
3585 d_a_work(ind+j)=d_a_old(j,i)
3590 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3592 dc_work(ind+j)=dc_old(j,i+nres)
3593 d_t_work(ind+j)=d_t_old(j,i+nres)
3594 d_a_work(ind+j)=d_a_old(j,i+nres)
3603 "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",&
3607 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
3608 vfric_mat(i,j),afric_mat(i,j),&
3609 prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3617 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) &
3618 +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3619 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3620 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3622 d_t_work_new(i)=ddt1+0.5d0*ddt2
3623 d_t_work(i)=ddt1+ddt2
3628 d_t(j,0)=d_t_work(j)
3633 dc(j,i)=dc_work(ind+j)
3634 d_t(j,i)=d_t_work(ind+j)
3640 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3641 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3642 .and.(mnum.ne.5)) then
3645 dc(j,inres)=dc_work(ind+j)
3646 d_t(j,inres)=d_t_work(ind+j)
3652 end subroutine sd_verlet1_ciccotti
3653 !-----------------------------------------------------------------------------
3654 subroutine sd_verlet2_ciccotti
3656 ! Calculating the adjusted velocities for accelerations
3658 ! implicit real*8 (a-h,o-z)
3659 ! include 'DIMENSIONS'
3660 ! include 'COMMON.CONTROL'
3661 ! include 'COMMON.VAR'
3662 ! include 'COMMON.MD'
3664 ! include 'COMMON.LANGEVIN'
3666 ! include 'COMMON.LANGEVIN.lang0'
3668 ! include 'COMMON.CHAIN'
3669 ! include 'COMMON.DERIV'
3670 ! include 'COMMON.GEO'
3671 ! include 'COMMON.LOCAL'
3672 ! include 'COMMON.INTERACT'
3673 ! include 'COMMON.IOUNITS'
3674 ! include 'COMMON.NAMES'
3675 !el real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV !(MAXRES6) maxres6=6*maxres
3676 real(kind=8),dimension(6*nres) :: stochforcvecV !(MAXRES6) maxres6=6*maxres
3677 !el common /stochcalc/ stochforcvec
3678 real(kind=8) :: ddt1,ddt2
3679 integer :: i,j,ind,inres
3681 ! Compute the stochastic forces which contribute to velocity change
3683 call stochastic_force(stochforcvecV)
3690 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3691 ! ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
3692 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
3694 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
3698 d_t(j,0)=d_t_work(j)
3703 d_t(j,i)=d_t_work(ind+j)
3709 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3711 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3714 d_t(j,inres)=d_t_work(ind+j)
3720 end subroutine sd_verlet2_ciccotti
3722 !-----------------------------------------------------------------------------
3724 !-----------------------------------------------------------------------------
3725 subroutine inertia_tensor
3727 ! Calculating the intertia tensor for the entire protein in order to
3728 ! remove the perpendicular components of velocity matrix which cause
3729 ! the molecule to rotate.
3732 ! implicit real*8 (a-h,o-z)
3733 ! include 'DIMENSIONS'
3734 ! include 'COMMON.CONTROL'
3735 ! include 'COMMON.VAR'
3736 ! include 'COMMON.MD'
3737 ! include 'COMMON.CHAIN'
3738 ! include 'COMMON.DERIV'
3739 ! include 'COMMON.GEO'
3740 ! include 'COMMON.LOCAL'
3741 ! include 'COMMON.INTERACT'
3742 ! include 'COMMON.IOUNITS'
3743 ! include 'COMMON.NAMES'
3745 real(kind=8),dimension(3,3) :: Im,Imcp,eigvec,Id
3746 real(kind=8),dimension(3) :: pr,eigval,L,vp,vrot
3747 real(kind=8) :: M_SC,mag,mag2,M_PEP
3748 real(kind=8),dimension(3,0:nres) :: vpp !(3,0:MAXRES)
3749 real(kind=8),dimension(3) :: vs_p,pp,incr,v
3750 real(kind=8),dimension(3,3) :: pr1,pr2
3752 !el common /gucio/ cm
3753 integer :: iti,inres,i,j,k,mnum
3764 ! calculating the center of the mass of the protein
3768 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
3769 if (itype(i,mnum).eq.ntyp1_molec(mnum)) cycle
3770 M_PEP=M_PEP+mp(mnum)
3772 cm(j)=cm(j)+(c(j,i)+0.5d0*dc(j,i))*mp(mnum)
3781 if (mnum.eq.5) cycle
3782 iti=iabs(itype(i,mnum))
3783 M_SC=M_SC+msc(iabs(iti),mnum)
3786 cm(j)=cm(j)+msc(iabs(iti),mnum)*c(j,inres)
3790 cm(j)=cm(j)/(M_SC+M_PEP)
3795 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
3797 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
3799 Im(1,1)=Im(1,1)+mp(mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
3800 Im(1,2)=Im(1,2)-mp(mnum)*pr(1)*pr(2)
3801 Im(1,3)=Im(1,3)-mp(mnum)*pr(1)*pr(3)
3802 Im(2,3)=Im(2,3)-mp(mnum)*pr(2)*pr(3)
3803 Im(2,2)=Im(2,2)+mp(mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
3804 Im(3,3)=Im(3,3)+mp(mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
3809 iti=iabs(itype(i,mnum))
3810 if (mnum.eq.5) cycle
3813 pr(j)=c(j,inres)-cm(j)
3815 Im(1,1)=Im(1,1)+msc(iabs(iti),mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
3816 Im(1,2)=Im(1,2)-msc(iabs(iti),mnum)*pr(1)*pr(2)
3817 Im(1,3)=Im(1,3)-msc(iabs(iti),mnum)*pr(1)*pr(3)
3818 Im(2,3)=Im(2,3)-msc(iabs(iti),mnum)*pr(2)*pr(3)
3819 Im(2,2)=Im(2,2)+msc(iabs(iti),mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
3820 Im(3,3)=Im(3,3)+msc(iabs(iti),mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
3825 Im(1,1)=Im(1,1)+Ip(mnum)*(1-dc_norm(1,i)*dc_norm(1,i))* &
3826 vbld(i+1)*vbld(i+1)*0.25d0
3827 Im(1,2)=Im(1,2)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(2,i))* &
3828 vbld(i+1)*vbld(i+1)*0.25d0
3829 Im(1,3)=Im(1,3)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(3,i))* &
3830 vbld(i+1)*vbld(i+1)*0.25d0
3831 Im(2,3)=Im(2,3)+Ip(mnum)*(-dc_norm(2,i)*dc_norm(3,i))* &
3832 vbld(i+1)*vbld(i+1)*0.25d0
3833 Im(2,2)=Im(2,2)+Ip(mnum)*(1-dc_norm(2,i)*dc_norm(2,i))* &
3834 vbld(i+1)*vbld(i+1)*0.25d0
3835 Im(3,3)=Im(3,3)+Ip(mnum)*(1-dc_norm(3,i)*dc_norm(3,i))* &
3836 vbld(i+1)*vbld(i+1)*0.25d0
3842 ! if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
3843 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3844 .and.(mnum.ne.5)) then
3845 iti=iabs(itype(i,mnum))
3847 Im(1,1)=Im(1,1)+Isc(iti,mnum)*(1-dc_norm(1,inres)* &
3848 dc_norm(1,inres))*vbld(inres)*vbld(inres)
3849 Im(1,2)=Im(1,2)-Isc(iti,mnum)*(dc_norm(1,inres)* &
3850 dc_norm(2,inres))*vbld(inres)*vbld(inres)
3851 Im(1,3)=Im(1,3)-Isc(iti,mnum)*(dc_norm(1,inres)* &
3852 dc_norm(3,inres))*vbld(inres)*vbld(inres)
3853 Im(2,3)=Im(2,3)-Isc(iti,mnum)*(dc_norm(2,inres)* &
3854 dc_norm(3,inres))*vbld(inres)*vbld(inres)
3855 Im(2,2)=Im(2,2)+Isc(iti,mnum)*(1-dc_norm(2,inres)* &
3856 dc_norm(2,inres))*vbld(inres)*vbld(inres)
3857 Im(3,3)=Im(3,3)+Isc(iti,mnum)*(1-dc_norm(3,inres)* &
3858 dc_norm(3,inres))*vbld(inres)*vbld(inres)
3863 ! write(iout,*) "The angular momentum before adjustment:"
3864 ! write(iout,*) (L(j),j=1,3)
3870 ! Copying the Im matrix for the djacob subroutine
3878 ! Finding the eigenvectors and eignvalues of the inertia tensor
3879 call djacob(3,3,10000,1.0d-10,Imcp,eigvec,eigval)
3880 ! write (iout,*) "Eigenvalues & Eigenvectors"
3881 ! write (iout,'(5x,3f10.5)') (eigval(i),i=1,3)
3884 ! write (iout,'(i5,3f10.5)') i,(eigvec(i,j),j=1,3)
3886 ! Constructing the diagonalized matrix
3888 if (dabs(eigval(i)).gt.1.0d-15) then
3889 Id(i,i)=1.0d0/eigval(i)
3896 Imcp(i,j)=eigvec(j,i)
3902 pr1(i,j)=pr1(i,j)+Id(i,k)*Imcp(k,j)
3909 pr2(i,j)=pr2(i,j)+eigvec(i,k)*pr1(k,j)
3913 ! Calculating the total rotational velocity of the molecule
3916 vrot(i)=vrot(i)+pr2(i,j)*L(j)
3919 ! Resetting the velocities
3921 call vecpr(vrot(1),dc(1,i),vp)
3923 d_t(j,i)=d_t(j,i)-vp(j)
3928 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3929 .and.(mnum.ne.5)) then
3930 ! if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
3931 ! if(itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3933 call vecpr(vrot(1),dc(1,inres),vp)
3935 d_t(j,inres)=d_t(j,inres)-vp(j)
3940 ! write(iout,*) "The angular momentum after adjustment:"
3941 ! write(iout,*) (L(j),j=1,3)
3944 end subroutine inertia_tensor
3945 !-----------------------------------------------------------------------------
3946 subroutine angmom(cm,L)
3949 ! implicit real*8 (a-h,o-z)
3950 ! include 'DIMENSIONS'
3951 ! include 'COMMON.CONTROL'
3952 ! include 'COMMON.VAR'
3953 ! include 'COMMON.MD'
3954 ! include 'COMMON.CHAIN'
3955 ! include 'COMMON.DERIV'
3956 ! include 'COMMON.GEO'
3957 ! include 'COMMON.LOCAL'
3958 ! include 'COMMON.INTERACT'
3959 ! include 'COMMON.IOUNITS'
3960 ! include 'COMMON.NAMES'
3961 real(kind=8) :: mscab
3962 real(kind=8),dimension(3) :: L,cm,pr,vp,vrot,incr,v,pp
3963 integer :: iti,inres,i,j,mnum
3964 ! Calculate the angular momentum
3973 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
3975 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
3978 v(j)=incr(j)+0.5d0*d_t(j,i)
3981 incr(j)=incr(j)+d_t(j,i)
3983 call vecpr(pr(1),v(1),vp)
3985 L(j)=L(j)+mp(mnum)*vp(j)
3989 pp(j)=0.5d0*d_t(j,i)
3991 call vecpr(pr(1),pp(1),vp)
3993 L(j)=L(j)+Ip(mnum)*vp(j)
4001 iti=iabs(itype(i,mnum))
4009 pr(j)=c(j,inres)-cm(j)
4011 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4012 .and.(mnum.ne.5)) then
4014 v(j)=incr(j)+d_t(j,inres)
4021 call vecpr(pr(1),v(1),vp)
4022 ! write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),
4023 ! & " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
4025 L(j)=L(j)+mscab*vp(j)
4027 ! write (iout,*) "L",(l(j),j=1,3)
4028 if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4029 .and.(mnum.ne.5)) then
4031 v(j)=incr(j)+d_t(j,inres)
4033 call vecpr(dc(1,inres),d_t(1,inres),vp)
4035 L(j)=L(j)+Isc(iti,mnum)*vp(j)
4039 incr(j)=incr(j)+d_t(j,i)
4043 end subroutine angmom
4044 !-----------------------------------------------------------------------------
4045 subroutine vcm_vel(vcm)
4048 ! implicit real*8 (a-h,o-z)
4049 ! include 'DIMENSIONS'
4050 ! include 'COMMON.VAR'
4051 ! include 'COMMON.MD'
4052 ! include 'COMMON.CHAIN'
4053 ! include 'COMMON.DERIV'
4054 ! include 'COMMON.GEO'
4055 ! include 'COMMON.LOCAL'
4056 ! include 'COMMON.INTERACT'
4057 ! include 'COMMON.IOUNITS'
4058 real(kind=8),dimension(3) :: vcm,vv
4059 real(kind=8) :: summas,amas
4069 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
4071 summas=summas+mp(mnum)
4073 vcm(j)=vcm(j)+mp(mnum)*(vv(j)+0.5d0*d_t(j,i))
4077 amas=msc(iabs(itype(i,mnum)),mnum)
4082 if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4083 .and.(mnum.ne.5)) then
4085 vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres))
4089 vcm(j)=vcm(j)+amas*vv(j)
4093 vv(j)=vv(j)+d_t(j,i)
4096 ! write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas
4098 vcm(j)=vcm(j)/summas
4101 end subroutine vcm_vel
4102 !-----------------------------------------------------------------------------
4104 !-----------------------------------------------------------------------------
4106 ! RATTLE algorithm for velocity Verlet - step 1, UNRES
4110 ! implicit real*8 (a-h,o-z)
4111 ! include 'DIMENSIONS'
4113 ! include 'COMMON.CONTROL'
4114 ! include 'COMMON.VAR'
4115 ! include 'COMMON.MD'
4117 ! include 'COMMON.LANGEVIN'
4119 ! include 'COMMON.LANGEVIN.lang0'
4121 ! include 'COMMON.CHAIN'
4122 ! include 'COMMON.DERIV'
4123 ! include 'COMMON.GEO'
4124 ! include 'COMMON.LOCAL'
4125 ! include 'COMMON.INTERACT'
4126 ! include 'COMMON.IOUNITS'
4127 ! include 'COMMON.NAMES'
4128 ! include 'COMMON.TIME1'
4129 !el real(kind=8) :: gginv(2*nres,2*nres),&
4130 !el gdc(3,2*nres,2*nres)
4131 real(kind=8) :: dC_uncor(3,2*nres) !,&
4132 !el real(kind=8) :: Cmat(2*nres,2*nres)
4133 real(kind=8) :: x(2*nres),xcorr(3,2*nres) !maxres2=2*maxres
4134 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
4135 !el common /przechowalnia/ nbond
4136 integer :: max_rattle = 5
4137 logical :: lprn = .false., lprn1 = .false., not_done
4138 real(kind=8) :: tol_rattle = 1.0d-5
4140 integer :: ii,i,j,jj,l,ind,ind1,nres2
4143 !el /common/ przechowalnia
4145 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4146 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4147 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4149 if (lprn) write (iout,*) "RATTLE1"
4153 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4154 .and.(mnum.ne.5)) nbond=nbond+1
4156 ! Make a folded form of the Ginv-matrix
4169 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4174 if (itype(k,1).ne.10 .and. itype(k,mnum).ne.ntyp1_molec(mnum)&
4175 .and.(mnum.ne.5)) then
4179 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4187 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4198 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4202 if (itype(k,1).ne.10) then
4206 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4214 write (iout,*) "Matrix GGinv"
4215 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
4221 if (iter.gt.max_rattle) then
4222 write (iout,*) "Error - too many iterations in RATTLE."
4225 ! Calculate the matrix C = GG**(-1) dC_old o dC
4230 dC_uncor(j,ind1)=dC(j,i)
4234 if (itype(i,1).ne.10) then
4237 dC_uncor(j,ind1)=dC(j,i+nres)
4246 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
4250 if (itype(k,1).ne.10) then
4253 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
4258 ! Calculate deviations from standard virtual-bond lengths
4262 x(ind)=vbld(i+1)**2-vbl**2
4265 if (itype(i,1).ne.10) then
4267 x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4271 write (iout,*) "Coordinates and violations"
4273 write(iout,'(i5,3f10.5,5x,e15.5)') &
4274 i,(dC_uncor(j,i),j=1,3),x(i)
4276 write (iout,*) "Velocities and violations"
4280 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4281 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4285 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4286 .and.(mnum.ne.5)) then
4289 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4290 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4291 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4294 ! write (iout,*) "gdc"
4296 ! write (iout,*) "i",i
4298 ! write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4304 if (dabs(x(i)).gt.xmax) then
4308 if (xmax.lt.tol_rattle) then
4312 ! Calculate the matrix of the system of equations
4317 Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
4322 write (iout,*) "Matrix Cmat"
4323 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4325 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4326 ! Add constraint term to positions
4333 xx = xx+x(ii)*gdc(j,ind,ii)
4337 d_t_new(j,i)=d_t_new(j,i)-xx/d_time
4341 if (itype(i,1).ne.10) then
4346 xx = xx+x(ii)*gdc(j,ind,ii)
4349 dC(j,i+nres)=dC(j,i+nres)-xx
4350 d_t_new(j,i+nres)=d_t_new(j,i+nres)-xx/d_time
4354 ! Rebuild the chain using the new coordinates
4355 call chainbuild_cart
4357 write (iout,*) "New coordinates, Lagrange multipliers,",&
4358 " and differences between actual and standard bond lengths"
4362 xx=vbld(i+1)**2-vbl**2
4363 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4364 i,(dC(j,i),j=1,3),x(ind),xx
4368 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4371 xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4372 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4373 i,(dC(j,i+nres),j=1,3),x(ind),xx
4376 write (iout,*) "Velocities and violations"
4380 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4381 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4384 if (itype(i,1).ne.10) then
4386 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4387 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4388 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4395 10 write (iout,*) "Error - singularity in solving the system",&
4396 " of equations for Lagrange multipliers."
4400 "RATTLE inactive; use -DRATTLE switch at compile time."
4403 end subroutine rattle1
4404 !-----------------------------------------------------------------------------
4406 ! RATTLE algorithm for velocity Verlet - step 2, UNRES
4410 ! implicit real*8 (a-h,o-z)
4411 ! include 'DIMENSIONS'
4413 ! include 'COMMON.CONTROL'
4414 ! include 'COMMON.VAR'
4415 ! include 'COMMON.MD'
4417 ! include 'COMMON.LANGEVIN'
4419 ! include 'COMMON.LANGEVIN.lang0'
4421 ! include 'COMMON.CHAIN'
4422 ! include 'COMMON.DERIV'
4423 ! include 'COMMON.GEO'
4424 ! include 'COMMON.LOCAL'
4425 ! include 'COMMON.INTERACT'
4426 ! include 'COMMON.IOUNITS'
4427 ! include 'COMMON.NAMES'
4428 ! include 'COMMON.TIME1'
4429 !el real(kind=8) :: gginv(2*nres,2*nres),&
4430 !el gdc(3,2*nres,2*nres)
4431 real(kind=8) :: dC_uncor(3,2*nres) !,&
4432 !el Cmat(2*nres,2*nres)
4433 real(kind=8) :: x(2*nres) !maxres2=2*maxres
4434 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
4435 !el common /przechowalnia/ nbond
4436 integer :: max_rattle = 5
4437 logical :: lprn = .false., lprn1 = .false., not_done
4438 real(kind=8) :: tol_rattle = 1.0d-5
4442 !el /common/ przechowalnia
4443 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4444 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4445 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4447 if (lprn) write (iout,*) "RATTLE2"
4448 if (lprn) write (iout,*) "Velocity correction"
4449 ! Calculate the matrix G dC
4455 gdc(j,i,ind)=GGinv(i,ind)*dC(j,k)
4460 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4461 .and.(mnum.ne.5)) then
4464 gdc(j,i,ind)=GGinv(i,ind)*dC(j,k+nres)
4470 ! write (iout,*) "gdc"
4472 ! write (iout,*) "i",i
4474 ! write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4478 ! Calculate the matrix of the system of equations
4485 Cmat(ind,j)=Cmat(ind,j)+dC(k,i)*gdc(k,ind,j)
4491 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4492 .and.(mnum.ne.5)) then
4497 Cmat(ind,j)=Cmat(ind,j)+dC(k,i+nres)*gdc(k,ind,j)
4502 ! Calculate the scalar product dC o d_t_new
4506 x(ind)=scalar(d_t(1,i),dC(1,i))
4510 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4511 .and.(mnum.ne.5)) then
4513 x(ind)=scalar(d_t(1,i+nres),dC(1,i+nres))
4517 write (iout,*) "Velocities and violations"
4521 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4522 i,ind,(d_t(j,i),j=1,3),x(ind)
4526 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4527 .and.(mnum.ne.5)) then
4529 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4530 i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind)
4536 if (dabs(x(i)).gt.xmax) then
4540 if (xmax.lt.tol_rattle) then
4545 write (iout,*) "Matrix Cmat"
4546 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4548 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4549 ! Add constraint term to velocities
4556 xx = xx+x(ii)*gdc(j,ind,ii)
4558 d_t(j,i)=d_t(j,i)-xx
4563 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4564 .and.(mnum.ne.5)) then
4569 xx = xx+x(ii)*gdc(j,ind,ii)
4571 d_t(j,i+nres)=d_t(j,i+nres)-xx
4577 "New velocities, Lagrange multipliers violations"
4581 if (lprn) write (iout,'(2i5,3f10.5,5x,2e15.5)') &
4582 i,ind,(d_t(j,i),j=1,3),x(ind),scalar(d_t(1,i),dC(1,i))
4586 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4589 write (iout,'(2i5,3f10.5,5x,2e15.5)') &
4590 i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind),&
4591 scalar(d_t(1,i+nres),dC(1,i+nres))
4597 10 write (iout,*) "Error - singularity in solving the system",&
4598 " of equations for Lagrange multipliers."
4602 "RATTLE inactive; use -DRATTLE option at compile time."
4605 end subroutine rattle2
4606 !-----------------------------------------------------------------------------
4607 subroutine rattle_brown
4608 ! RATTLE/LINCS algorithm for Brownian dynamics, UNRES
4612 ! implicit real*8 (a-h,o-z)
4613 ! include 'DIMENSIONS'
4615 ! include 'COMMON.CONTROL'
4616 ! include 'COMMON.VAR'
4617 ! include 'COMMON.MD'
4619 ! include 'COMMON.LANGEVIN'
4621 ! include 'COMMON.LANGEVIN.lang0'
4623 ! include 'COMMON.CHAIN'
4624 ! include 'COMMON.DERIV'
4625 ! include 'COMMON.GEO'
4626 ! include 'COMMON.LOCAL'
4627 ! include 'COMMON.INTERACT'
4628 ! include 'COMMON.IOUNITS'
4629 ! include 'COMMON.NAMES'
4630 ! include 'COMMON.TIME1'
4631 !el real(kind=8) :: gginv(2*nres,2*nres),&
4632 !el gdc(3,2*nres,2*nres)
4633 real(kind=8) :: dC_uncor(3,2*nres) !,&
4634 !el real(kind=8) :: Cmat(2*nres,2*nres)
4635 real(kind=8) :: x(2*nres) !maxres2=2*maxres
4636 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
4637 !el common /przechowalnia/ nbond
4638 integer :: max_rattle = 5
4639 logical :: lprn = .false., lprn1 = .false., not_done
4640 real(kind=8) :: tol_rattle = 1.0d-5
4644 !el /common/ przechowalnia
4645 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4646 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4647 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4650 if (lprn) write (iout,*) "RATTLE_BROWN"
4653 if (itype(i,1).ne.10) nbond=nbond+1
4655 ! Make a folded form of the Ginv-matrix
4668 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4672 if (itype(k,1).ne.10) then
4676 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4683 if (itype(i,1).ne.10) then
4693 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4697 if (itype(k,1).ne.10) then
4701 if (j.eq.1 .and. l.eq.1)GGinv(ii,jj)=fricmat(ind,ind1)
4709 write (iout,*) "Matrix GGinv"
4710 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
4716 if (iter.gt.max_rattle) then
4717 write (iout,*) "Error - too many iterations in RATTLE."
4720 ! Calculate the matrix C = GG**(-1) dC_old o dC
4725 dC_uncor(j,ind1)=dC(j,i)
4729 if (itype(i,1).ne.10) then
4732 dC_uncor(j,ind1)=dC(j,i+nres)
4741 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
4745 if (itype(k,1).ne.10) then
4748 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
4753 ! Calculate deviations from standard virtual-bond lengths
4757 x(ind)=vbld(i+1)**2-vbl**2
4760 if (itype(i,1).ne.10) then
4762 x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4766 write (iout,*) "Coordinates and violations"
4768 write(iout,'(i5,3f10.5,5x,e15.5)') &
4769 i,(dC_uncor(j,i),j=1,3),x(i)
4771 write (iout,*) "Velocities and violations"
4775 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4776 i,ind,(d_t(j,i),j=1,3),scalar(d_t(1,i),dC_old(1,i))
4779 if (itype(i,1).ne.10) then
4781 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4782 i+nres,ind,(d_t(j,i+nres),j=1,3),&
4783 scalar(d_t(1,i+nres),dC_old(1,i+nres))
4786 write (iout,*) "gdc"
4788 write (iout,*) "i",i
4790 write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4796 if (dabs(x(i)).gt.xmax) then
4800 if (xmax.lt.tol_rattle) then
4804 ! Calculate the matrix of the system of equations
4809 Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
4814 write (iout,*) "Matrix Cmat"
4815 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4817 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4818 ! Add constraint term to positions
4825 xx = xx+x(ii)*gdc(j,ind,ii)
4828 d_t(j,i)=d_t(j,i)+xx/d_time
4833 if (itype(i,1).ne.10) then
4838 xx = xx+x(ii)*gdc(j,ind,ii)
4841 d_t(j,i+nres)=d_t(j,i+nres)+xx/d_time
4842 dC(j,i+nres)=dC(j,i+nres)+xx
4846 ! Rebuild the chain using the new coordinates
4847 call chainbuild_cart
4849 write (iout,*) "New coordinates, Lagrange multipliers,",&
4850 " and differences between actual and standard bond lengths"
4854 xx=vbld(i+1)**2-vbl**2
4855 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4856 i,(dC(j,i),j=1,3),x(ind),xx
4859 if (itype(i,1).ne.10) then
4861 xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4862 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4863 i,(dC(j,i+nres),j=1,3),x(ind),xx
4866 write (iout,*) "Velocities and violations"
4870 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4871 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4874 if (itype(i,1).ne.10) then
4876 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4877 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4878 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4885 10 write (iout,*) "Error - singularity in solving the system",&
4886 " of equations for Lagrange multipliers."
4890 "RATTLE inactive; use -DRATTLE option at compile time"
4893 end subroutine rattle_brown
4894 !-----------------------------------------------------------------------------
4896 !-----------------------------------------------------------------------------
4897 subroutine friction_force
4902 ! implicit real*8 (a-h,o-z)
4903 ! include 'DIMENSIONS'
4904 ! include 'COMMON.VAR'
4905 ! include 'COMMON.CHAIN'
4906 ! include 'COMMON.DERIV'
4907 ! include 'COMMON.GEO'
4908 ! include 'COMMON.LOCAL'
4909 ! include 'COMMON.INTERACT'
4910 ! include 'COMMON.MD'
4912 ! include 'COMMON.LANGEVIN'
4914 ! include 'COMMON.LANGEVIN.lang0'
4916 ! include 'COMMON.IOUNITS'
4917 !el real(kind=8),dimension(6*nres) :: gamvec !(MAXRES6) maxres6=6*maxres
4918 !el common /syfek/ gamvec
4919 real(kind=8) :: vv(3),vvtot(3,nres),v_work(6*nres) !,&
4920 !el ginvfric(2*nres,2*nres) !maxres2=2*maxres
4921 !el common /przechowalnia/ ginvfric
4923 logical :: lprn = .false., checkmode = .false.
4924 integer :: i,j,ind,k,nres2,nres6,mnum
4928 if(.not.allocated(gamvec)) allocate(gamvec(nres6)) !(MAXRES6)
4929 if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
4937 d_t_work(j)=d_t(j,0)
4942 d_t_work(ind+j)=d_t(j,i)
4948 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
4949 .and.(mnum.ne.5)) then
4951 d_t_work(ind+j)=d_t(j,i+nres)
4957 call fricmat_mult(d_t_work,fric_work)
4959 if (.not.checkmode) return
4962 write (iout,*) "d_t_work and fric_work"
4964 write (iout,'(i3,2e15.5)') i,d_t_work(i),fric_work(i)
4968 friction(j,0)=fric_work(j)
4973 friction(j,i)=fric_work(ind+j)
4979 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
4980 .and.(mnum.ne.5)) then
4981 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
4983 friction(j,i+nres)=fric_work(ind+j)
4989 write(iout,*) "Friction backbone"
4991 write(iout,'(i5,3e15.5,5x,3e15.5)') &
4992 i,(friction(j,i),j=1,3),(d_t(j,i),j=1,3)
4994 write(iout,*) "Friction side chain"
4996 write(iout,'(i5,3e15.5,5x,3e15.5)') &
4997 i,(friction(j,i+nres),j=1,3),(d_t(j,i+nres),j=1,3)
5006 vvtot(j,i)=vv(j)+0.5d0*d_t(j,i)
5007 vvtot(j,i+nres)=vv(j)+d_t(j,i+nres)
5008 vv(j)=vv(j)+d_t(j,i)
5011 write (iout,*) "vvtot backbone and sidechain"
5013 write (iout,'(i5,3e15.5,5x,3e15.5)') i,(vvtot(j,i),j=1,3),&
5014 (vvtot(j,i+nres),j=1,3)
5019 v_work(ind+j)=vvtot(j,i)
5025 v_work(ind+j)=vvtot(j,i+nres)
5029 write (iout,*) "v_work gamvec and site-based friction forces"
5031 write (iout,'(i5,3e15.5)') i,v_work(i),gamvec(i),&
5035 ! fric_work1(i)=0.0d0
5037 ! fric_work1(i)=fric_work1(i)-A(j,i)*gamvec(j)*v_work(j)
5040 ! write (iout,*) "fric_work and fric_work1"
5042 ! write (iout,'(i5,2e15.5)') i,fric_work(i),fric_work1(i)
5048 ginvfric(i,j)=ginvfric(i,j)+ginv(i,k)*fricmat(k,j)
5052 write (iout,*) "ginvfric"
5054 write (iout,'(i5,100f8.3)') i,(ginvfric(i,j),j=1,dimen)
5056 write (iout,*) "symmetry check"
5059 write (iout,*) i,j,ginvfric(i,j)-ginvfric(j,i)
5064 end subroutine friction_force
5065 !-----------------------------------------------------------------------------
5066 subroutine setup_fricmat
5070 use control_data, only:time_Bcast
5071 use control, only:tcpu
5073 ! implicit real*8 (a-h,o-z)
5077 real(kind=8) :: time00
5079 ! include 'DIMENSIONS'
5080 ! include 'COMMON.VAR'
5081 ! include 'COMMON.CHAIN'
5082 ! include 'COMMON.DERIV'
5083 ! include 'COMMON.GEO'
5084 ! include 'COMMON.LOCAL'
5085 ! include 'COMMON.INTERACT'
5086 ! include 'COMMON.MD'
5087 ! include 'COMMON.SETUP'
5088 ! include 'COMMON.TIME1'
5089 ! integer licznik /0/
5092 ! include 'COMMON.LANGEVIN'
5094 ! include 'COMMON.LANGEVIN.lang0'
5096 ! include 'COMMON.IOUNITS'
5098 integer :: i,j,ind,ind1,m
5099 logical :: lprn = .false.
5100 real(kind=8) :: dtdi !el ,gamvec(2*nres)
5101 !el real(kind=8),dimension(2*nres,2*nres) :: ginvfric,fcopy
5102 ! real(kind=8),allocatable,dimension(:,:) :: fcopy
5103 !el real(kind=8),dimension(2*nres*(2*nres+1)/2) :: Ghalf !(mmaxres2) (mmaxres2=(maxres2*(maxres2+1)/2))
5104 !el common /syfek/ gamvec
5105 real(kind=8) :: work(8*2*nres)
5106 integer :: iwork(2*nres)
5107 !el common /przechowalnia/ ginvfric,Ghalf,fcopy
5108 integer :: ii,iti,k,l,nzero,nres2,nres6,ierr,mnum
5112 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
5113 if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
5114 if (fg_rank.ne.king) goto 10
5119 if(.not.allocated(gamvec)) allocate(gamvec(nres2)) !(MAXRES2)
5120 if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
5121 if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
5122 !el allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
5123 if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !maxres2=2*maxres
5125 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
5126 ! Zeroing out fricmat
5132 ! Load the friction coefficients corresponding to peptide groups
5137 gamvec(ind1)=gamp(mnum)
5139 ! Load the friction coefficients corresponding to side chains
5143 gamsc(ntyp1_molec(j),j)=1.0d0
5150 gamvec(ii)=gamsc(iabs(iti),mnum)
5152 if (surfarea) call sdarea(gamvec)
5154 ! write (iout,*) "Matrix A and vector gamma"
5156 ! write (iout,'(i2,$)') i
5158 ! write (iout,'(f4.1,$)') A(i,j)
5160 ! write (iout,'(f8.3)') gamvec(i)
5164 write (iout,*) "Vector gamvec"
5166 write (iout,'(i5,f10.5)') i, gamvec(i)
5170 ! The friction matrix
5175 dtdi=dtdi+A(j,k)*A(j,i)*gamvec(j)
5182 write (iout,'(//a)') "Matrix fricmat"
5183 call matout2(dimen,dimen,nres2,nres2,fricmat)
5185 if (lang.eq.2 .or. lang.eq.3) then
5186 ! Mass-scale the friction matrix if non-direct integration will be performed
5192 Ginvfric(i,j)=Ginvfric(i,j)+ &
5193 Gsqrm(i,k)*Gsqrm(l,j)*fricmat(k,l)
5198 ! Diagonalize the friction matrix
5203 Ghalf(ind)=Ginvfric(i,j)
5206 call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
5209 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
5210 " mass-scaled friction matrix"
5211 call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
5213 ! Precompute matrices for tinker stochastic integrator
5220 mt1(i,j)=mt1(i,j)+fricvec(k,i)*gsqrm(k,j)
5221 mt2(i,j)=mt2(i,j)+fricvec(k,i)*gsqrp(k,j)
5227 else if (lang.eq.4) then
5228 ! Diagonalize the friction matrix
5233 Ghalf(ind)=fricmat(i,j)
5236 call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
5239 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
5241 call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
5243 ! Determine the number of zero eigenvalues of the friction matrix
5244 nzero=max0(dimen-dimen1,0)
5245 ! do while (fricgam(nzero+1).le.1.0d-5 .and. nzero.lt.dimen)
5248 write (iout,*) "Number of zero eigenvalues:",nzero
5253 fricmat(i,j)=fricmat(i,j) &
5254 +fricvec(i,k)*fricvec(j,k)/fricgam(k)
5259 write (iout,'(//a)') "Generalized inverse of fricmat"
5260 call matout(dimen,dimen,nres6,nres6,fricmat)
5265 if (nfgtasks.gt.1) then
5266 if (fg_rank.eq.0) then
5267 ! The matching BROADCAST for fg processors is called in ERGASTULUM
5273 call MPI_Bcast(10,1,MPI_INTEGER,king,FG_COMM,IERROR)
5275 time_Bcast=time_Bcast+MPI_Wtime()-time00
5277 time_Bcast=time_Bcast+tcpu()-time00
5279 ! print *,"Processor",myrank,
5280 ! & " BROADCAST iorder in SETUP_FRICMAT"
5283 write (iout,*) "setup_fricmat licznik"!,licznik !sp
5289 ! Scatter the friction matrix
5290 call MPI_Scatterv(fricmat(1,1),nginv_counts(0),&
5291 nginv_start(0),MPI_DOUBLE_PRECISION,fcopy(1,1),&
5292 myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
5295 time_scatter=time_scatter+MPI_Wtime()-time00
5296 time_scatter_fmat=time_scatter_fmat+MPI_Wtime()-time00
5298 time_scatter=time_scatter+tcpu()-time00
5299 time_scatter_fmat=time_scatter_fmat+tcpu()-time00
5303 do j=1,2*my_ng_count
5304 fricmat(j,i)=fcopy(i,j)
5307 ! write (iout,*) "My chunk of fricmat"
5308 ! call MATOUT2(my_ng_count,dimen,maxres2,maxres2,fcopy)
5312 end subroutine setup_fricmat
5313 !-----------------------------------------------------------------------------
5314 subroutine stochastic_force(stochforcvec)
5317 use random, only:anorm_distr
5318 ! implicit real*8 (a-h,o-z)
5319 ! include 'DIMENSIONS'
5320 use control, only: tcpu
5325 ! include 'COMMON.VAR'
5326 ! include 'COMMON.CHAIN'
5327 ! include 'COMMON.DERIV'
5328 ! include 'COMMON.GEO'
5329 ! include 'COMMON.LOCAL'
5330 ! include 'COMMON.INTERACT'
5331 ! include 'COMMON.MD'
5332 ! include 'COMMON.TIME1'
5334 ! include 'COMMON.LANGEVIN'
5336 ! include 'COMMON.LANGEVIN.lang0'
5338 ! include 'COMMON.IOUNITS'
5340 real(kind=8) :: x,sig,lowb,highb
5341 real(kind=8) :: ff(3),force(3,0:2*nres),zeta2,lowb2
5342 real(kind=8) :: highb2,sig2,forcvec(6*nres),stochforcvec(6*nres)
5343 real(kind=8) :: time00
5344 logical :: lprn = .false.
5345 integer :: i,j,ind,mnum
5349 stochforc(j,i)=0.0d0
5359 ! Compute the stochastic forces acting on bodies. Store in force.
5365 force(j,i)=anorm_distr(x,sig,lowb,highb)
5373 force(j,i+nres)=anorm_distr(x,sig2,lowb2,highb2)
5377 time_fsample=time_fsample+MPI_Wtime()-time00
5379 time_fsample=time_fsample+tcpu()-time00
5381 ! Compute the stochastic forces acting on virtual-bond vectors.
5387 stochforc(j,i)=ff(j)+0.5d0*force(j,i)
5390 ff(j)=ff(j)+force(j,i)
5392 ! if (itype(i+1,1).ne.ntyp1) then
5394 if (itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
5396 stochforc(j,i)=stochforc(j,i)+force(j,i+nres+1)
5397 ff(j)=ff(j)+force(j,i+nres+1)
5402 stochforc(j,0)=ff(j)+force(j,nnt+nres)
5406 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5407 .and.(mnum.ne.5)) then
5408 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
5410 stochforc(j,i+nres)=force(j,i+nres)
5416 stochforcvec(j)=stochforc(j,0)
5421 stochforcvec(ind+j)=stochforc(j,i)
5427 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5428 .and.(mnum.ne.5)) then
5429 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
5431 stochforcvec(ind+j)=stochforc(j,i+nres)
5437 write (iout,*) "stochforcvec"
5439 write(iout,'(i5,e15.5)') i,stochforcvec(i)
5441 write(iout,*) "Stochastic forces backbone"
5443 write(iout,'(i5,3e15.5)') i,(stochforc(j,i),j=1,3)
5445 write(iout,*) "Stochastic forces side chain"
5447 write(iout,'(i5,3e15.5)') &
5448 i,(stochforc(j,i+nres),j=1,3)
5456 write (iout,*) i,ind
5458 forcvec(ind+j)=force(j,i)
5463 write (iout,*) i,ind
5465 forcvec(j+ind)=force(j,i+nres)
5470 write (iout,*) "forcvec"
5474 write (iout,'(2i3,2f10.5)') i,j,force(j,i),&
5481 write (iout,'(2i3,2f10.5)') i,j,force(j,i+nres),&
5490 end subroutine stochastic_force
5491 !-----------------------------------------------------------------------------
5492 subroutine sdarea(gamvec)
5494 ! Scale the friction coefficients according to solvent accessible surface areas
5495 ! Code adapted from TINKER
5499 ! implicit real*8 (a-h,o-z)
5500 ! include 'DIMENSIONS'
5501 ! include 'COMMON.CONTROL'
5502 ! include 'COMMON.VAR'
5503 ! include 'COMMON.MD'
5505 ! include 'COMMON.LANGEVIN'
5507 ! include 'COMMON.LANGEVIN.lang0'
5509 ! include 'COMMON.CHAIN'
5510 ! include 'COMMON.DERIV'
5511 ! include 'COMMON.GEO'
5512 ! include 'COMMON.LOCAL'
5513 ! include 'COMMON.INTERACT'
5514 ! include 'COMMON.IOUNITS'
5515 ! include 'COMMON.NAMES'
5516 real(kind=8),dimension(2*nres) :: radius,gamvec !(maxres2)
5517 real(kind=8),parameter :: twosix = 1.122462048309372981d0
5518 logical :: lprn = .false.
5519 real(kind=8) :: probe,area,ratio
5520 integer :: i,j,ind,iti,mnum
5522 ! determine new friction coefficients every few SD steps
5524 ! set the atomic radii to estimates of sigma values
5526 ! print *,"Entered sdarea"
5532 ! Load peptide group radii
5535 radius(i)=pstok(mnum)
5537 ! Load side chain radii
5541 radius(i+nres)=restok(iti,mnum)
5544 ! write (iout,*) "i",i," radius",radius(i)
5547 radius(i) = radius(i) / twosix
5548 if (radius(i) .ne. 0.0d0) radius(i) = radius(i) + probe
5551 ! scale atomic friction coefficients by accessible area
5553 if (lprn) write (iout,*) &
5554 "Original gammas, surface areas, scaling factors, new gammas, ",&
5555 "std's of stochastic forces"
5558 if (radius(i).gt.0.0d0) then
5559 call surfatom (i,area,radius)
5560 ratio = dmax1(area/(4.0d0*pi*radius(i)**2),1.0d-1)
5561 if (lprn) write (iout,'(i5,3f10.5,$)') &
5562 i,gamvec(ind+1),area,ratio
5565 gamvec(ind) = ratio * gamvec(ind)
5568 stdforcp(i)=stdfp(mnum)*dsqrt(gamvec(ind))
5569 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcp(i)
5573 if (radius(i+nres).gt.0.0d0) then
5574 call surfatom (i+nres,area,radius)
5575 ratio = dmax1(area/(4.0d0*pi*radius(i+nres)**2),1.0d-1)
5576 if (lprn) write (iout,'(i5,3f10.5,$)') &
5577 i,gamvec(ind+1),area,ratio
5580 gamvec(ind) = ratio * gamvec(ind)
5583 stdforcsc(i)=stdfsc(itype(i,mnum),mnum)*dsqrt(gamvec(ind))
5584 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcsc(i)
5589 end subroutine sdarea
5590 !-----------------------------------------------------------------------------
5592 !-----------------------------------------------------------------------------
5595 ! ###################################################
5596 ! ## COPYRIGHT (C) 1996 by Jay William Ponder ##
5597 ! ## All Rights Reserved ##
5598 ! ###################################################
5600 ! ################################################################
5602 ! ## subroutine surfatom -- exposed surface area of an atom ##
5604 ! ################################################################
5607 ! "surfatom" performs an analytical computation of the surface
5608 ! area of a specified atom; a simplified version of "surface"
5610 ! literature references:
5612 ! T. J. Richmond, "Solvent Accessible Surface Area and
5613 ! Excluded Volume in Proteins", Journal of Molecular Biology,
5616 ! L. Wesson and D. Eisenberg, "Atomic Solvation Parameters
5617 ! Applied to Molecular Dynamics of Proteins in Solution",
5618 ! Protein Science, 1, 227-235 (1992)
5620 ! variables and parameters:
5622 ! ir number of atom for which area is desired
5623 ! area accessible surface area of the atom
5624 ! radius radii of each of the individual atoms
5627 subroutine surfatom(ir,area,radius)
5629 ! implicit real*8 (a-h,o-z)
5630 ! include 'DIMENSIONS'
5632 ! include 'COMMON.GEO'
5633 ! include 'COMMON.IOUNITS'
5635 integer :: nsup,nstart_sup
5636 ! double precision c,dc,dc_old,d_c_work,xloc,xrot,dc_norm
5637 ! common /chain/ c(3,maxres2+2),dc(3,0:maxres2),dc_old(3,0:maxres2),
5638 ! & xloc(3,maxres),xrot(3,maxres),dc_norm(3,0:maxres2),
5639 ! & dc_work(MAXRES6),nres,nres0
5640 integer,parameter :: maxarc=300
5644 integer :: mi,ni,narc
5645 integer :: key(maxarc)
5646 integer :: intag(maxarc)
5647 integer :: intag1(maxarc)
5648 real(kind=8) :: area,arcsum
5649 real(kind=8) :: arclen,exang
5650 real(kind=8) :: delta,delta2
5651 real(kind=8) :: eps,rmove
5652 real(kind=8) :: xr,yr,zr
5653 real(kind=8) :: rr,rrsq
5654 real(kind=8) :: rplus,rminus
5655 real(kind=8) :: axx,axy,axz
5656 real(kind=8) :: ayx,ayy
5657 real(kind=8) :: azx,azy,azz
5658 real(kind=8) :: uxj,uyj,uzj
5659 real(kind=8) :: tx,ty,tz
5660 real(kind=8) :: txb,tyb,td
5661 real(kind=8) :: tr2,tr,txr,tyr
5662 real(kind=8) :: tk1,tk2
5663 real(kind=8) :: thec,the,t,tb
5664 real(kind=8) :: txk,tyk,tzk
5665 real(kind=8) :: t1,ti,tf,tt
5666 real(kind=8) :: txj,tyj,tzj
5667 real(kind=8) :: ccsq,cc,xysq
5668 real(kind=8) :: bsqk,bk,cosine
5669 real(kind=8) :: dsqj,gi,pix2
5670 real(kind=8) :: therk,dk,gk
5671 real(kind=8) :: risqk,rik
5672 real(kind=8) :: radius(2*nres) !(maxatm) (maxatm=maxres2)
5673 real(kind=8) :: ri(maxarc),risq(maxarc)
5674 real(kind=8) :: ux(maxarc),uy(maxarc),uz(maxarc)
5675 real(kind=8) :: xc(maxarc),yc(maxarc),zc(maxarc)
5676 real(kind=8) :: xc1(maxarc),yc1(maxarc),zc1(maxarc)
5677 real(kind=8) :: dsq(maxarc),bsq(maxarc)
5678 real(kind=8) :: dsq1(maxarc),bsq1(maxarc)
5679 real(kind=8) :: arci(maxarc),arcf(maxarc)
5680 real(kind=8) :: ex(maxarc),lt(maxarc),gr(maxarc)
5681 real(kind=8) :: b(maxarc),b1(maxarc),bg(maxarc)
5682 real(kind=8) :: kent(maxarc),kout(maxarc)
5683 real(kind=8) :: ther(maxarc)
5684 logical :: moved,top
5685 logical :: omit(maxarc)
5688 ! maxatm = 2*nres !maxres2 maxres2=2*maxres
5689 ! maxlight = 8*maxatm
5692 ! maxtors = 4*maxatm
5694 ! zero out the surface area for the sphere of interest
5697 ! write (2,*) "ir",ir," radius",radius(ir)
5698 if (radius(ir) .eq. 0.0d0) return
5700 ! set the overlap significance and connectivity shift
5704 delta2 = delta * delta
5709 ! store coordinates and radius of the sphere of interest
5717 ! initialize values of some counters and summations
5726 ! test each sphere to see if it overlaps the sphere of interest
5729 if (i.eq.ir .or. radius(i).eq.0.0d0) goto 30
5730 rplus = rr + radius(i)
5732 if (abs(tx) .ge. rplus) goto 30
5734 if (abs(ty) .ge. rplus) goto 30
5736 if (abs(tz) .ge. rplus) goto 30
5738 ! check for sphere overlap by testing distance against radii
5740 xysq = tx*tx + ty*ty
5741 if (xysq .lt. delta2) then
5748 if (rplus-cc .le. delta) goto 30
5749 rminus = rr - radius(i)
5751 ! check to see if sphere of interest is completely buried
5753 if (cc-abs(rminus) .le. delta) then
5754 if (rminus .le. 0.0d0) goto 170
5758 ! check for too many overlaps with sphere of interest
5760 if (io .ge. maxarc) then
5762 20 format (/,' SURFATOM -- Increase the Value of MAXARC')
5766 ! get overlap between current sphere and sphere of interest
5775 gr(io) = (ccsq+rplus*rminus) / (2.0d0*rr*b1(io))
5781 ! case where no other spheres overlap the sphere of interest
5784 area = 4.0d0 * pi * rrsq
5788 ! case where only one sphere overlaps the sphere of interest
5791 area = pix2 * (1.0d0 + gr(1))
5792 area = mod(area,4.0d0*pi) * rrsq
5796 ! case where many spheres intersect the sphere of interest;
5797 ! sort the intersecting spheres by their degree of overlap
5799 call sort2 (io,gr,key)
5802 intag(i) = intag1(k)
5811 ! get radius of each overlap circle on surface of the sphere
5816 risq(i) = rrsq - gi*gi
5817 ri(i) = sqrt(risq(i))
5818 ther(i) = 0.5d0*pi - asin(min(1.0d0,max(-1.0d0,gr(i))))
5821 ! find boundary of inaccessible area on sphere of interest
5824 if (.not. omit(k)) then
5831 ! check to see if J circle is intersecting K circle;
5832 ! get distance between circle centers and sum of radii
5835 if (omit(j)) goto 60
5836 cc = (txk*xc(j)+tyk*yc(j)+tzk*zc(j))/(bk*b(j))
5837 cc = acos(min(1.0d0,max(-1.0d0,cc)))
5838 td = therk + ther(j)
5840 ! check to see if circles enclose separate regions
5842 if (cc .ge. td) goto 60
5844 ! check for circle J completely inside circle K
5846 if (cc+ther(j) .lt. therk) goto 40
5848 ! check for circles that are essentially parallel
5850 if (cc .gt. delta) goto 50
5855 ! check to see if sphere of interest is completely buried
5858 if (pix2-cc .le. td) goto 170
5864 ! find T value of circle intersections
5867 if (omit(k)) goto 110
5882 ! rotation matrix elements
5894 if (.not. omit(j)) then
5899 ! rotate spheres so K vector colinear with z-axis
5901 uxj = txj*axx + tyj*axy - tzj*axz
5902 uyj = tyj*ayy - txj*ayx
5903 uzj = txj*azx + tyj*azy + tzj*azz
5904 cosine = min(1.0d0,max(-1.0d0,uzj/b(j)))
5905 if (acos(cosine) .lt. therk+ther(j)) then
5906 dsqj = uxj*uxj + uyj*uyj
5911 tr2 = risqk*dsqj - tb*tb
5917 ! get T values of intersection for K circle
5920 tb = min(1.0d0,max(-1.0d0,tb))
5922 if (tyb-txr .lt. 0.0d0) tk1 = pix2 - tk1
5924 tb = min(1.0d0,max(-1.0d0,tb))
5926 if (tyb+txr .lt. 0.0d0) tk2 = pix2 - tk2
5927 thec = (rrsq*uzj-gk*bg(j)) / (rik*ri(j)*b(j))
5928 if (abs(thec) .lt. 1.0d0) then
5930 else if (thec .ge. 1.0d0) then
5932 else if (thec .le. -1.0d0) then
5936 ! see if "tk1" is entry or exit point; check t=0 point;
5937 ! "ti" is exit point, "tf" is entry point
5939 cosine = min(1.0d0,max(-1.0d0, &
5940 (uzj*gk-uxj*rik)/(b(j)*rr)))
5941 if ((acos(cosine)-ther(j))*(tk2-tk1) .le. 0.0d0) then
5949 if (narc .ge. maxarc) then
5951 70 format (/,' SURFATOM -- Increase the Value',&
5955 if (tf .le. ti) then
5976 ! special case; K circle without intersections
5978 if (narc .le. 0) goto 90
5980 ! general case; sum up arclength and set connectivity code
5982 call sort2 (narc,arci,key)
5987 if (narc .gt. 1) then
5990 if (t .lt. arci(j)) then
5991 arcsum = arcsum + arci(j) - t
5992 exang = exang + ex(ni)
5994 if (jb .ge. maxarc) then
5996 80 format (/,' SURFATOM -- Increase the Value',&
6001 kent(jb) = maxarc*i + k
6003 kout(jb) = maxarc*k + i
6012 arcsum = arcsum + pix2 - t
6014 exang = exang + ex(ni)
6017 kent(jb) = maxarc*i + k
6019 kout(jb) = maxarc*k + i
6026 arclen = arclen + gr(k)*arcsum
6029 if (arclen .eq. 0.0d0) goto 170
6030 if (jb .eq. 0) goto 150
6032 ! find number of independent boundaries and check connectivity
6036 if (kout(k) .ne. 0) then
6043 if (m .eq. kent(ii)) then
6046 if (j .eq. jb) goto 150
6058 ! attempt to fix connectivity error by moving atom slightly
6062 140 format (/,' SURFATOM -- Connectivity Error at Atom',i6)
6071 ! compute the exposed surface area for the sphere of interest
6074 area = ib*pix2 + exang + arclen
6075 area = mod(area,4.0d0*pi) * rrsq
6077 ! attempt to fix negative area by moving atom slightly
6079 if (area .lt. 0.0d0) then
6082 160 format (/,' SURFATOM -- Negative Area at Atom',i6)
6093 end subroutine surfatom
6094 !----------------------------------------------------------------
6095 !----------------------------------------------------------------
6096 subroutine alloc_MD_arrays
6097 !EL Allocation of arrays used by MD module
6099 integer :: nres2,nres6
6102 !----------------------
6106 allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
6107 allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
6108 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
6109 allocate(fricvec(nres2,nres2))
6110 allocate(pfric_mat(nres2,nres2),vfric_mat(nres2,nres2))
6111 allocate(afric_mat(nres2,nres2),prand_mat(nres2,nres2))
6112 allocate(vrand_mat1(nres2,nres2),vrand_mat2(nres2,nres2)) !(MAXRES2,MAXRES2)
6113 allocate(pfric0_mat(nres2,nres2,0:maxflag_stoch))
6114 allocate(afric0_mat(nres2,nres2,0:maxflag_stoch))
6115 allocate(vfric0_mat(nres2,nres2,0:maxflag_stoch))
6116 allocate(prand0_mat(nres2,nres2,0:maxflag_stoch))
6117 allocate(vrand0_mat1(nres2,nres2,0:maxflag_stoch))
6118 allocate(vrand0_mat2(nres2,nres2,0:maxflag_stoch)) !(MAXRES2,MAXRES2,0:maxflag_stoch)
6119 allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
6121 allocate(mt1(nres2,nres2),mt2(nres2,nres2),mt3(nres2,nres2)) !(maxres2,maxres2)
6122 !----------------------
6124 ! commom.langevin.lang0
6126 allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
6127 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
6128 allocate(fricvec(nres2,nres2)) !(MAXRES2,MAXRES2)
6129 allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
6130 allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
6133 if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2))
6134 !----------------------
6135 ! commom.hairpin in CSA module
6136 !----------------------
6137 ! common.mce in MCM_MD module
6138 !----------------------
6140 ! common /mdgrad/ in module.energy
6141 ! common /back_constr/ in module.energy
6142 ! common /qmeas/ in module.energy
6145 allocate(potEcomp(0:n_ene+4)) !(0:n_ene+4)
6147 allocate(d_t(3,0:nres2),d_a(3,0:nres2),d_t_old(3,0:nres2)) !(3,0:MAXRES2)
6148 allocate(d_a_work(nres6)) !(6*MAXRES)
6150 allocate(DM(nres2),DU1(nres2),DU2(nres2))
6151 allocate(DMorig(nres2),DU1orig(nres2),DU2orig(nres2))
6153 allocate(Gmat(nres2,nres2),A(nres2,nres2))
6154 if(.not.allocated(Ginv)) allocate(Ginv(nres2,nres2)) !in control: ergastulum
6155 allocate(Gsqrp(nres2,nres2),Gsqrm(nres2,nres2),Gvec(nres2,nres2)) !(maxres2,maxres2)
6157 allocate(Geigen(nres2)) !(maxres2)
6158 if(.not.allocated(vtot)) allocate(vtot(nres2)) !(maxres2)
6159 ! common /inertia/ in io_conf: parmread
6160 ! real(kind=8),dimension(:),allocatable :: ISC,msc !(ntyp+1)
6161 ! common /langevin/in io read_MDpar
6162 ! real(kind=8),dimension(:),allocatable :: gamsc !(ntyp1)
6163 ! real(kind=8),dimension(:),allocatable :: stdfsc !(ntyp)
6164 ! in io_conf: parmread
6165 ! real(kind=8),dimension(:),allocatable :: restok !(ntyp+1)
6166 ! common /mdpmpi/ in control: ergastulum
6167 if(.not.allocated(ng_start)) allocate(ng_start(0:nfgtasks-1))
6168 if(.not.allocated(ng_counts)) allocate(ng_counts(0:nfgtasks-1))
6169 if(.not.allocated(nginv_counts)) allocate(nginv_counts(0:nfgtasks-1)) !(0:MaxProcs-1)
6170 if(.not.allocated(nginv_start)) allocate(nginv_start(0:nfgtasks)) !(0:MaxProcs)
6171 !----------------------
6172 ! common.muca in read_muca
6173 ! common /double_muca/
6174 ! real(kind=8) :: elow,ehigh,factor,hbin,factor_min
6175 ! real(kind=8),dimension(:),allocatable :: emuca,nemuca,&
6176 ! nemuca2,hist !(4*maxres)
6177 ! common /integer_muca/
6178 ! integer :: nmuca,imtime,muca_smooth
6180 ! real(kind=8),dimension(:),allocatable :: elowi,ehighi !(maxprocs)
6181 !----------------------
6183 ! common /mdgrad/ in module.energy
6184 ! common /back_constr/ in module.energy
6185 ! common /qmeas/ in module.energy
6189 allocate(d_t_work(nres6),d_t_work_new(nres6),d_af_work(nres6))
6190 allocate(d_as_work(nres6),kinetic_force(nres6)) !(MAXRES6)
6191 allocate(d_t_new(3,0:nres2),d_a_old(3,0:nres2),d_a_short(3,0:nres2)) !,d_a !(3,0:MAXRES2)
6192 allocate(stdforcp(nres),stdforcsc(nres)) !(MAXRES)
6193 !----------------------
6195 allocate(D_ban(nres6)) !(MAXRES6) maxres6=6*maxres
6196 ! common /stochcalc/ stochforcvec
6197 allocate(stochforcvec(nres6)) !(MAXRES6) maxres6=6*maxres
6198 !----------------------
6200 end subroutine alloc_MD_arrays
6201 !-----------------------------------------------------------------------------
6202 !-----------------------------------------------------------------------------