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) then
1004 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1.and.mnum.ne.5) then
1007 v_work(ind)=d_t(j,i+nres)
1012 write (66,'(80f10.5)') &
1013 ((d_t(j,i),j=1,3),i=0,nres-1),((d_t(j,i+nres),j=1,3),i=1,nres)
1017 v_transf(i)=v_transf(i)+gvec(j,i)*v_work(j)
1019 v_transf(i)= v_transf(i)*dsqrt(geigen(i))
1021 write (67,'(80f10.5)') (v_transf(i),i=1,ind)
1024 if (mod(itime,ntwx).eq.0) then
1026 write (tytul,'("time",f8.2)') totT
1028 call hairpin(.true.,nharp,iharp)
1029 call secondary2(.true.)
1030 call pdbout(potE,tytul,ipdb)
1035 if (rstcount.eq.1000.or.itime.eq.n_timestep) then
1036 open(irest2,file=rest2name,status='unknown')
1037 write(irest2,*) totT,EK,potE,totE,t_bath
1039 ! AL 4/17/17: Now writing d_t(0,:) too
1041 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
1043 ! AL 4/17/17: Now writing d_c(0,:) too
1045 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
1053 t_MD=MPI_Wtime()-tt0
1057 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))') &
1059 'MD calculations setup:',t_MDsetup,&
1060 'Energy & gradient evaluation:',t_enegrad,&
1061 'Stochastic MD setup:',t_langsetup,&
1062 'Stochastic MD step setup:',t_sdsetup,&
1064 write (iout,'(/28(1h=),a25,27(1h=))') &
1065 ' End of MD calculation '
1067 write (iout,*) "time for etotal",t_etotal," elong",t_elong,&
1069 write (iout,*) "time_fric",time_fric," time_stoch",time_stoch,&
1070 " time_fricmatmult",time_fricmatmult," time_fsample ",&
1075 !-----------------------------------------------------------------------------
1076 subroutine velverlet_step(itime)
1077 !-------------------------------------------------------------------------------
1078 ! Perform a single velocity Verlet step; the time step can be rescaled if
1079 ! increments in accelerations exceed the threshold
1080 !-------------------------------------------------------------------------------
1081 ! implicit real*8 (a-h,o-z)
1082 ! include 'DIMENSIONS'
1084 use control, only:tcpu
1088 integer :: ierror,ierrcode
1089 real(kind=8) :: errcode
1091 ! include 'COMMON.SETUP'
1092 ! include 'COMMON.VAR'
1093 ! include 'COMMON.MD'
1095 ! include 'COMMON.LANGEVIN'
1097 ! include 'COMMON.LANGEVIN.lang0'
1099 ! include 'COMMON.CHAIN'
1100 ! include 'COMMON.DERIV'
1101 ! include 'COMMON.GEO'
1102 ! include 'COMMON.LOCAL'
1103 ! include 'COMMON.INTERACT'
1104 ! include 'COMMON.IOUNITS'
1105 ! include 'COMMON.NAMES'
1106 ! include 'COMMON.TIME1'
1107 ! include 'COMMON.MUCA'
1108 real(kind=8),dimension(3) :: vcm,incr
1109 real(kind=8),dimension(3) :: L
1110 integer :: count,rstcount !ilen,
1112 character(len=50) :: tytul
1113 integer :: maxcount_scale = 20
1114 !el common /gucio/ cm
1115 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1116 !el common /stochcalc/ stochforcvec
1117 integer :: itime,icount_scale,itime_scal,i,j,ifac_time
1119 real(kind=8) :: epdrift,tt0,fac_time
1121 if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres)) !(MAXRES6) maxres6=6*maxres
1127 else if (lang.eq.2 .or. lang.eq.3) then
1129 call stochastic_force(stochforcvec)
1132 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
1134 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1141 icount_scale=icount_scale+1
1142 if (icount_scale.gt.maxcount_scale) then
1144 "ERROR: too many attempts at scaling down the time step. ",&
1145 "amax=",amax,"epdrift=",epdrift,&
1146 "damax=",damax,"edriftmax=",edriftmax,&
1150 call MPI_Abort(MPI_COMM_WORLD,IERROR,IERRCODE)
1154 ! First step of the velocity Verlet algorithm
1159 else if (lang.eq.3) then
1161 call sd_verlet1_ciccotti
1163 else if (lang.eq.1) then
1168 ! Build the chain from the newly calculated coordinates
1169 call chainbuild_cart
1170 if (rattle) call rattle1
1172 if (large.and. mod(itime,ntwe).eq.0) then
1173 write (iout,*) "Cartesian and internal coordinates: step 1"
1178 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),&
1179 (dc(j,i+nres),j=1,3)
1181 write (iout,*) "Accelerations"
1183 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1184 (d_a(j,i+nres),j=1,3)
1186 write (iout,*) "Velocities, step 1"
1188 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1189 (d_t(j,i+nres),j=1,3)
1198 ! Calculate energy and forces
1200 call etotal(potEcomp)
1201 ! AL 4/17/17: Reduce the steps if NaNs occurred.
1202 if (potEcomp(0).gt.0.99e20 .or. isnan(potEcomp(0)).gt.0) then
1207 if (large.and. mod(itime,ntwe).eq.0) &
1208 call enerprint(potEcomp)
1211 t_etotal=t_etotal+MPI_Wtime()-tt0
1213 t_etotal=t_etotal+tcpu()-tt0
1216 potE=potEcomp(0)-potEcomp(20)
1218 ! Get the new accelerations
1221 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1223 t_enegrad=t_enegrad+tcpu()-tt0
1225 ! Determine maximum acceleration and scale down the timestep if needed
1227 amax=amax/(itime_scal+1)**2
1228 call predict_edrift(epdrift)
1229 if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then
1230 ! Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
1232 ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax)) &
1234 itime_scal=itime_scal+ifac_time
1235 ! fac_time=dmin1(damax/amax,0.5d0)
1236 fac_time=0.5d0**ifac_time
1237 d_time=d_time*fac_time
1238 if (lang.eq.2 .or. lang.eq.3) then
1240 ! write (iout,*) "Calling sd_verlet_setup: 1"
1241 ! Rescale the stochastic forces and recalculate or restore
1242 ! the matrices of tinker integrator
1243 if (itime_scal.gt.maxflag_stoch) then
1244 if (large) write (iout,'(a,i5,a)') &
1245 "Calculate matrices for stochastic step;",&
1246 " itime_scal ",itime_scal
1248 call sd_verlet_p_setup
1250 call sd_verlet_ciccotti_setup
1252 write (iout,'(2a,i3,a,i3,1h.)') &
1253 "Warning: cannot store matrices for stochastic",&
1254 " integration because the index",itime_scal,&
1255 " is greater than",maxflag_stoch
1256 write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",&
1257 " integration Langevin algorithm for better efficiency."
1258 else if (flag_stoch(itime_scal)) then
1259 if (large) write (iout,'(a,i5,a,l1)') &
1260 "Restore matrices for stochastic step; itime_scal ",&
1261 itime_scal," flag ",flag_stoch(itime_scal)
1264 pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
1265 afric_mat(i,j)=afric0_mat(i,j,itime_scal)
1266 vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
1267 prand_mat(i,j)=prand0_mat(i,j,itime_scal)
1268 vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
1269 vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
1273 if (large) write (iout,'(2a,i5,a,l1)') &
1274 "Calculate & store matrices for stochastic step;",&
1275 " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
1277 call sd_verlet_p_setup
1279 call sd_verlet_ciccotti_setup
1281 flag_stoch(ifac_time)=.true.
1284 pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
1285 afric0_mat(i,j,itime_scal)=afric_mat(i,j)
1286 vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
1287 prand0_mat(i,j,itime_scal)=prand_mat(i,j)
1288 vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
1289 vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
1293 fac_time=1.0d0/dsqrt(fac_time)
1295 stochforcvec(i)=fac_time*stochforcvec(i)
1298 else if (lang.eq.1) then
1299 ! Rescale the accelerations due to stochastic forces
1300 fac_time=1.0d0/dsqrt(fac_time)
1302 d_as_work(i)=d_as_work(i)*fac_time
1305 if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)') &
1306 "itime",itime," Timestep scaled down to ",&
1307 d_time," ifac_time",ifac_time," itime_scal",itime_scal
1309 ! Second step of the velocity Verlet algorithm
1314 else if (lang.eq.3) then
1316 call sd_verlet2_ciccotti
1318 else if (lang.eq.1) then
1323 if (rattle) call rattle2
1326 if (d_time.ne.d_time0) then
1329 if (lang.eq.2 .or. lang.eq.3) then
1330 if (large) write (iout,'(a)') &
1331 "Restore original matrices for stochastic step"
1332 ! write (iout,*) "Calling sd_verlet_setup: 2"
1333 ! Restore the matrices of tinker integrator if the time step has been restored
1336 pfric_mat(i,j)=pfric0_mat(i,j,0)
1337 afric_mat(i,j)=afric0_mat(i,j,0)
1338 vfric_mat(i,j)=vfric0_mat(i,j,0)
1339 prand_mat(i,j)=prand0_mat(i,j,0)
1340 vrand_mat1(i,j)=vrand0_mat1(i,j,0)
1341 vrand_mat2(i,j)=vrand0_mat2(i,j,0)
1350 ! Calculate the kinetic and the total energy and the kinetic temperature
1354 ! call kinetic1(EK1)
1355 ! write (iout,*) "step",itime," EK",EK," EK1",EK1
1357 ! Couple the system to Berendsen bath if needed
1358 if (tbf .and. lang.eq.0) then
1361 kinetic_T=2.0d0/(dimen3*Rb)*EK
1362 ! Backup the coordinates, velocities, and accelerations
1366 d_t_old(j,i)=d_t(j,i)
1367 d_a_old(j,i)=d_a(j,i)
1371 if (mod(itime,ntwe).eq.0 .and. large) then
1372 write (iout,*) "Velocities, step 2"
1374 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1375 (d_t(j,i+nres),j=1,3)
1380 end subroutine velverlet_step
1381 !-----------------------------------------------------------------------------
1382 subroutine RESPA_step(itime)
1383 !-------------------------------------------------------------------------------
1384 ! Perform a single RESPA step.
1385 !-------------------------------------------------------------------------------
1386 ! implicit real*8 (a-h,o-z)
1387 ! include 'DIMENSIONS'
1391 use control, only:tcpu
1393 ! use io_conf, only:cartprint
1396 integer :: IERROR,ERRCODE
1398 ! include 'COMMON.SETUP'
1399 ! include 'COMMON.CONTROL'
1400 ! include 'COMMON.VAR'
1401 ! include 'COMMON.MD'
1403 ! include 'COMMON.LANGEVIN'
1405 ! include 'COMMON.LANGEVIN.lang0'
1407 ! include 'COMMON.CHAIN'
1408 ! include 'COMMON.DERIV'
1409 ! include 'COMMON.GEO'
1410 ! include 'COMMON.LOCAL'
1411 ! include 'COMMON.INTERACT'
1412 ! include 'COMMON.IOUNITS'
1413 ! include 'COMMON.NAMES'
1414 ! include 'COMMON.TIME1'
1415 real(kind=8),dimension(0:n_ene) :: energia_short,energia_long
1416 real(kind=8),dimension(3) :: L,vcm,incr
1417 real(kind=8),dimension(3,0:2*nres) :: dc_old0,d_t_old0,d_a_old0 !(3,0:maxres2) maxres2=2*maxres
1418 logical :: PRINT_AMTS_MSG = .false.
1419 integer :: count,rstcount !ilen,
1421 character(len=50) :: tytul
1422 integer :: maxcount_scale = 10
1423 !el common /gucio/ cm
1424 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1425 !el common /stochcalc/ stochforcvec
1426 integer :: itime,itt,i,j,itsplit
1428 !el common /cipiszcze/ itt
1430 real(kind=8) :: epdrift,tt0,epdriftmax
1433 if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres)) !(MAXRES6) maxres6=6*maxres
1437 if (large.and. mod(itime,ntwe).eq.0) then
1438 write (iout,*) "***************** RESPA itime",itime
1439 write (iout,*) "Cartesian and internal coordinates: step 0"
1441 call pdbout(0.0d0,"cipiszcze",iout)
1443 write (iout,*) "Accelerations from long-range forces"
1445 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1446 (d_a(j,i+nres),j=1,3)
1448 write (iout,*) "Velocities, step 0"
1450 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1451 (d_t(j,i+nres),j=1,3)
1456 ! Perform the initial RESPA step (increment velocities)
1457 ! write (iout,*) "*********************** RESPA ini"
1460 if (mod(itime,ntwe).eq.0 .and. large) then
1461 write (iout,*) "Velocities, end"
1463 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1464 (d_t(j,i+nres),j=1,3)
1468 ! Compute the short-range forces
1474 ! 7/2/2009 commented out
1476 ! call etotal_short(energia_short)
1479 ! 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
1482 d_a(j,i)=d_a_short(j,i)
1486 if (large.and. mod(itime,ntwe).eq.0) then
1487 write (iout,*) "energia_short",energia_short(0)
1488 write (iout,*) "Accelerations from short-range forces"
1490 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1491 (d_a(j,i+nres),j=1,3)
1496 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1498 t_enegrad=t_enegrad+tcpu()-tt0
1503 d_t_old(j,i)=d_t(j,i)
1504 d_a_old(j,i)=d_a(j,i)
1507 ! 6/30/08 A-MTS: attempt at increasing the split number
1510 dc_old0(j,i)=dc_old(j,i)
1511 d_t_old0(j,i)=d_t_old(j,i)
1512 d_a_old0(j,i)=d_a_old(j,i)
1515 if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
1516 if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
1523 ! write (iout,*) "itime",itime," ntime_split",ntime_split
1524 ! Split the time step
1525 d_time=d_time0/ntime_split
1526 ! Perform the short-range RESPA steps (velocity Verlet increments of
1527 ! positions and velocities using short-range forces)
1528 ! write (iout,*) "*********************** RESPA split"
1529 do itsplit=1,ntime_split
1532 else if (lang.eq.2 .or. lang.eq.3) then
1534 call stochastic_force(stochforcvec)
1537 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
1539 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1544 ! First step of the velocity Verlet algorithm
1549 else if (lang.eq.3) then
1551 call sd_verlet1_ciccotti
1553 else if (lang.eq.1) then
1558 ! Build the chain from the newly calculated coordinates
1559 call chainbuild_cart
1560 if (rattle) call rattle1
1562 if (large.and. mod(itime,ntwe).eq.0) then
1563 write (iout,*) "***** ITSPLIT",itsplit
1564 write (iout,*) "Cartesian and internal coordinates: step 1"
1565 call pdbout(0.0d0,"cipiszcze",iout)
1568 write (iout,*) "Velocities, step 1"
1570 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1571 (d_t(j,i+nres),j=1,3)
1580 ! Calculate energy and forces
1582 call etotal_short(energia_short)
1583 ! AL 4/17/17: Exit itime_split loop when energy goes infinite
1584 if (energia_short(0).gt.0.99e20 .or. isnan(energia_short(0)) ) then
1585 if (PRINT_AMTS_MSG) &
1586 write (iout,*) "Infinities/NaNs in energia_short",energia_short(0),"; increasing ntime_split to",ntime_split
1587 ntime_split=ntime_split*2
1588 if (ntime_split.gt.maxtime_split) then
1591 "Cannot rescue the run; aborting job. Retry with a smaller time step"
1593 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1596 "Cannot rescue the run; terminating. Retry with a smaller time step"
1602 if (large.and. mod(itime,ntwe).eq.0) &
1603 call enerprint(energia_short)
1606 t_eshort=t_eshort+MPI_Wtime()-tt0
1608 t_eshort=t_eshort+tcpu()-tt0
1612 ! Get the new accelerations
1614 ! 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
1617 d_a_short(j,i)=d_a(j,i)
1621 if (large.and. mod(itime,ntwe).eq.0) then
1622 write (iout,*)"energia_short",energia_short(0)
1623 write (iout,*) "Accelerations from short-range forces"
1625 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1626 (d_a(j,i+nres),j=1,3)
1631 ! Determine maximum acceleration and scale down the timestep if needed
1633 amax=amax/ntime_split**2
1634 call predict_edrift(epdrift)
1635 if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0) &
1636 write (iout,*) "amax",amax," damax",damax,&
1637 " epdrift",epdrift," epdriftmax",epdriftmax
1638 ! Exit loop and try with increased split number if the change of
1639 ! acceleration is too big
1640 if (amax.gt.damax .or. epdrift.gt.edriftmax) then
1641 if (ntime_split.lt.maxtime_split) then
1643 ntime_split=ntime_split*2
1644 ! AL 4/17/17: We should exit the itime_split loop when acceleration change is too big
1648 dc_old(j,i)=dc_old0(j,i)
1649 d_t_old(j,i)=d_t_old0(j,i)
1650 d_a_old(j,i)=d_a_old0(j,i)
1653 if (PRINT_AMTS_MSG) then
1654 write (iout,*) "acceleration/energy drift too large",amax,&
1655 epdrift," split increased to ",ntime_split," itime",itime,&
1661 "Uh-hu. Bumpy landscape. Maximum splitting number",&
1663 " already reached!!! Trying to carry on!"
1667 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1669 t_enegrad=t_enegrad+tcpu()-tt0
1671 ! Second step of the velocity Verlet algorithm
1676 else if (lang.eq.3) then
1678 call sd_verlet2_ciccotti
1680 else if (lang.eq.1) then
1685 if (rattle) call rattle2
1686 ! Backup the coordinates, velocities, and accelerations
1690 d_t_old(j,i)=d_t(j,i)
1691 d_a_old(j,i)=d_a(j,i)
1698 ! Restore the time step
1700 ! Compute long-range forces
1707 call etotal_long(energia_long)
1708 if (energia_long(0).gt.0.99e20 .or. isnan(energia_long(0))) then
1711 "Infinitied/NaNs in energia_long, Aborting MPI job."
1713 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1715 write (iout,*) "Infinitied/NaNs in energia_long, terminating."
1719 if (large.and. mod(itime,ntwe).eq.0) &
1720 call enerprint(energia_long)
1723 t_elong=t_elong+MPI_Wtime()-tt0
1725 t_elong=t_elong+tcpu()-tt0
1731 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1733 t_enegrad=t_enegrad+tcpu()-tt0
1735 ! Compute accelerations from long-range forces
1737 if (large.and. mod(itime,ntwe).eq.0) then
1738 write (iout,*) "energia_long",energia_long(0)
1739 write (iout,*) "Cartesian and internal coordinates: step 2"
1741 call pdbout(0.0d0,"cipiszcze",iout)
1743 write (iout,*) "Accelerations from long-range forces"
1745 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1746 (d_a(j,i+nres),j=1,3)
1748 write (iout,*) "Velocities, step 2"
1750 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1751 (d_t(j,i+nres),j=1,3)
1755 ! Compute the final RESPA step (increment velocities)
1756 ! write (iout,*) "*********************** RESPA fin"
1758 ! Compute the complete potential energy
1760 potEcomp(i)=energia_short(i)+energia_long(i)
1762 potE=potEcomp(0)-potEcomp(20)
1763 ! potE=energia_short(0)+energia_long(0)
1766 ! Calculate the kinetic and the total energy and the kinetic temperature
1769 ! Couple the system to Berendsen bath if needed
1770 if (tbf .and. lang.eq.0) then
1773 kinetic_T=2.0d0/(dimen3*Rb)*EK
1774 ! Backup the coordinates, velocities, and accelerations
1776 if (mod(itime,ntwe).eq.0 .and. large) then
1777 write (iout,*) "Velocities, end"
1779 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1780 (d_t(j,i+nres),j=1,3)
1785 end subroutine RESPA_step
1786 !-----------------------------------------------------------------------------
1787 subroutine RESPA_vel
1788 ! First and last RESPA step (incrementing velocities using long-range
1791 ! implicit real*8 (a-h,o-z)
1792 ! include 'DIMENSIONS'
1793 ! include 'COMMON.CONTROL'
1794 ! include 'COMMON.VAR'
1795 ! include 'COMMON.MD'
1796 ! include 'COMMON.CHAIN'
1797 ! include 'COMMON.DERIV'
1798 ! include 'COMMON.GEO'
1799 ! include 'COMMON.LOCAL'
1800 ! include 'COMMON.INTERACT'
1801 ! include 'COMMON.IOUNITS'
1802 ! include 'COMMON.NAMES'
1803 integer :: i,j,inres,mnum
1806 d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
1810 d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
1815 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
1816 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
1817 .and.(mnum.ne.5)) then
1820 d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
1825 end subroutine RESPA_vel
1826 !-----------------------------------------------------------------------------
1828 ! Applying velocity Verlet algorithm - step 1 to coordinates
1830 ! implicit real*8 (a-h,o-z)
1831 ! include 'DIMENSIONS'
1832 ! include 'COMMON.CONTROL'
1833 ! include 'COMMON.VAR'
1834 ! include 'COMMON.MD'
1835 ! include 'COMMON.CHAIN'
1836 ! include 'COMMON.DERIV'
1837 ! include 'COMMON.GEO'
1838 ! include 'COMMON.LOCAL'
1839 ! include 'COMMON.INTERACT'
1840 ! include 'COMMON.IOUNITS'
1841 ! include 'COMMON.NAMES'
1842 real(kind=8) :: adt,adt2
1843 integer :: i,j,inres,mnum
1846 write (iout,*) "VELVERLET1 START: DC"
1848 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),&
1849 (dc(j,i+nres),j=1,3)
1853 adt=d_a_old(j,0)*d_time
1855 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1856 d_t_new(j,0)=d_t_old(j,0)+adt2
1857 d_t(j,0)=d_t_old(j,0)+adt
1861 adt=d_a_old(j,i)*d_time
1863 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1864 d_t_new(j,i)=d_t_old(j,i)+adt2
1865 d_t(j,i)=d_t_old(j,i)+adt
1870 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
1871 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
1872 .and.(mnum.ne.5)) then
1875 adt=d_a_old(j,inres)*d_time
1877 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1878 d_t_new(j,inres)=d_t_old(j,inres)+adt2
1879 d_t(j,inres)=d_t_old(j,inres)+adt
1884 write (iout,*) "VELVERLET1 END: DC"
1886 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),&
1887 (dc(j,i+nres),j=1,3)
1891 end subroutine verlet1
1892 !-----------------------------------------------------------------------------
1894 ! Step 2 of the velocity Verlet algorithm: update velocities
1896 ! implicit real*8 (a-h,o-z)
1897 ! include 'DIMENSIONS'
1898 ! include 'COMMON.CONTROL'
1899 ! include 'COMMON.VAR'
1900 ! include 'COMMON.MD'
1901 ! include 'COMMON.CHAIN'
1902 ! include 'COMMON.DERIV'
1903 ! include 'COMMON.GEO'
1904 ! include 'COMMON.LOCAL'
1905 ! include 'COMMON.INTERACT'
1906 ! include 'COMMON.IOUNITS'
1907 ! include 'COMMON.NAMES'
1908 integer :: i,j,inres,mnum
1911 d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
1915 d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
1920 ! iti=iabs(itype(i,mnum))
1921 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
1922 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
1923 .and.(mnum.ne.5)) then
1926 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
1931 end subroutine verlet2
1932 !-----------------------------------------------------------------------------
1933 subroutine sddir_precalc
1934 ! Applying velocity Verlet algorithm - step 1 to coordinates
1935 ! implicit real*8 (a-h,o-z)
1936 ! include 'DIMENSIONS'
1942 ! include 'COMMON.CONTROL'
1943 ! include 'COMMON.VAR'
1944 ! include 'COMMON.MD'
1946 ! include 'COMMON.LANGEVIN'
1948 ! include 'COMMON.LANGEVIN.lang0'
1950 ! include 'COMMON.CHAIN'
1951 ! include 'COMMON.DERIV'
1952 ! include 'COMMON.GEO'
1953 ! include 'COMMON.LOCAL'
1954 ! include 'COMMON.INTERACT'
1955 ! include 'COMMON.IOUNITS'
1956 ! include 'COMMON.NAMES'
1957 ! include 'COMMON.TIME1'
1958 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1959 !el common /stochcalc/ stochforcvec
1960 real(kind=8) :: time00
1962 ! Compute friction and stochastic forces
1967 time_fric=time_fric+MPI_Wtime()-time00
1969 call stochastic_force(stochforcvec)
1970 time_stoch=time_stoch+MPI_Wtime()-time00
1973 ! Compute the acceleration due to friction forces (d_af_work) and stochastic
1974 ! forces (d_as_work)
1976 call ginv_mult(fric_work, d_af_work)
1977 call ginv_mult(stochforcvec, d_as_work)
1979 end subroutine sddir_precalc
1980 !-----------------------------------------------------------------------------
1981 subroutine sddir_verlet1
1982 ! Applying velocity Verlet algorithm - step 1 to velocities
1985 ! implicit real*8 (a-h,o-z)
1986 ! include 'DIMENSIONS'
1987 ! include 'COMMON.CONTROL'
1988 ! include 'COMMON.VAR'
1989 ! include 'COMMON.MD'
1991 ! include 'COMMON.LANGEVIN'
1993 ! include 'COMMON.LANGEVIN.lang0'
1995 ! include 'COMMON.CHAIN'
1996 ! include 'COMMON.DERIV'
1997 ! include 'COMMON.GEO'
1998 ! include 'COMMON.LOCAL'
1999 ! include 'COMMON.INTERACT'
2000 ! include 'COMMON.IOUNITS'
2001 ! include 'COMMON.NAMES'
2002 ! Revised 3/31/05 AL: correlation between random contributions to
2003 ! position and velocity increments included.
2004 real(kind=8) :: sqrt13 = 0.57735026918962576451d0 ! 1/sqrt(3)
2005 real(kind=8) :: adt,adt2
2006 integer :: i,j,ind,inres,mnum
2008 ! Add the contribution from BOTH friction and stochastic force to the
2009 ! coordinates, but ONLY the contribution from the friction forces to velocities
2012 adt=(d_a_old(j,0)+d_af_work(j))*d_time
2013 adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
2014 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
2015 d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
2016 d_t(j,0)=d_t_old(j,0)+adt
2021 adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
2022 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
2023 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
2024 d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
2025 d_t(j,i)=d_t_old(j,i)+adt
2031 ! iti=iabs(itype(i,mnum))
2032 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2033 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2034 .and.(mnum.ne.5)) then
2037 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
2038 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
2039 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
2040 d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
2041 d_t(j,inres)=d_t_old(j,inres)+adt
2047 end subroutine sddir_verlet1
2048 !-----------------------------------------------------------------------------
2049 subroutine sddir_verlet2
2050 ! Calculating the adjusted velocities for accelerations
2053 ! implicit real*8 (a-h,o-z)
2054 ! include 'DIMENSIONS'
2055 ! include 'COMMON.CONTROL'
2056 ! include 'COMMON.VAR'
2057 ! include 'COMMON.MD'
2059 ! include 'COMMON.LANGEVIN'
2061 ! include 'COMMON.LANGEVIN.lang0'
2063 ! include 'COMMON.CHAIN'
2064 ! include 'COMMON.DERIV'
2065 ! include 'COMMON.GEO'
2066 ! include 'COMMON.LOCAL'
2067 ! include 'COMMON.INTERACT'
2068 ! include 'COMMON.IOUNITS'
2069 ! include 'COMMON.NAMES'
2070 real(kind=8),dimension(6*nres) :: stochforcvec,d_as_work1 !(MAXRES6) maxres6=6*maxres
2071 real(kind=8) :: cos60 = 0.5d0, sin60 = 0.86602540378443864676d0
2072 integer :: i,j,ind,inres,mnum
2073 ! Revised 3/31/05 AL: correlation between random contributions to
2074 ! position and velocity increments included.
2075 ! The correlation coefficients are calculated at low-friction limit.
2076 ! Also, friction forces are now not calculated with new velocities.
2078 ! call friction_force
2079 call stochastic_force(stochforcvec)
2081 ! Compute the acceleration due to friction forces (d_af_work) and stochastic
2082 ! forces (d_as_work)
2084 call ginv_mult(stochforcvec, d_as_work1)
2090 d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j)) &
2091 +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
2096 d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j)) &
2097 +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
2103 ! iti=iabs(itype(i,mnum))
2104 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2105 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2106 .and.(mnum.ne.5)) then
2109 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres) &
2110 +d_af_work(ind+j))+sin60*d_as_work(ind+j) &
2111 +cos60*d_as_work1(ind+j))*d_time
2117 end subroutine sddir_verlet2
2118 !-----------------------------------------------------------------------------
2119 subroutine max_accel
2121 ! Find the maximum difference in the accelerations of the the sites
2122 ! at the beginning and the end of the time step.
2125 ! implicit real*8 (a-h,o-z)
2126 ! include 'DIMENSIONS'
2127 ! include 'COMMON.CONTROL'
2128 ! include 'COMMON.VAR'
2129 ! include 'COMMON.MD'
2130 ! include 'COMMON.CHAIN'
2131 ! include 'COMMON.DERIV'
2132 ! include 'COMMON.GEO'
2133 ! include 'COMMON.LOCAL'
2134 ! include 'COMMON.INTERACT'
2135 ! include 'COMMON.IOUNITS'
2136 real(kind=8),dimension(3) :: aux,accel,accel_old
2137 real(kind=8) :: dacc
2141 ! aux(j)=d_a(j,0)-d_a_old(j,0)
2142 accel_old(j)=d_a_old(j,0)
2149 ! 7/3/08 changed to asymmetric difference
2151 ! accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
2152 accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
2153 accel(j)=accel(j)+0.5d0*d_a(j,i)
2154 ! if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
2155 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
2156 dacc=dabs(accel(j)-accel_old(j))
2157 ! write (iout,*) i,dacc
2158 if (dacc.gt.amax) amax=dacc
2166 accel_old(j)=d_a_old(j,0)
2171 accel_old(j)=accel_old(j)+d_a_old(j,1)
2172 accel(j)=accel(j)+d_a(j,1)
2177 ! iti=iabs(itype(i,mnum))
2178 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2179 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2180 .and.(mnum.ne.5)) then
2182 ! accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
2183 accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
2184 accel(j)=accel(j)+d_a(j,i+nres)
2188 ! if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
2189 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
2190 dacc=dabs(accel(j)-accel_old(j))
2191 ! write (iout,*) "side-chain",i,dacc
2192 if (dacc.gt.amax) amax=dacc
2196 accel_old(j)=accel_old(j)+d_a_old(j,i)
2197 accel(j)=accel(j)+d_a(j,i)
2198 ! aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
2202 end subroutine max_accel
2203 !-----------------------------------------------------------------------------
2204 subroutine predict_edrift(epdrift)
2206 ! Predict the drift of the potential energy
2209 use control_data, only: lmuca
2210 ! implicit real*8 (a-h,o-z)
2211 ! include 'DIMENSIONS'
2212 ! include 'COMMON.CONTROL'
2213 ! include 'COMMON.VAR'
2214 ! include 'COMMON.MD'
2215 ! include 'COMMON.CHAIN'
2216 ! include 'COMMON.DERIV'
2217 ! include 'COMMON.GEO'
2218 ! include 'COMMON.LOCAL'
2219 ! include 'COMMON.INTERACT'
2220 ! include 'COMMON.IOUNITS'
2221 ! include 'COMMON.MUCA'
2222 real(kind=8) :: epdrift,epdriftij
2224 ! Drift of the potential energy
2230 epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
2231 if (lmuca) epdriftij=epdriftij*factor
2232 ! write (iout,*) "back",i,j,epdriftij
2233 if (epdriftij.gt.epdrift) epdrift=epdriftij
2237 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1.and.&
2238 molnum(i).ne.5) then
2241 dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
2242 if (lmuca) epdriftij=epdriftij*factor
2243 ! write (iout,*) "side",i,j,epdriftij
2244 if (epdriftij.gt.epdrift) epdrift=epdriftij
2248 epdrift=0.5d0*epdrift*d_time*d_time
2249 ! write (iout,*) "epdrift",epdrift
2251 end subroutine predict_edrift
2252 !-----------------------------------------------------------------------------
2253 subroutine verlet_bath
2255 ! Coupling to the thermostat by using the Berendsen algorithm
2258 ! implicit real*8 (a-h,o-z)
2259 ! include 'DIMENSIONS'
2260 ! include 'COMMON.CONTROL'
2261 ! include 'COMMON.VAR'
2262 ! include 'COMMON.MD'
2263 ! include 'COMMON.CHAIN'
2264 ! include 'COMMON.DERIV'
2265 ! include 'COMMON.GEO'
2266 ! include 'COMMON.LOCAL'
2267 ! include 'COMMON.INTERACT'
2268 ! include 'COMMON.IOUNITS'
2269 ! include 'COMMON.NAMES'
2270 real(kind=8) :: T_half,fact
2271 integer :: i,j,inres,mnum
2273 T_half=2.0d0/(dimen3*Rb)*EK
2274 fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
2275 ! write(iout,*) "T_half", T_half
2276 ! write(iout,*) "EK", EK
2277 ! write(iout,*) "fact", fact
2279 d_t(j,0)=fact*d_t(j,0)
2283 d_t(j,i)=fact*d_t(j,i)
2288 ! iti=iabs(itype(i,mnum))
2289 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2290 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2291 .and.(mnum.ne.5)) then
2294 d_t(j,inres)=fact*d_t(j,inres)
2299 end subroutine verlet_bath
2300 !-----------------------------------------------------------------------------
2302 ! Set up the initial conditions of a MD simulation
2305 use control, only:tcpu
2306 !el use io_basic, only:ilen
2309 use minimm, only:minim_dc,minimize,sc_move
2310 use io_config, only:readrst
2311 use io, only:statout
2312 ! implicit real*8 (a-h,o-z)
2313 ! include 'DIMENSIONS'
2316 character(len=16) :: form
2317 integer :: IERROR,ERRCODE
2319 ! include 'COMMON.SETUP'
2320 ! include 'COMMON.CONTROL'
2321 ! include 'COMMON.VAR'
2322 ! include 'COMMON.MD'
2324 ! include 'COMMON.LANGEVIN'
2326 ! include 'COMMON.LANGEVIN.lang0'
2328 ! include 'COMMON.CHAIN'
2329 ! include 'COMMON.DERIV'
2330 ! include 'COMMON.GEO'
2331 ! include 'COMMON.LOCAL'
2332 ! include 'COMMON.INTERACT'
2333 ! include 'COMMON.IOUNITS'
2334 ! include 'COMMON.NAMES'
2335 ! include 'COMMON.REMD'
2336 real(kind=8),dimension(0:n_ene) :: energia_long,energia_short
2337 real(kind=8),dimension(3) :: vcm,incr,L
2338 real(kind=8) :: xv,sigv,lowb,highb
2339 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
2340 character(len=256) :: qstr
2343 character(len=50) :: tytul
2344 logical :: file_exist
2345 !el common /gucio/ cm
2346 integer :: i,j,ipos,iq,iw,nft_sc,iretcode,nfun,itime,ierr,mnum
2347 real(kind=8) :: etot,tt0
2351 ! write(iout,*) "d_time", d_time
2352 ! Compute the standard deviations of stochastic forces for Langevin dynamics
2353 ! if the friction coefficients do not depend on surface area
2354 if (lang.gt.0 .and. .not.surfarea) then
2357 stdforcp(i)=stdfp(mnum)*dsqrt(gamp(mnum))
2361 stdforcsc(i)=stdfsc(iabs(itype(i,mnum)),mnum) &
2362 *dsqrt(gamsc(iabs(itype(i,mnum)),mnum))
2365 ! Open the pdb file for snapshotshots
2368 if (ilen(tmpdir).gt.0) &
2369 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2370 liczba(:ilen(liczba))//".pdb")
2372 file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2376 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) &
2377 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2378 liczba(:ilen(liczba))//".x")
2379 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2382 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) &
2383 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2384 liczba(:ilen(liczba))//".cx")
2385 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2391 if (ilen(tmpdir).gt.0) &
2392 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
2393 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
2395 if (ilen(tmpdir).gt.0) &
2396 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
2397 cartname=prefix(:ilen(prefix))//"_MD.cx"
2401 write (qstr,'(256(1h ))')
2404 iq = qinfrag(i,iset)*10
2405 iw = wfrag(i,iset)/100
2407 if(me.eq.king.or..not.out1file) &
2408 write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
2409 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
2414 iq = qinpair(i,iset)*10
2415 iw = wpair(i,iset)/100
2417 if(me.eq.king.or..not.out1file) &
2418 write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
2419 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
2423 ! pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
2425 ! cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
2427 ! cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
2429 ! statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
2433 if (restart1file) then
2435 inquire(file=mremd_rst_name,exist=file_exist)
2436 write (*,*) me," Before broadcast: file_exist",file_exist
2438 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,&
2441 write (*,*) me," After broadcast: file_exist",file_exist
2442 ! inquire(file=mremd_rst_name,exist=file_exist)
2443 if(me.eq.king.or..not.out1file) &
2444 write(iout,*) "Initial state read by master and distributed"
2446 if (ilen(tmpdir).gt.0) &
2447 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_' &
2448 //liczba(:ilen(liczba))//'.rst')
2449 inquire(file=rest2name,exist=file_exist)
2452 if(.not.restart1file) then
2453 if(me.eq.king.or..not.out1file) &
2454 write(iout,*) "Initial state will be read from file ",&
2455 rest2name(:ilen(rest2name))
2458 call rescale_weights(t_bath)
2460 if(me.eq.king.or..not.out1file)then
2461 if (restart1file) then
2462 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),&
2465 write(iout,*) "File ",rest2name(:ilen(rest2name)),&
2468 write(iout,*) "Initial velocities randomly generated"
2475 ! Generate initial velocities
2476 if(me.eq.king.or..not.out1file) &
2477 write(iout,*) "Initial velocities randomly generated"
2482 ! rest2name = prefix(:ilen(prefix))//'.rst'
2483 if(me.eq.king.or..not.out1file)then
2484 write (iout,*) "Initial velocities"
2486 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2487 (d_t(j,i+nres),j=1,3)
2489 ! Zeroing the total angular momentum of the system
2490 write(iout,*) "Calling the zero-angular momentum subroutine"
2493 ! Getting the potential energy and forces and velocities and accelerations
2495 ! write (iout,*) "velocity of the center of the mass:"
2496 ! write (iout,*) (vcm(j),j=1,3)
2498 d_t(j,0)=d_t(j,0)-vcm(j)
2500 ! Removing the velocity of the center of mass
2502 if(me.eq.king.or..not.out1file)then
2503 write (iout,*) "vcm right after adjustment:"
2504 write (iout,*) (vcm(j),j=1,3)
2506 if ((.not.rest).and.(indpdb.eq.0)) then
2508 if(iranconf.ne.0) then
2510 print *, 'Calling OVERLAP_SC'
2511 call overlap_sc(fail)
2514 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
2515 print *,'SC_move',nft_sc,etot
2516 if(me.eq.king.or..not.out1file) &
2517 write(iout,*) 'SC_move',nft_sc,etot
2521 print *, 'Calling MINIM_DC'
2522 call minim_dc(etot,iretcode,nfun)
2524 call geom_to_var(nvar,varia)
2525 print *,'Calling MINIMIZE.'
2526 call minimize(etot,varia,iretcode,nfun)
2527 call var_to_geom(nvar,varia)
2529 if(me.eq.king.or..not.out1file) &
2530 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2533 call chainbuild_cart
2538 kinetic_T=2.0d0/(dimen3*Rb)*EK
2539 if(me.eq.king.or..not.out1file)then
2549 call etotal(potEcomp)
2550 if (large) call enerprint(potEcomp)
2553 t_etotal=t_etotal+MPI_Wtime()-tt0
2555 t_etotal=t_etotal+tcpu()-tt0
2562 if (amax*d_time .gt. dvmax) then
2563 d_time=d_time*dvmax/amax
2564 if(me.eq.king.or..not.out1file) write (iout,*) &
2565 "Time step reduced to",d_time,&
2566 " because of too large initial acceleration."
2568 if(me.eq.king.or..not.out1file)then
2569 write(iout,*) "Potential energy and its components"
2570 call enerprint(potEcomp)
2571 ! write(iout,*) (potEcomp(i),i=0,n_ene)
2573 potE=potEcomp(0)-potEcomp(20)
2576 if (ntwe.ne.0) call statout(itime)
2577 if(me.eq.king.or..not.out1file) &
2578 write (iout,'(/a/3(a25,1pe14.5/))') "Initial:", &
2579 " Kinetic energy",EK," Potential energy",potE, &
2580 " Total energy",totE," Maximum acceleration ", &
2583 write (iout,*) "Initial coordinates"
2585 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),&
2588 write (iout,*) "Initial dC"
2590 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),&
2591 (dc(j,i+nres),j=1,3)
2593 write (iout,*) "Initial velocities"
2594 write (iout,"(13x,' backbone ',23x,' side chain')")
2596 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2597 (d_t(j,i+nres),j=1,3)
2599 write (iout,*) "Initial accelerations"
2601 ! write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2602 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),&
2603 (d_a(j,i+nres),j=1,3)
2609 d_t_old(j,i)=d_t(j,i)
2610 d_a_old(j,i)=d_a(j,i)
2612 ! write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
2621 call etotal_short(energia_short)
2622 if (large) call enerprint(potEcomp)
2625 t_eshort=t_eshort+MPI_Wtime()-tt0
2627 t_eshort=t_eshort+tcpu()-tt0
2632 if(.not.out1file .and. large) then
2633 write (iout,*) "energia_long",energia_long(0),&
2634 " energia_short",energia_short(0),&
2635 " total",energia_long(0)+energia_short(0)
2636 write (iout,*) "Initial fast-force accelerations"
2638 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
2639 (d_a(j,i+nres),j=1,3)
2642 ! 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
2645 d_a_short(j,i)=d_a(j,i)
2654 call etotal_long(energia_long)
2655 if (large) call enerprint(potEcomp)
2658 t_elong=t_elong+MPI_Wtime()-tt0
2660 t_elong=t_elong+tcpu()-tt0
2665 if(.not.out1file .and. large) then
2666 write (iout,*) "energia_long",energia_long(0)
2667 write (iout,*) "Initial slow-force accelerations"
2669 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
2670 (d_a(j,i+nres),j=1,3)
2674 t_enegrad=t_enegrad+MPI_Wtime()-tt0
2676 t_enegrad=t_enegrad+tcpu()-tt0
2680 end subroutine init_MD
2681 !-----------------------------------------------------------------------------
2682 subroutine random_vel
2684 ! implicit real*8 (a-h,o-z)
2686 use random, only:anorm_distr
2688 ! include 'DIMENSIONS'
2689 ! include 'COMMON.CONTROL'
2690 ! include 'COMMON.VAR'
2691 ! include 'COMMON.MD'
2693 ! include 'COMMON.LANGEVIN'
2695 ! include 'COMMON.LANGEVIN.lang0'
2697 ! include 'COMMON.CHAIN'
2698 ! include 'COMMON.DERIV'
2699 ! include 'COMMON.GEO'
2700 ! include 'COMMON.LOCAL'
2701 ! include 'COMMON.INTERACT'
2702 ! include 'COMMON.IOUNITS'
2703 ! include 'COMMON.NAMES'
2704 ! include 'COMMON.TIME1'
2705 real(kind=8) :: xv,sigv,lowb,highb ,Ek1
2707 real(kind=8) ,allocatable, dimension(:) :: DDU1,DDU2,DL2,DL1,xsolv,DML,rs
2708 real(kind=8) :: sumx
2710 real(kind=8) ,allocatable, dimension(:) :: rsold
2711 real (kind=8),allocatable,dimension(:,:) :: matold
2714 integer :: i,j,ii,k,ind,mark,imark,mnum
2715 ! Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
2716 ! First generate velocities in the eigenspace of the G matrix
2717 ! write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
2724 sigv=dsqrt((Rb*t_bath)/geigen(i))
2727 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
2729 write (iout,*) "i",i," ii",ii," geigen",geigen(i),&
2730 " d_t_work_new",d_t_work_new(ii)
2741 Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
2744 write (iout,*) "Ek from eigenvectors",Ek1
2745 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
2749 ! Transform velocities to UNRES coordinate space
2750 allocate (DL1(2*nres))
2751 allocate (DDU1(2*nres))
2752 allocate (DL2(2*nres))
2753 allocate (DDU2(2*nres))
2754 allocate (xsolv(2*nres))
2755 allocate (DML(2*nres))
2756 allocate (rs(2*nres))
2758 allocate (rsold(2*nres))
2759 allocate (matold(2*nres,2*nres))
2761 matold(1,1)=DMorig(1)
2762 matold(1,2)=DU1orig(1)
2763 matold(1,3)=DU2orig(1)
2764 write (*,*) DMorig(1),DU1orig(1),DU2orig(1)
2769 matold(i,j)=DMorig(i)
2770 matold(i,j-1)=DU1orig(i-1)
2772 matold(i,j-2)=DU2orig(i-2)
2780 matold(i,j+1)=DU1orig(i)
2786 matold(i,j+2)=DU2orig(i)
2790 matold(dimen,dimen)=DMorig(dimen)
2791 matold(dimen,dimen-1)=DU1orig(dimen-1)
2792 matold(dimen,dimen-2)=DU2orig(dimen-2)
2793 write(iout,*) "old gmatrix"
2794 call matout(dimen,dimen,2*nres,2*nres,matold)
2798 ! Find the ith eigenvector of the pentadiagonal inertiq matrix
2802 DML(j)=DMorig(j)-geigen(i)
2805 DML(j-1)=DMorig(j)-geigen(i)
2810 DDU1(imark-1)=DU2orig(imark-1)
2811 do j=imark+1,dimen-1
2812 DDU1(j-1)=DU1orig(j)
2820 DDU2(j)=DU2orig(j+1)
2829 write (iout,*) "DL2,DL1,DML,DDU1,DDU2"
2830 write(iout,'(10f10.5)') (DL2(k),k=3,dimen-1)
2831 write(iout,'(10f10.5)') (DL1(k),k=2,dimen-1)
2832 write(iout,'(10f10.5)') (DML(k),k=1,dimen-1)
2833 write(iout,'(10f10.5)') (DDU1(k),k=1,dimen-2)
2834 write(iout,'(10f10.5)') (DDU2(k),k=1,dimen-3)
2837 if (imark.gt.2) rs(imark-2)=-DU2orig(imark-2)
2838 if (imark.gt.1) rs(imark-1)=-DU1orig(imark-1)
2839 if (imark.lt.dimen) rs(imark)=-DU1orig(imark)
2840 if (imark.lt.dimen-1) rs(imark+1)=-DU2orig(imark)
2844 ! write (iout,*) "Vector rs"
2846 ! write (iout,*) j,rs(j)
2849 call FDIAG(dimen-1,DL2,DL1,DML,DDU1,DDU2,rs,xsolv,mark)
2856 sumx=-geigen(i)*xsolv(j)
2858 sumx=sumx+matold(j,k)*xsolv(k)
2861 sumx=sumx+matold(j,k)*xsolv(k-1)
2863 write(iout,'(i5,3f10.5)') j,sumx,rsold(j),sumx-rsold(j)
2866 sumx=-geigen(i)*xsolv(j-1)
2868 sumx=sumx+matold(j,k)*xsolv(k)
2871 sumx=sumx+matold(j,k)*xsolv(k-1)
2873 write(iout,'(i5,3f10.5)') j-1,sumx,rsold(j-1),sumx-rsold(j-1)
2877 "Solution of equations system",i," for eigenvalue",geigen(i)
2879 write(iout,'(i5,f10.5)') j,xsolv(j)
2882 do j=dimen-1,imark,-1
2887 write (iout,*) "Un-normalized eigenvector",i," for eigenvalue",geigen(i)
2889 write(iout,'(i5,f10.5)') j,xsolv(j)
2892 ! Normalize ith eigenvector
2895 sumx=sumx+xsolv(j)**2
2899 xsolv(j)=xsolv(j)/sumx
2902 write (iout,*) "Eigenvector",i," for eigenvalue",geigen(i)
2904 write(iout,'(i5,3f10.5)') j,xsolv(j)
2907 ! All done at this point for eigenvector i, exit loop
2915 write (iout,*) "Unable to find eigenvector",i
2918 ! write (iout,*) "i=",i
2920 ! write (iout,*) "k=",k
2923 ! write(iout,*) "ind",ind," ind1",3*(i-1)+k
2924 d_t_work(ind)=d_t_work(ind) &
2925 +xsolv(j)*d_t_work_new(3*(i-1)+k)
2928 enddo ! i (loop over the eigenvectors)
2931 write (iout,*) "d_t_work"
2933 write (iout,"(i5,f10.5)") i,d_t_work(i)
2938 ! if (itype(i,1).eq.10) then
2940 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
2941 iti=iabs(itype(i,mnum))
2942 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2943 if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
2944 .or.(mnum.eq.5)) then
2951 ! write (iout,*) "k",k," ii+k",ii+k," ii+j+k",ii+j+k,"EK1",Ek1
2952 Ek1=Ek1+0.5d0*mp(mnum)*((d_t_work(ii+j+k)+d_t_work_new(ii+k))/2)**2+&
2953 0.5d0*0.25d0*IP(mnum)*(d_t_work(ii+j+k)-d_t_work_new(ii+k))**2
2956 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2957 .and.(mnum.ne.5)) ii=ii+3
2958 write (iout,*) "i",i," itype",itype(i,mnum)," mass",msc(itype(i,mnum),mnum)
2959 write (iout,*) "ii",ii
2962 write (iout,*) "k",k," ii",ii,"EK1",EK1
2963 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2965 Ek1=Ek1+0.5d0*Isc(iabs(itype(i,mnum),mnum))*(d_t_work(ii)-d_t_work(ii-3))**2
2966 Ek1=Ek1+0.5d0*msc(iabs(itype(i,mnum)),mnum)*d_t_work(ii)**2
2968 write (iout,*) "i",i," ii",ii
2970 write (iout,*) "Ek from d_t_work",Ek1
2971 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
2987 d_t(k,j)=d_t_work(ind)
2991 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2992 .and.(mnum.ne.5)) then
2994 d_t(k,j+nres)=d_t_work(ind)
3000 write (iout,*) "Random velocities in the Calpha,SC space"
3002 write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3)
3005 write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3)
3012 ! if (itype(i,1).eq.10) then
3014 if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
3015 .or.(mnum.eq.5)) then
3017 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
3021 d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
3022 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
3027 write (iout,*)"Random velocities in the virtual-bond-vector space"
3029 write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3)
3032 write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3)
3035 write (iout,*) "Ek from d_t_work",Ek1
3036 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
3044 d_t_work(ind)=d_t_work(ind) &
3045 +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
3047 ! write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
3051 ! Transfer to the d_t vector
3053 d_t(j,0)=d_t_work(j)
3059 d_t(j,i)=d_t_work(ind)
3064 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3065 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3066 .and.(mnum.ne.5)) then
3069 d_t(j,i+nres)=d_t_work(ind)
3075 ! write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",&
3076 ! 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
3080 end subroutine random_vel
3081 !-----------------------------------------------------------------------------
3083 subroutine sd_verlet_p_setup
3084 ! Sets up the parameters of stochastic Verlet algorithm
3085 ! implicit real*8 (a-h,o-z)
3086 ! include 'DIMENSIONS'
3087 use control, only: tcpu
3092 ! include 'COMMON.CONTROL'
3093 ! include 'COMMON.VAR'
3094 ! include 'COMMON.MD'
3096 ! include 'COMMON.LANGEVIN'
3098 ! include 'COMMON.LANGEVIN.lang0'
3100 ! include 'COMMON.CHAIN'
3101 ! include 'COMMON.DERIV'
3102 ! include 'COMMON.GEO'
3103 ! include 'COMMON.LOCAL'
3104 ! include 'COMMON.INTERACT'
3105 ! include 'COMMON.IOUNITS'
3106 ! include 'COMMON.NAMES'
3107 ! include 'COMMON.TIME1'
3108 real(kind=8),dimension(6*nres) :: emgdt !(MAXRES6) maxres6=6*maxres
3109 real(kind=8) :: pterm,vterm,rho,rhoc,vsig
3110 real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
3111 prand_vec,vrand_vec1,vrand_vec2 !(MAXRES6) maxres6=6*maxres
3112 logical :: lprn = .false.
3113 real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
3114 real(kind=8) :: ktm,gdt,egdt,gdt2,gdt3,gdt4,gdt5,gdt6,gdt7,gdt8,&
3116 integer :: i,maxres2
3123 ! AL 8/17/04 Code adapted from tinker
3125 ! Get the frictional and random terms for stochastic dynamics in the
3126 ! eigenspace of mass-scaled UNRES friction matrix
3130 gdt = fricgam(i) * d_time
3132 ! Stochastic dynamics reduces to simple MD for zero friction
3134 if (gdt .le. zero) then
3135 pfric_vec(i) = 1.0d0
3136 vfric_vec(i) = d_time
3137 afric_vec(i) = 0.5d0 * d_time * d_time
3138 prand_vec(i) = 0.0d0
3139 vrand_vec1(i) = 0.0d0
3140 vrand_vec2(i) = 0.0d0
3142 ! Analytical expressions when friction coefficient is large
3145 if (gdt .ge. gdt_radius) then
3148 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
3149 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
3150 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
3151 vterm = 1.0d0 - egdt**2
3152 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
3154 ! Use series expansions when friction coefficient is small
3165 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0 &
3166 - gdt5/120.0d0 + gdt6/720.0d0 &
3167 - gdt7/5040.0d0 + gdt8/40320.0d0 &
3168 - gdt9/362880.0d0) / fricgam(i)**2
3169 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
3170 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
3171 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0 &
3172 + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0 &
3173 + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0 &
3174 + 127.0d0*gdt9/90720.0d0
3175 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0 &
3176 - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0 &
3177 - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0 &
3178 - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
3179 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0 &
3180 - 17.0d0*gdt2/1280.0d0 &
3181 + 17.0d0*gdt3/6144.0d0 &
3182 + 40967.0d0*gdt4/34406400.0d0 &
3183 - 57203.0d0*gdt5/275251200.0d0 &
3184 - 1429487.0d0*gdt6/13212057600.0d0)
3187 ! Compute the scaling factors of random terms for the nonzero friction case
3189 ktm = 0.5d0*d_time/fricgam(i)
3190 psig = dsqrt(ktm*pterm) / fricgam(i)
3191 vsig = dsqrt(ktm*vterm)
3192 rhoc = dsqrt(1.0d0 - rho*rho)
3194 vrand_vec1(i) = vsig * rho
3195 vrand_vec2(i) = vsig * rhoc
3200 "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",&
3203 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),&
3204 afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
3208 ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
3211 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
3212 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
3213 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
3214 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
3215 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
3216 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
3219 t_sdsetup=t_sdsetup+MPI_Wtime()
3221 t_sdsetup=t_sdsetup+tcpu()-tt0
3224 end subroutine sd_verlet_p_setup
3225 !-----------------------------------------------------------------------------
3226 subroutine eigtransf1(n,ndim,ab,d,c)
3230 real(kind=8) :: ab(ndim,ndim,n),c(ndim,n),d(ndim)
3236 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
3241 end subroutine eigtransf1
3242 !-----------------------------------------------------------------------------
3243 subroutine eigtransf(n,ndim,a,b,d,c)
3247 real(kind=8) :: a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
3253 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
3258 end subroutine eigtransf
3259 !-----------------------------------------------------------------------------
3260 subroutine sd_verlet1
3262 ! Applying stochastic velocity Verlet algorithm - step 1 to velocities
3264 ! implicit real*8 (a-h,o-z)
3265 ! include 'DIMENSIONS'
3266 ! include 'COMMON.CONTROL'
3267 ! include 'COMMON.VAR'
3268 ! include 'COMMON.MD'
3270 ! include 'COMMON.LANGEVIN'
3272 ! include 'COMMON.LANGEVIN.lang0'
3274 ! include 'COMMON.CHAIN'
3275 ! include 'COMMON.DERIV'
3276 ! include 'COMMON.GEO'
3277 ! include 'COMMON.LOCAL'
3278 ! include 'COMMON.INTERACT'
3279 ! include 'COMMON.IOUNITS'
3280 ! include 'COMMON.NAMES'
3281 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
3282 !el common /stochcalc/ stochforcvec
3283 logical :: lprn = .false.
3284 real(kind=8) :: ddt1,ddt2
3285 integer :: i,j,ind,inres
3287 ! write (iout,*) "dc_old"
3289 ! write (iout,'(i5,3f10.5,5x,3f10.5)')
3290 ! & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3293 dc_work(j)=dc_old(j,0)
3294 d_t_work(j)=d_t_old(j,0)
3295 d_a_work(j)=d_a_old(j,0)
3300 dc_work(ind+j)=dc_old(j,i)
3301 d_t_work(ind+j)=d_t_old(j,i)
3302 d_a_work(ind+j)=d_a_old(j,i)
3308 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3309 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3310 .and.(mnum.ne.5)) then
3312 dc_work(ind+j)=dc_old(j,i+nres)
3313 d_t_work(ind+j)=d_t_old(j,i+nres)
3314 d_a_work(ind+j)=d_a_old(j,i+nres)
3322 "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",&
3326 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
3327 vfric_mat(i,j),afric_mat(i,j),&
3328 prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3336 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) &
3337 +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3338 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3339 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3341 d_t_work_new(i)=ddt1+0.5d0*ddt2
3342 d_t_work(i)=ddt1+ddt2
3347 d_t(j,0)=d_t_work(j)
3352 dc(j,i)=dc_work(ind+j)
3353 d_t(j,i)=d_t_work(ind+j)
3359 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3360 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3361 .and.(mnum.ne.5)) then
3364 dc(j,inres)=dc_work(ind+j)
3365 d_t(j,inres)=d_t_work(ind+j)
3371 end subroutine sd_verlet1
3372 !-----------------------------------------------------------------------------
3373 subroutine sd_verlet2
3375 ! Calculating the adjusted velocities for accelerations
3377 ! implicit real*8 (a-h,o-z)
3378 ! include 'DIMENSIONS'
3379 ! include 'COMMON.CONTROL'
3380 ! include 'COMMON.VAR'
3381 ! include 'COMMON.MD'
3383 ! include 'COMMON.LANGEVIN'
3385 ! include 'COMMON.LANGEVIN.lang0'
3387 ! include 'COMMON.CHAIN'
3388 ! include 'COMMON.DERIV'
3389 ! include 'COMMON.GEO'
3390 ! include 'COMMON.LOCAL'
3391 ! include 'COMMON.INTERACT'
3392 ! include 'COMMON.IOUNITS'
3393 ! include 'COMMON.NAMES'
3394 !el real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV !(MAXRES6) maxres6=6*maxres
3395 real(kind=8),dimension(6*nres) :: stochforcvecV !(MAXRES6) maxres6=6*maxres
3396 !el common /stochcalc/ stochforcvec
3398 real(kind=8) :: ddt1,ddt2
3399 integer :: i,j,ind,inres
3400 ! Compute the stochastic forces which contribute to velocity change
3402 call stochastic_force(stochforcvecV)
3409 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3410 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+ &
3411 vrand_mat2(i,j)*stochforcvecV(j)
3413 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
3417 d_t(j,0)=d_t_work(j)
3422 d_t(j,i)=d_t_work(ind+j)
3428 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3429 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3430 .and.(mnum.ne.5)) then
3433 d_t(j,inres)=d_t_work(ind+j)
3439 end subroutine sd_verlet2
3440 !-----------------------------------------------------------------------------
3441 subroutine sd_verlet_ciccotti_setup
3443 ! Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
3445 ! implicit real*8 (a-h,o-z)
3446 ! include 'DIMENSIONS'
3447 use control, only: tcpu
3452 ! include 'COMMON.CONTROL'
3453 ! include 'COMMON.VAR'
3454 ! include 'COMMON.MD'
3456 ! include 'COMMON.LANGEVIN'
3458 ! include 'COMMON.LANGEVIN.lang0'
3460 ! include 'COMMON.CHAIN'
3461 ! include 'COMMON.DERIV'
3462 ! include 'COMMON.GEO'
3463 ! include 'COMMON.LOCAL'
3464 ! include 'COMMON.INTERACT'
3465 ! include 'COMMON.IOUNITS'
3466 ! include 'COMMON.NAMES'
3467 ! include 'COMMON.TIME1'
3468 real(kind=8),dimension(6*nres) :: emgdt !(MAXRES6) maxres6=6*maxres
3469 real(kind=8) :: pterm,vterm,rho,rhoc,vsig
3470 real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
3471 prand_vec,vrand_vec1,vrand_vec2 !(MAXRES6) maxres6=6*maxres
3472 logical :: lprn = .false.
3473 real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
3474 real(kind=8) :: ktm,gdt,egdt,tt0
3475 integer :: i,maxres2
3482 ! AL 8/17/04 Code adapted from tinker
3484 ! Get the frictional and random terms for stochastic dynamics in the
3485 ! eigenspace of mass-scaled UNRES friction matrix
3489 write (iout,*) "i",i," fricgam",fricgam(i)
3490 gdt = fricgam(i) * d_time
3492 ! Stochastic dynamics reduces to simple MD for zero friction
3494 if (gdt .le. zero) then
3495 pfric_vec(i) = 1.0d0
3496 vfric_vec(i) = d_time
3497 afric_vec(i) = 0.5d0*d_time*d_time
3498 prand_vec(i) = afric_vec(i)
3499 vrand_vec2(i) = vfric_vec(i)
3501 ! Analytical expressions when friction coefficient is large
3506 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
3507 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
3508 prand_vec(i) = afric_vec(i)
3509 vrand_vec2(i) = vfric_vec(i)
3511 ! Compute the scaling factors of random terms for the nonzero friction case
3513 ! ktm = 0.5d0*d_time/fricgam(i)
3514 ! psig = dsqrt(ktm*pterm) / fricgam(i)
3515 ! vsig = dsqrt(ktm*vterm)
3516 ! prand_vec(i) = psig*afric_vec(i)
3517 ! vrand_vec2(i) = vsig*vfric_vec(i)
3522 "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",&
3525 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),&
3526 afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
3530 ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
3532 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
3533 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
3534 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
3535 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
3536 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
3538 t_sdsetup=t_sdsetup+MPI_Wtime()
3540 t_sdsetup=t_sdsetup+tcpu()-tt0
3543 end subroutine sd_verlet_ciccotti_setup
3544 !-----------------------------------------------------------------------------
3545 subroutine sd_verlet1_ciccotti
3547 ! Applying stochastic velocity Verlet algorithm - step 1 to velocities
3548 ! implicit real*8 (a-h,o-z)
3550 ! include 'DIMENSIONS'
3554 ! include 'COMMON.CONTROL'
3555 ! include 'COMMON.VAR'
3556 ! include 'COMMON.MD'
3558 ! include 'COMMON.LANGEVIN'
3560 ! include 'COMMON.LANGEVIN.lang0'
3562 ! include 'COMMON.CHAIN'
3563 ! include 'COMMON.DERIV'
3564 ! include 'COMMON.GEO'
3565 ! include 'COMMON.LOCAL'
3566 ! include 'COMMON.INTERACT'
3567 ! include 'COMMON.IOUNITS'
3568 ! include 'COMMON.NAMES'
3569 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
3570 !el common /stochcalc/ stochforcvec
3571 logical :: lprn = .false.
3572 real(kind=8) :: ddt1,ddt2
3573 integer :: i,j,ind,inres
3574 ! write (iout,*) "dc_old"
3576 ! write (iout,'(i5,3f10.5,5x,3f10.5)')
3577 ! & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3580 dc_work(j)=dc_old(j,0)
3581 d_t_work(j)=d_t_old(j,0)
3582 d_a_work(j)=d_a_old(j,0)
3587 dc_work(ind+j)=dc_old(j,i)
3588 d_t_work(ind+j)=d_t_old(j,i)
3589 d_a_work(ind+j)=d_a_old(j,i)
3594 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3596 dc_work(ind+j)=dc_old(j,i+nres)
3597 d_t_work(ind+j)=d_t_old(j,i+nres)
3598 d_a_work(ind+j)=d_a_old(j,i+nres)
3607 "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",&
3611 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
3612 vfric_mat(i,j),afric_mat(i,j),&
3613 prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3621 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) &
3622 +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3623 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3624 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3626 d_t_work_new(i)=ddt1+0.5d0*ddt2
3627 d_t_work(i)=ddt1+ddt2
3632 d_t(j,0)=d_t_work(j)
3637 dc(j,i)=dc_work(ind+j)
3638 d_t(j,i)=d_t_work(ind+j)
3644 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3645 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3646 .and.(mnum.ne.5)) then
3649 dc(j,inres)=dc_work(ind+j)
3650 d_t(j,inres)=d_t_work(ind+j)
3656 end subroutine sd_verlet1_ciccotti
3657 !-----------------------------------------------------------------------------
3658 subroutine sd_verlet2_ciccotti
3660 ! Calculating the adjusted velocities for accelerations
3662 ! implicit real*8 (a-h,o-z)
3663 ! include 'DIMENSIONS'
3664 ! include 'COMMON.CONTROL'
3665 ! include 'COMMON.VAR'
3666 ! include 'COMMON.MD'
3668 ! include 'COMMON.LANGEVIN'
3670 ! include 'COMMON.LANGEVIN.lang0'
3672 ! include 'COMMON.CHAIN'
3673 ! include 'COMMON.DERIV'
3674 ! include 'COMMON.GEO'
3675 ! include 'COMMON.LOCAL'
3676 ! include 'COMMON.INTERACT'
3677 ! include 'COMMON.IOUNITS'
3678 ! include 'COMMON.NAMES'
3679 !el real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV !(MAXRES6) maxres6=6*maxres
3680 real(kind=8),dimension(6*nres) :: stochforcvecV !(MAXRES6) maxres6=6*maxres
3681 !el common /stochcalc/ stochforcvec
3682 real(kind=8) :: ddt1,ddt2
3683 integer :: i,j,ind,inres
3685 ! Compute the stochastic forces which contribute to velocity change
3687 call stochastic_force(stochforcvecV)
3694 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3695 ! ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
3696 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
3698 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
3702 d_t(j,0)=d_t_work(j)
3707 d_t(j,i)=d_t_work(ind+j)
3713 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3715 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3718 d_t(j,inres)=d_t_work(ind+j)
3724 end subroutine sd_verlet2_ciccotti
3726 !-----------------------------------------------------------------------------
3728 !-----------------------------------------------------------------------------
3729 subroutine inertia_tensor
3731 ! Calculating the intertia tensor for the entire protein in order to
3732 ! remove the perpendicular components of velocity matrix which cause
3733 ! the molecule to rotate.
3736 ! implicit real*8 (a-h,o-z)
3737 ! include 'DIMENSIONS'
3738 ! include 'COMMON.CONTROL'
3739 ! include 'COMMON.VAR'
3740 ! include 'COMMON.MD'
3741 ! include 'COMMON.CHAIN'
3742 ! include 'COMMON.DERIV'
3743 ! include 'COMMON.GEO'
3744 ! include 'COMMON.LOCAL'
3745 ! include 'COMMON.INTERACT'
3746 ! include 'COMMON.IOUNITS'
3747 ! include 'COMMON.NAMES'
3749 real(kind=8),dimension(3,3) :: Im,Imcp,eigvec,Id
3750 real(kind=8),dimension(3) :: pr,eigval,L,vp,vrot
3751 real(kind=8) :: M_SC,mag,mag2,M_PEP
3752 real(kind=8),dimension(3,0:nres) :: vpp !(3,0:MAXRES)
3753 real(kind=8),dimension(3) :: vs_p,pp,incr,v
3754 real(kind=8),dimension(3,3) :: pr1,pr2
3756 !el common /gucio/ cm
3757 integer :: iti,inres,i,j,k,mnum
3768 ! calculating the center of the mass of the protein
3772 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
3773 if (itype(i,mnum).eq.ntyp1_molec(mnum)) cycle
3774 M_PEP=M_PEP+mp(mnum)
3776 cm(j)=cm(j)+(c(j,i)+0.5d0*dc(j,i))*mp(mnum)
3785 if (mnum.eq.5) cycle
3786 iti=iabs(itype(i,mnum))
3787 M_SC=M_SC+msc(iabs(iti),mnum)
3790 cm(j)=cm(j)+msc(iabs(iti),mnum)*c(j,inres)
3794 cm(j)=cm(j)/(M_SC+M_PEP)
3799 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
3801 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
3803 Im(1,1)=Im(1,1)+mp(mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
3804 Im(1,2)=Im(1,2)-mp(mnum)*pr(1)*pr(2)
3805 Im(1,3)=Im(1,3)-mp(mnum)*pr(1)*pr(3)
3806 Im(2,3)=Im(2,3)-mp(mnum)*pr(2)*pr(3)
3807 Im(2,2)=Im(2,2)+mp(mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
3808 Im(3,3)=Im(3,3)+mp(mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
3813 iti=iabs(itype(i,mnum))
3814 if (mnum.eq.5) cycle
3817 pr(j)=c(j,inres)-cm(j)
3819 Im(1,1)=Im(1,1)+msc(iabs(iti),mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
3820 Im(1,2)=Im(1,2)-msc(iabs(iti),mnum)*pr(1)*pr(2)
3821 Im(1,3)=Im(1,3)-msc(iabs(iti),mnum)*pr(1)*pr(3)
3822 Im(2,3)=Im(2,3)-msc(iabs(iti),mnum)*pr(2)*pr(3)
3823 Im(2,2)=Im(2,2)+msc(iabs(iti),mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
3824 Im(3,3)=Im(3,3)+msc(iabs(iti),mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
3829 Im(1,1)=Im(1,1)+Ip(mnum)*(1-dc_norm(1,i)*dc_norm(1,i))* &
3830 vbld(i+1)*vbld(i+1)*0.25d0
3831 Im(1,2)=Im(1,2)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(2,i))* &
3832 vbld(i+1)*vbld(i+1)*0.25d0
3833 Im(1,3)=Im(1,3)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(3,i))* &
3834 vbld(i+1)*vbld(i+1)*0.25d0
3835 Im(2,3)=Im(2,3)+Ip(mnum)*(-dc_norm(2,i)*dc_norm(3,i))* &
3836 vbld(i+1)*vbld(i+1)*0.25d0
3837 Im(2,2)=Im(2,2)+Ip(mnum)*(1-dc_norm(2,i)*dc_norm(2,i))* &
3838 vbld(i+1)*vbld(i+1)*0.25d0
3839 Im(3,3)=Im(3,3)+Ip(mnum)*(1-dc_norm(3,i)*dc_norm(3,i))* &
3840 vbld(i+1)*vbld(i+1)*0.25d0
3846 ! if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
3847 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3848 .and.(mnum.ne.5)) then
3849 iti=iabs(itype(i,mnum))
3851 Im(1,1)=Im(1,1)+Isc(iti,mnum)*(1-dc_norm(1,inres)* &
3852 dc_norm(1,inres))*vbld(inres)*vbld(inres)
3853 Im(1,2)=Im(1,2)-Isc(iti,mnum)*(dc_norm(1,inres)* &
3854 dc_norm(2,inres))*vbld(inres)*vbld(inres)
3855 Im(1,3)=Im(1,3)-Isc(iti,mnum)*(dc_norm(1,inres)* &
3856 dc_norm(3,inres))*vbld(inres)*vbld(inres)
3857 Im(2,3)=Im(2,3)-Isc(iti,mnum)*(dc_norm(2,inres)* &
3858 dc_norm(3,inres))*vbld(inres)*vbld(inres)
3859 Im(2,2)=Im(2,2)+Isc(iti,mnum)*(1-dc_norm(2,inres)* &
3860 dc_norm(2,inres))*vbld(inres)*vbld(inres)
3861 Im(3,3)=Im(3,3)+Isc(iti,mnum)*(1-dc_norm(3,inres)* &
3862 dc_norm(3,inres))*vbld(inres)*vbld(inres)
3867 ! write(iout,*) "The angular momentum before adjustment:"
3868 ! write(iout,*) (L(j),j=1,3)
3874 ! Copying the Im matrix for the djacob subroutine
3882 ! Finding the eigenvectors and eignvalues of the inertia tensor
3883 call djacob(3,3,10000,1.0d-10,Imcp,eigvec,eigval)
3884 ! write (iout,*) "Eigenvalues & Eigenvectors"
3885 ! write (iout,'(5x,3f10.5)') (eigval(i),i=1,3)
3888 ! write (iout,'(i5,3f10.5)') i,(eigvec(i,j),j=1,3)
3890 ! Constructing the diagonalized matrix
3892 if (dabs(eigval(i)).gt.1.0d-15) then
3893 Id(i,i)=1.0d0/eigval(i)
3900 Imcp(i,j)=eigvec(j,i)
3906 pr1(i,j)=pr1(i,j)+Id(i,k)*Imcp(k,j)
3913 pr2(i,j)=pr2(i,j)+eigvec(i,k)*pr1(k,j)
3917 ! Calculating the total rotational velocity of the molecule
3920 vrot(i)=vrot(i)+pr2(i,j)*L(j)
3923 ! Resetting the velocities
3925 call vecpr(vrot(1),dc(1,i),vp)
3927 d_t(j,i)=d_t(j,i)-vp(j)
3932 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3933 .and.(mnum.ne.5)) then
3934 ! if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
3935 ! if(itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3937 call vecpr(vrot(1),dc(1,inres),vp)
3939 d_t(j,inres)=d_t(j,inres)-vp(j)
3944 ! write(iout,*) "The angular momentum after adjustment:"
3945 ! write(iout,*) (L(j),j=1,3)
3948 end subroutine inertia_tensor
3949 !-----------------------------------------------------------------------------
3950 subroutine angmom(cm,L)
3953 ! implicit real*8 (a-h,o-z)
3954 ! include 'DIMENSIONS'
3955 ! include 'COMMON.CONTROL'
3956 ! include 'COMMON.VAR'
3957 ! include 'COMMON.MD'
3958 ! include 'COMMON.CHAIN'
3959 ! include 'COMMON.DERIV'
3960 ! include 'COMMON.GEO'
3961 ! include 'COMMON.LOCAL'
3962 ! include 'COMMON.INTERACT'
3963 ! include 'COMMON.IOUNITS'
3964 ! include 'COMMON.NAMES'
3965 real(kind=8) :: mscab
3966 real(kind=8),dimension(3) :: L,cm,pr,vp,vrot,incr,v,pp
3967 integer :: iti,inres,i,j,mnum
3968 ! Calculate the angular momentum
3977 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
3979 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
3982 v(j)=incr(j)+0.5d0*d_t(j,i)
3985 incr(j)=incr(j)+d_t(j,i)
3987 call vecpr(pr(1),v(1),vp)
3989 L(j)=L(j)+mp(mnum)*vp(j)
3993 pp(j)=0.5d0*d_t(j,i)
3995 call vecpr(pr(1),pp(1),vp)
3997 L(j)=L(j)+Ip(mnum)*vp(j)
4005 iti=iabs(itype(i,mnum))
4013 pr(j)=c(j,inres)-cm(j)
4015 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4016 .and.(mnum.ne.5)) then
4018 v(j)=incr(j)+d_t(j,inres)
4025 call vecpr(pr(1),v(1),vp)
4026 ! write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),
4027 ! & " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
4029 L(j)=L(j)+mscab*vp(j)
4031 ! write (iout,*) "L",(l(j),j=1,3)
4032 if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4033 .and.(mnum.ne.5)) then
4035 v(j)=incr(j)+d_t(j,inres)
4037 call vecpr(dc(1,inres),d_t(1,inres),vp)
4039 L(j)=L(j)+Isc(iti,mnum)*vp(j)
4043 incr(j)=incr(j)+d_t(j,i)
4047 end subroutine angmom
4048 !-----------------------------------------------------------------------------
4049 subroutine vcm_vel(vcm)
4052 ! implicit real*8 (a-h,o-z)
4053 ! include 'DIMENSIONS'
4054 ! include 'COMMON.VAR'
4055 ! include 'COMMON.MD'
4056 ! include 'COMMON.CHAIN'
4057 ! include 'COMMON.DERIV'
4058 ! include 'COMMON.GEO'
4059 ! include 'COMMON.LOCAL'
4060 ! include 'COMMON.INTERACT'
4061 ! include 'COMMON.IOUNITS'
4062 real(kind=8),dimension(3) :: vcm,vv
4063 real(kind=8) :: summas,amas
4073 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
4075 summas=summas+mp(mnum)
4077 vcm(j)=vcm(j)+mp(mnum)*(vv(j)+0.5d0*d_t(j,i))
4081 amas=msc(iabs(itype(i,mnum)),mnum)
4086 if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4087 .and.(mnum.ne.5)) then
4089 vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres))
4093 vcm(j)=vcm(j)+amas*vv(j)
4097 vv(j)=vv(j)+d_t(j,i)
4100 ! write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas
4102 vcm(j)=vcm(j)/summas
4105 end subroutine vcm_vel
4106 !-----------------------------------------------------------------------------
4108 !-----------------------------------------------------------------------------
4110 ! RATTLE algorithm for velocity Verlet - step 1, UNRES
4114 ! implicit real*8 (a-h,o-z)
4115 ! include 'DIMENSIONS'
4117 ! include 'COMMON.CONTROL'
4118 ! include 'COMMON.VAR'
4119 ! include 'COMMON.MD'
4121 ! include 'COMMON.LANGEVIN'
4123 ! include 'COMMON.LANGEVIN.lang0'
4125 ! include 'COMMON.CHAIN'
4126 ! include 'COMMON.DERIV'
4127 ! include 'COMMON.GEO'
4128 ! include 'COMMON.LOCAL'
4129 ! include 'COMMON.INTERACT'
4130 ! include 'COMMON.IOUNITS'
4131 ! include 'COMMON.NAMES'
4132 ! include 'COMMON.TIME1'
4133 !el real(kind=8) :: gginv(2*nres,2*nres),&
4134 !el gdc(3,2*nres,2*nres)
4135 real(kind=8) :: dC_uncor(3,2*nres) !,&
4136 !el real(kind=8) :: Cmat(2*nres,2*nres)
4137 real(kind=8) :: x(2*nres),xcorr(3,2*nres) !maxres2=2*maxres
4138 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
4139 !el common /przechowalnia/ nbond
4140 integer :: max_rattle = 5
4141 logical :: lprn = .false., lprn1 = .false., not_done
4142 real(kind=8) :: tol_rattle = 1.0d-5
4144 integer :: ii,i,j,jj,l,ind,ind1,nres2
4147 !el /common/ przechowalnia
4149 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4150 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4151 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4153 if (lprn) write (iout,*) "RATTLE1"
4157 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4158 .and.(mnum.ne.5)) nbond=nbond+1
4160 ! Make a folded form of the Ginv-matrix
4173 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4178 if (itype(k,1).ne.10 .and. itype(k,mnum).ne.ntyp1_molec(mnum)&
4179 .and.(mnum.ne.5)) then
4183 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4191 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4202 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4206 if (itype(k,1).ne.10) then
4210 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4218 write (iout,*) "Matrix GGinv"
4219 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
4225 if (iter.gt.max_rattle) then
4226 write (iout,*) "Error - too many iterations in RATTLE."
4229 ! Calculate the matrix C = GG**(-1) dC_old o dC
4234 dC_uncor(j,ind1)=dC(j,i)
4238 if (itype(i,1).ne.10) then
4241 dC_uncor(j,ind1)=dC(j,i+nres)
4250 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
4254 if (itype(k,1).ne.10) then
4257 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
4262 ! Calculate deviations from standard virtual-bond lengths
4266 x(ind)=vbld(i+1)**2-vbl**2
4269 if (itype(i,1).ne.10) then
4271 x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4275 write (iout,*) "Coordinates and violations"
4277 write(iout,'(i5,3f10.5,5x,e15.5)') &
4278 i,(dC_uncor(j,i),j=1,3),x(i)
4280 write (iout,*) "Velocities and violations"
4284 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4285 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4289 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4290 .and.(mnum.ne.5)) then
4293 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4294 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4295 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4298 ! write (iout,*) "gdc"
4300 ! write (iout,*) "i",i
4302 ! write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4308 if (dabs(x(i)).gt.xmax) then
4312 if (xmax.lt.tol_rattle) then
4316 ! Calculate the matrix of the system of equations
4321 Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
4326 write (iout,*) "Matrix Cmat"
4327 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4329 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4330 ! Add constraint term to positions
4337 xx = xx+x(ii)*gdc(j,ind,ii)
4341 d_t_new(j,i)=d_t_new(j,i)-xx/d_time
4345 if (itype(i,1).ne.10) then
4350 xx = xx+x(ii)*gdc(j,ind,ii)
4353 dC(j,i+nres)=dC(j,i+nres)-xx
4354 d_t_new(j,i+nres)=d_t_new(j,i+nres)-xx/d_time
4358 ! Rebuild the chain using the new coordinates
4359 call chainbuild_cart
4361 write (iout,*) "New coordinates, Lagrange multipliers,",&
4362 " and differences between actual and standard bond lengths"
4366 xx=vbld(i+1)**2-vbl**2
4367 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4368 i,(dC(j,i),j=1,3),x(ind),xx
4372 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4375 xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4376 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4377 i,(dC(j,i+nres),j=1,3),x(ind),xx
4380 write (iout,*) "Velocities and violations"
4384 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4385 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4388 if (itype(i,1).ne.10) then
4390 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4391 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4392 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4399 10 write (iout,*) "Error - singularity in solving the system",&
4400 " of equations for Lagrange multipliers."
4404 "RATTLE inactive; use -DRATTLE switch at compile time."
4407 end subroutine rattle1
4408 !-----------------------------------------------------------------------------
4410 ! RATTLE algorithm for velocity Verlet - step 2, UNRES
4414 ! implicit real*8 (a-h,o-z)
4415 ! include 'DIMENSIONS'
4417 ! include 'COMMON.CONTROL'
4418 ! include 'COMMON.VAR'
4419 ! include 'COMMON.MD'
4421 ! include 'COMMON.LANGEVIN'
4423 ! include 'COMMON.LANGEVIN.lang0'
4425 ! include 'COMMON.CHAIN'
4426 ! include 'COMMON.DERIV'
4427 ! include 'COMMON.GEO'
4428 ! include 'COMMON.LOCAL'
4429 ! include 'COMMON.INTERACT'
4430 ! include 'COMMON.IOUNITS'
4431 ! include 'COMMON.NAMES'
4432 ! include 'COMMON.TIME1'
4433 !el real(kind=8) :: gginv(2*nres,2*nres),&
4434 !el gdc(3,2*nres,2*nres)
4435 real(kind=8) :: dC_uncor(3,2*nres) !,&
4436 !el Cmat(2*nres,2*nres)
4437 real(kind=8) :: x(2*nres) !maxres2=2*maxres
4438 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
4439 !el common /przechowalnia/ nbond
4440 integer :: max_rattle = 5
4441 logical :: lprn = .false., lprn1 = .false., not_done
4442 real(kind=8) :: tol_rattle = 1.0d-5
4446 !el /common/ przechowalnia
4447 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4448 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4449 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4451 if (lprn) write (iout,*) "RATTLE2"
4452 if (lprn) write (iout,*) "Velocity correction"
4453 ! Calculate the matrix G dC
4459 gdc(j,i,ind)=GGinv(i,ind)*dC(j,k)
4464 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4465 .and.(mnum.ne.5)) then
4468 gdc(j,i,ind)=GGinv(i,ind)*dC(j,k+nres)
4474 ! write (iout,*) "gdc"
4476 ! write (iout,*) "i",i
4478 ! write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4482 ! Calculate the matrix of the system of equations
4489 Cmat(ind,j)=Cmat(ind,j)+dC(k,i)*gdc(k,ind,j)
4495 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4496 .and.(mnum.ne.5)) then
4501 Cmat(ind,j)=Cmat(ind,j)+dC(k,i+nres)*gdc(k,ind,j)
4506 ! Calculate the scalar product dC o d_t_new
4510 x(ind)=scalar(d_t(1,i),dC(1,i))
4514 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4515 .and.(mnum.ne.5)) then
4517 x(ind)=scalar(d_t(1,i+nres),dC(1,i+nres))
4521 write (iout,*) "Velocities and violations"
4525 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4526 i,ind,(d_t(j,i),j=1,3),x(ind)
4530 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4531 .and.(mnum.ne.5)) then
4533 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4534 i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind)
4540 if (dabs(x(i)).gt.xmax) then
4544 if (xmax.lt.tol_rattle) then
4549 write (iout,*) "Matrix Cmat"
4550 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4552 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4553 ! Add constraint term to velocities
4560 xx = xx+x(ii)*gdc(j,ind,ii)
4562 d_t(j,i)=d_t(j,i)-xx
4567 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4568 .and.(mnum.ne.5)) then
4573 xx = xx+x(ii)*gdc(j,ind,ii)
4575 d_t(j,i+nres)=d_t(j,i+nres)-xx
4581 "New velocities, Lagrange multipliers violations"
4585 if (lprn) write (iout,'(2i5,3f10.5,5x,2e15.5)') &
4586 i,ind,(d_t(j,i),j=1,3),x(ind),scalar(d_t(1,i),dC(1,i))
4590 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4593 write (iout,'(2i5,3f10.5,5x,2e15.5)') &
4594 i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind),&
4595 scalar(d_t(1,i+nres),dC(1,i+nres))
4601 10 write (iout,*) "Error - singularity in solving the system",&
4602 " of equations for Lagrange multipliers."
4606 "RATTLE inactive; use -DRATTLE option at compile time."
4609 end subroutine rattle2
4610 !-----------------------------------------------------------------------------
4611 subroutine rattle_brown
4612 ! RATTLE/LINCS algorithm for Brownian dynamics, UNRES
4616 ! implicit real*8 (a-h,o-z)
4617 ! include 'DIMENSIONS'
4619 ! include 'COMMON.CONTROL'
4620 ! include 'COMMON.VAR'
4621 ! include 'COMMON.MD'
4623 ! include 'COMMON.LANGEVIN'
4625 ! include 'COMMON.LANGEVIN.lang0'
4627 ! include 'COMMON.CHAIN'
4628 ! include 'COMMON.DERIV'
4629 ! include 'COMMON.GEO'
4630 ! include 'COMMON.LOCAL'
4631 ! include 'COMMON.INTERACT'
4632 ! include 'COMMON.IOUNITS'
4633 ! include 'COMMON.NAMES'
4634 ! include 'COMMON.TIME1'
4635 !el real(kind=8) :: gginv(2*nres,2*nres),&
4636 !el gdc(3,2*nres,2*nres)
4637 real(kind=8) :: dC_uncor(3,2*nres) !,&
4638 !el real(kind=8) :: Cmat(2*nres,2*nres)
4639 real(kind=8) :: x(2*nres) !maxres2=2*maxres
4640 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
4641 !el common /przechowalnia/ nbond
4642 integer :: max_rattle = 5
4643 logical :: lprn = .false., lprn1 = .false., not_done
4644 real(kind=8) :: tol_rattle = 1.0d-5
4648 !el /common/ przechowalnia
4649 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4650 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4651 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4654 if (lprn) write (iout,*) "RATTLE_BROWN"
4657 if (itype(i,1).ne.10) nbond=nbond+1
4659 ! Make a folded form of the Ginv-matrix
4672 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4676 if (itype(k,1).ne.10) then
4680 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4687 if (itype(i,1).ne.10) then
4697 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4701 if (itype(k,1).ne.10) then
4705 if (j.eq.1 .and. l.eq.1)GGinv(ii,jj)=fricmat(ind,ind1)
4713 write (iout,*) "Matrix GGinv"
4714 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
4720 if (iter.gt.max_rattle) then
4721 write (iout,*) "Error - too many iterations in RATTLE."
4724 ! Calculate the matrix C = GG**(-1) dC_old o dC
4729 dC_uncor(j,ind1)=dC(j,i)
4733 if (itype(i,1).ne.10) then
4736 dC_uncor(j,ind1)=dC(j,i+nres)
4745 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
4749 if (itype(k,1).ne.10) then
4752 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
4757 ! Calculate deviations from standard virtual-bond lengths
4761 x(ind)=vbld(i+1)**2-vbl**2
4764 if (itype(i,1).ne.10) then
4766 x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4770 write (iout,*) "Coordinates and violations"
4772 write(iout,'(i5,3f10.5,5x,e15.5)') &
4773 i,(dC_uncor(j,i),j=1,3),x(i)
4775 write (iout,*) "Velocities and violations"
4779 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4780 i,ind,(d_t(j,i),j=1,3),scalar(d_t(1,i),dC_old(1,i))
4783 if (itype(i,1).ne.10) then
4785 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4786 i+nres,ind,(d_t(j,i+nres),j=1,3),&
4787 scalar(d_t(1,i+nres),dC_old(1,i+nres))
4790 write (iout,*) "gdc"
4792 write (iout,*) "i",i
4794 write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4800 if (dabs(x(i)).gt.xmax) then
4804 if (xmax.lt.tol_rattle) then
4808 ! Calculate the matrix of the system of equations
4813 Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
4818 write (iout,*) "Matrix Cmat"
4819 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4821 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4822 ! Add constraint term to positions
4829 xx = xx+x(ii)*gdc(j,ind,ii)
4832 d_t(j,i)=d_t(j,i)+xx/d_time
4837 if (itype(i,1).ne.10) then
4842 xx = xx+x(ii)*gdc(j,ind,ii)
4845 d_t(j,i+nres)=d_t(j,i+nres)+xx/d_time
4846 dC(j,i+nres)=dC(j,i+nres)+xx
4850 ! Rebuild the chain using the new coordinates
4851 call chainbuild_cart
4853 write (iout,*) "New coordinates, Lagrange multipliers,",&
4854 " and differences between actual and standard bond lengths"
4858 xx=vbld(i+1)**2-vbl**2
4859 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4860 i,(dC(j,i),j=1,3),x(ind),xx
4863 if (itype(i,1).ne.10) then
4865 xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4866 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4867 i,(dC(j,i+nres),j=1,3),x(ind),xx
4870 write (iout,*) "Velocities and violations"
4874 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4875 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4878 if (itype(i,1).ne.10) then
4880 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4881 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4882 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4889 10 write (iout,*) "Error - singularity in solving the system",&
4890 " of equations for Lagrange multipliers."
4894 "RATTLE inactive; use -DRATTLE option at compile time"
4897 end subroutine rattle_brown
4898 !-----------------------------------------------------------------------------
4900 !-----------------------------------------------------------------------------
4901 subroutine friction_force
4906 ! implicit real*8 (a-h,o-z)
4907 ! include 'DIMENSIONS'
4908 ! include 'COMMON.VAR'
4909 ! include 'COMMON.CHAIN'
4910 ! include 'COMMON.DERIV'
4911 ! include 'COMMON.GEO'
4912 ! include 'COMMON.LOCAL'
4913 ! include 'COMMON.INTERACT'
4914 ! include 'COMMON.MD'
4916 ! include 'COMMON.LANGEVIN'
4918 ! include 'COMMON.LANGEVIN.lang0'
4920 ! include 'COMMON.IOUNITS'
4921 !el real(kind=8),dimension(6*nres) :: gamvec !(MAXRES6) maxres6=6*maxres
4922 !el common /syfek/ gamvec
4923 real(kind=8) :: vv(3),vvtot(3,nres),v_work(6*nres) !,&
4924 !el ginvfric(2*nres,2*nres) !maxres2=2*maxres
4925 !el common /przechowalnia/ ginvfric
4927 logical :: lprn = .false., checkmode = .false.
4928 integer :: i,j,ind,k,nres2,nres6,mnum
4932 if(.not.allocated(gamvec)) allocate(gamvec(nres6)) !(MAXRES6)
4933 if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
4941 d_t_work(j)=d_t(j,0)
4946 d_t_work(ind+j)=d_t(j,i)
4952 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
4953 .and.(mnum.ne.5)) then
4955 d_t_work(ind+j)=d_t(j,i+nres)
4961 call fricmat_mult(d_t_work,fric_work)
4963 if (.not.checkmode) return
4966 write (iout,*) "d_t_work and fric_work"
4968 write (iout,'(i3,2e15.5)') i,d_t_work(i),fric_work(i)
4972 friction(j,0)=fric_work(j)
4977 friction(j,i)=fric_work(ind+j)
4983 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
4984 .and.(mnum.ne.5)) then
4985 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
4987 friction(j,i+nres)=fric_work(ind+j)
4993 write(iout,*) "Friction backbone"
4995 write(iout,'(i5,3e15.5,5x,3e15.5)') &
4996 i,(friction(j,i),j=1,3),(d_t(j,i),j=1,3)
4998 write(iout,*) "Friction side chain"
5000 write(iout,'(i5,3e15.5,5x,3e15.5)') &
5001 i,(friction(j,i+nres),j=1,3),(d_t(j,i+nres),j=1,3)
5010 vvtot(j,i)=vv(j)+0.5d0*d_t(j,i)
5011 vvtot(j,i+nres)=vv(j)+d_t(j,i+nres)
5012 vv(j)=vv(j)+d_t(j,i)
5015 write (iout,*) "vvtot backbone and sidechain"
5017 write (iout,'(i5,3e15.5,5x,3e15.5)') i,(vvtot(j,i),j=1,3),&
5018 (vvtot(j,i+nres),j=1,3)
5023 v_work(ind+j)=vvtot(j,i)
5029 v_work(ind+j)=vvtot(j,i+nres)
5033 write (iout,*) "v_work gamvec and site-based friction forces"
5035 write (iout,'(i5,3e15.5)') i,v_work(i),gamvec(i),&
5039 ! fric_work1(i)=0.0d0
5041 ! fric_work1(i)=fric_work1(i)-A(j,i)*gamvec(j)*v_work(j)
5044 ! write (iout,*) "fric_work and fric_work1"
5046 ! write (iout,'(i5,2e15.5)') i,fric_work(i),fric_work1(i)
5052 ginvfric(i,j)=ginvfric(i,j)+ginv(i,k)*fricmat(k,j)
5056 write (iout,*) "ginvfric"
5058 write (iout,'(i5,100f8.3)') i,(ginvfric(i,j),j=1,dimen)
5060 write (iout,*) "symmetry check"
5063 write (iout,*) i,j,ginvfric(i,j)-ginvfric(j,i)
5068 end subroutine friction_force
5069 !-----------------------------------------------------------------------------
5070 subroutine setup_fricmat
5074 use control_data, only:time_Bcast
5075 use control, only:tcpu
5077 ! implicit real*8 (a-h,o-z)
5081 real(kind=8) :: time00
5083 ! include 'DIMENSIONS'
5084 ! include 'COMMON.VAR'
5085 ! include 'COMMON.CHAIN'
5086 ! include 'COMMON.DERIV'
5087 ! include 'COMMON.GEO'
5088 ! include 'COMMON.LOCAL'
5089 ! include 'COMMON.INTERACT'
5090 ! include 'COMMON.MD'
5091 ! include 'COMMON.SETUP'
5092 ! include 'COMMON.TIME1'
5093 ! integer licznik /0/
5096 ! include 'COMMON.LANGEVIN'
5098 ! include 'COMMON.LANGEVIN.lang0'
5100 ! include 'COMMON.IOUNITS'
5102 integer :: i,j,ind,ind1,m
5103 logical :: lprn = .false.
5104 real(kind=8) :: dtdi !el ,gamvec(2*nres)
5105 !el real(kind=8),dimension(2*nres,2*nres) :: ginvfric,fcopy
5106 ! real(kind=8),allocatable,dimension(:,:) :: fcopy
5107 !el real(kind=8),dimension(2*nres*(2*nres+1)/2) :: Ghalf !(mmaxres2) (mmaxres2=(maxres2*(maxres2+1)/2))
5108 !el common /syfek/ gamvec
5109 real(kind=8) :: work(8*2*nres)
5110 integer :: iwork(2*nres)
5111 !el common /przechowalnia/ ginvfric,Ghalf,fcopy
5112 integer :: ii,iti,k,l,nzero,nres2,nres6,ierr,mnum
5116 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
5117 if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
5118 if (fg_rank.ne.king) goto 10
5123 if(.not.allocated(gamvec)) allocate(gamvec(nres2)) !(MAXRES2)
5124 if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
5125 if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
5126 !el allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
5127 if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !maxres2=2*maxres
5129 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
5130 ! Zeroing out fricmat
5136 ! Load the friction coefficients corresponding to peptide groups
5141 gamvec(ind1)=gamp(mnum)
5143 ! Load the friction coefficients corresponding to side chains
5147 gamsc(ntyp1_molec(j),j)=1.0d0
5154 gamvec(ii)=gamsc(iabs(iti),mnum)
5156 if (surfarea) call sdarea(gamvec)
5158 ! write (iout,*) "Matrix A and vector gamma"
5160 ! write (iout,'(i2,$)') i
5162 ! write (iout,'(f4.1,$)') A(i,j)
5164 ! write (iout,'(f8.3)') gamvec(i)
5168 write (iout,*) "Vector gamvec"
5170 write (iout,'(i5,f10.5)') i, gamvec(i)
5174 ! The friction matrix
5179 dtdi=dtdi+A(j,k)*A(j,i)*gamvec(j)
5186 write (iout,'(//a)') "Matrix fricmat"
5187 call matout2(dimen,dimen,nres2,nres2,fricmat)
5189 if (lang.eq.2 .or. lang.eq.3) then
5190 ! Mass-scale the friction matrix if non-direct integration will be performed
5196 Ginvfric(i,j)=Ginvfric(i,j)+ &
5197 Gsqrm(i,k)*Gsqrm(l,j)*fricmat(k,l)
5202 ! Diagonalize the friction matrix
5207 Ghalf(ind)=Ginvfric(i,j)
5210 call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
5213 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
5214 " mass-scaled friction matrix"
5215 call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
5217 ! Precompute matrices for tinker stochastic integrator
5224 mt1(i,j)=mt1(i,j)+fricvec(k,i)*gsqrm(k,j)
5225 mt2(i,j)=mt2(i,j)+fricvec(k,i)*gsqrp(k,j)
5231 else if (lang.eq.4) then
5232 ! Diagonalize the friction matrix
5237 Ghalf(ind)=fricmat(i,j)
5240 call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
5243 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
5245 call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
5247 ! Determine the number of zero eigenvalues of the friction matrix
5248 nzero=max0(dimen-dimen1,0)
5249 ! do while (fricgam(nzero+1).le.1.0d-5 .and. nzero.lt.dimen)
5252 write (iout,*) "Number of zero eigenvalues:",nzero
5257 fricmat(i,j)=fricmat(i,j) &
5258 +fricvec(i,k)*fricvec(j,k)/fricgam(k)
5263 write (iout,'(//a)') "Generalized inverse of fricmat"
5264 call matout(dimen,dimen,nres6,nres6,fricmat)
5269 if (nfgtasks.gt.1) then
5270 if (fg_rank.eq.0) then
5271 ! The matching BROADCAST for fg processors is called in ERGASTULUM
5277 call MPI_Bcast(10,1,MPI_INTEGER,king,FG_COMM,IERROR)
5279 time_Bcast=time_Bcast+MPI_Wtime()-time00
5281 time_Bcast=time_Bcast+tcpu()-time00
5283 ! print *,"Processor",myrank,
5284 ! & " BROADCAST iorder in SETUP_FRICMAT"
5287 write (iout,*) "setup_fricmat licznik"!,licznik !sp
5293 ! Scatter the friction matrix
5294 call MPI_Scatterv(fricmat(1,1),nginv_counts(0),&
5295 nginv_start(0),MPI_DOUBLE_PRECISION,fcopy(1,1),&
5296 myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
5299 time_scatter=time_scatter+MPI_Wtime()-time00
5300 time_scatter_fmat=time_scatter_fmat+MPI_Wtime()-time00
5302 time_scatter=time_scatter+tcpu()-time00
5303 time_scatter_fmat=time_scatter_fmat+tcpu()-time00
5307 do j=1,2*my_ng_count
5308 fricmat(j,i)=fcopy(i,j)
5311 ! write (iout,*) "My chunk of fricmat"
5312 ! call MATOUT2(my_ng_count,dimen,maxres2,maxres2,fcopy)
5316 end subroutine setup_fricmat
5317 !-----------------------------------------------------------------------------
5318 subroutine stochastic_force(stochforcvec)
5321 use random, only:anorm_distr
5322 ! implicit real*8 (a-h,o-z)
5323 ! include 'DIMENSIONS'
5324 use control, only: tcpu
5329 ! include 'COMMON.VAR'
5330 ! include 'COMMON.CHAIN'
5331 ! include 'COMMON.DERIV'
5332 ! include 'COMMON.GEO'
5333 ! include 'COMMON.LOCAL'
5334 ! include 'COMMON.INTERACT'
5335 ! include 'COMMON.MD'
5336 ! include 'COMMON.TIME1'
5338 ! include 'COMMON.LANGEVIN'
5340 ! include 'COMMON.LANGEVIN.lang0'
5342 ! include 'COMMON.IOUNITS'
5344 real(kind=8) :: x,sig,lowb,highb
5345 real(kind=8) :: ff(3),force(3,0:2*nres),zeta2,lowb2
5346 real(kind=8) :: highb2,sig2,forcvec(6*nres),stochforcvec(6*nres)
5347 real(kind=8) :: time00
5348 logical :: lprn = .false.
5349 integer :: i,j,ind,mnum
5353 stochforc(j,i)=0.0d0
5363 ! Compute the stochastic forces acting on bodies. Store in force.
5369 force(j,i)=anorm_distr(x,sig,lowb,highb)
5377 force(j,i+nres)=anorm_distr(x,sig2,lowb2,highb2)
5381 time_fsample=time_fsample+MPI_Wtime()-time00
5383 time_fsample=time_fsample+tcpu()-time00
5385 ! Compute the stochastic forces acting on virtual-bond vectors.
5391 stochforc(j,i)=ff(j)+0.5d0*force(j,i)
5394 ff(j)=ff(j)+force(j,i)
5396 ! if (itype(i+1,1).ne.ntyp1) then
5398 if (itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
5400 stochforc(j,i)=stochforc(j,i)+force(j,i+nres+1)
5401 ff(j)=ff(j)+force(j,i+nres+1)
5406 stochforc(j,0)=ff(j)+force(j,nnt+nres)
5410 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5411 .and.(mnum.ne.5)) then
5412 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
5414 stochforc(j,i+nres)=force(j,i+nres)
5420 stochforcvec(j)=stochforc(j,0)
5425 stochforcvec(ind+j)=stochforc(j,i)
5431 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5432 .and.(mnum.ne.5)) then
5433 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
5435 stochforcvec(ind+j)=stochforc(j,i+nres)
5441 write (iout,*) "stochforcvec"
5443 write(iout,'(i5,e15.5)') i,stochforcvec(i)
5445 write(iout,*) "Stochastic forces backbone"
5447 write(iout,'(i5,3e15.5)') i,(stochforc(j,i),j=1,3)
5449 write(iout,*) "Stochastic forces side chain"
5451 write(iout,'(i5,3e15.5)') &
5452 i,(stochforc(j,i+nres),j=1,3)
5460 write (iout,*) i,ind
5462 forcvec(ind+j)=force(j,i)
5467 write (iout,*) i,ind
5469 forcvec(j+ind)=force(j,i+nres)
5474 write (iout,*) "forcvec"
5478 write (iout,'(2i3,2f10.5)') i,j,force(j,i),&
5485 write (iout,'(2i3,2f10.5)') i,j,force(j,i+nres),&
5494 end subroutine stochastic_force
5495 !-----------------------------------------------------------------------------
5496 subroutine sdarea(gamvec)
5498 ! Scale the friction coefficients according to solvent accessible surface areas
5499 ! Code adapted from TINKER
5503 ! implicit real*8 (a-h,o-z)
5504 ! include 'DIMENSIONS'
5505 ! include 'COMMON.CONTROL'
5506 ! include 'COMMON.VAR'
5507 ! include 'COMMON.MD'
5509 ! include 'COMMON.LANGEVIN'
5511 ! include 'COMMON.LANGEVIN.lang0'
5513 ! include 'COMMON.CHAIN'
5514 ! include 'COMMON.DERIV'
5515 ! include 'COMMON.GEO'
5516 ! include 'COMMON.LOCAL'
5517 ! include 'COMMON.INTERACT'
5518 ! include 'COMMON.IOUNITS'
5519 ! include 'COMMON.NAMES'
5520 real(kind=8),dimension(2*nres) :: radius,gamvec !(maxres2)
5521 real(kind=8),parameter :: twosix = 1.122462048309372981d0
5522 logical :: lprn = .false.
5523 real(kind=8) :: probe,area,ratio
5524 integer :: i,j,ind,iti,mnum
5526 ! determine new friction coefficients every few SD steps
5528 ! set the atomic radii to estimates of sigma values
5530 ! print *,"Entered sdarea"
5536 ! Load peptide group radii
5539 radius(i)=pstok(mnum)
5541 ! Load side chain radii
5545 radius(i+nres)=restok(iti,mnum)
5548 ! write (iout,*) "i",i," radius",radius(i)
5551 radius(i) = radius(i) / twosix
5552 if (radius(i) .ne. 0.0d0) radius(i) = radius(i) + probe
5555 ! scale atomic friction coefficients by accessible area
5557 if (lprn) write (iout,*) &
5558 "Original gammas, surface areas, scaling factors, new gammas, ",&
5559 "std's of stochastic forces"
5562 if (radius(i).gt.0.0d0) then
5563 call surfatom (i,area,radius)
5564 ratio = dmax1(area/(4.0d0*pi*radius(i)**2),1.0d-1)
5565 if (lprn) write (iout,'(i5,3f10.5,$)') &
5566 i,gamvec(ind+1),area,ratio
5569 gamvec(ind) = ratio * gamvec(ind)
5572 stdforcp(i)=stdfp(mnum)*dsqrt(gamvec(ind))
5573 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcp(i)
5577 if (radius(i+nres).gt.0.0d0) then
5578 call surfatom (i+nres,area,radius)
5579 ratio = dmax1(area/(4.0d0*pi*radius(i+nres)**2),1.0d-1)
5580 if (lprn) write (iout,'(i5,3f10.5,$)') &
5581 i,gamvec(ind+1),area,ratio
5584 gamvec(ind) = ratio * gamvec(ind)
5587 stdforcsc(i)=stdfsc(itype(i,mnum),mnum)*dsqrt(gamvec(ind))
5588 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcsc(i)
5593 end subroutine sdarea
5594 !-----------------------------------------------------------------------------
5596 !-----------------------------------------------------------------------------
5599 ! ###################################################
5600 ! ## COPYRIGHT (C) 1996 by Jay William Ponder ##
5601 ! ## All Rights Reserved ##
5602 ! ###################################################
5604 ! ################################################################
5606 ! ## subroutine surfatom -- exposed surface area of an atom ##
5608 ! ################################################################
5611 ! "surfatom" performs an analytical computation of the surface
5612 ! area of a specified atom; a simplified version of "surface"
5614 ! literature references:
5616 ! T. J. Richmond, "Solvent Accessible Surface Area and
5617 ! Excluded Volume in Proteins", Journal of Molecular Biology,
5620 ! L. Wesson and D. Eisenberg, "Atomic Solvation Parameters
5621 ! Applied to Molecular Dynamics of Proteins in Solution",
5622 ! Protein Science, 1, 227-235 (1992)
5624 ! variables and parameters:
5626 ! ir number of atom for which area is desired
5627 ! area accessible surface area of the atom
5628 ! radius radii of each of the individual atoms
5631 subroutine surfatom(ir,area,radius)
5633 ! implicit real*8 (a-h,o-z)
5634 ! include 'DIMENSIONS'
5636 ! include 'COMMON.GEO'
5637 ! include 'COMMON.IOUNITS'
5639 integer :: nsup,nstart_sup
5640 ! double precision c,dc,dc_old,d_c_work,xloc,xrot,dc_norm
5641 ! common /chain/ c(3,maxres2+2),dc(3,0:maxres2),dc_old(3,0:maxres2),
5642 ! & xloc(3,maxres),xrot(3,maxres),dc_norm(3,0:maxres2),
5643 ! & dc_work(MAXRES6),nres,nres0
5644 integer,parameter :: maxarc=300
5648 integer :: mi,ni,narc
5649 integer :: key(maxarc)
5650 integer :: intag(maxarc)
5651 integer :: intag1(maxarc)
5652 real(kind=8) :: area,arcsum
5653 real(kind=8) :: arclen,exang
5654 real(kind=8) :: delta,delta2
5655 real(kind=8) :: eps,rmove
5656 real(kind=8) :: xr,yr,zr
5657 real(kind=8) :: rr,rrsq
5658 real(kind=8) :: rplus,rminus
5659 real(kind=8) :: axx,axy,axz
5660 real(kind=8) :: ayx,ayy
5661 real(kind=8) :: azx,azy,azz
5662 real(kind=8) :: uxj,uyj,uzj
5663 real(kind=8) :: tx,ty,tz
5664 real(kind=8) :: txb,tyb,td
5665 real(kind=8) :: tr2,tr,txr,tyr
5666 real(kind=8) :: tk1,tk2
5667 real(kind=8) :: thec,the,t,tb
5668 real(kind=8) :: txk,tyk,tzk
5669 real(kind=8) :: t1,ti,tf,tt
5670 real(kind=8) :: txj,tyj,tzj
5671 real(kind=8) :: ccsq,cc,xysq
5672 real(kind=8) :: bsqk,bk,cosine
5673 real(kind=8) :: dsqj,gi,pix2
5674 real(kind=8) :: therk,dk,gk
5675 real(kind=8) :: risqk,rik
5676 real(kind=8) :: radius(2*nres) !(maxatm) (maxatm=maxres2)
5677 real(kind=8) :: ri(maxarc),risq(maxarc)
5678 real(kind=8) :: ux(maxarc),uy(maxarc),uz(maxarc)
5679 real(kind=8) :: xc(maxarc),yc(maxarc),zc(maxarc)
5680 real(kind=8) :: xc1(maxarc),yc1(maxarc),zc1(maxarc)
5681 real(kind=8) :: dsq(maxarc),bsq(maxarc)
5682 real(kind=8) :: dsq1(maxarc),bsq1(maxarc)
5683 real(kind=8) :: arci(maxarc),arcf(maxarc)
5684 real(kind=8) :: ex(maxarc),lt(maxarc),gr(maxarc)
5685 real(kind=8) :: b(maxarc),b1(maxarc),bg(maxarc)
5686 real(kind=8) :: kent(maxarc),kout(maxarc)
5687 real(kind=8) :: ther(maxarc)
5688 logical :: moved,top
5689 logical :: omit(maxarc)
5692 ! maxatm = 2*nres !maxres2 maxres2=2*maxres
5693 ! maxlight = 8*maxatm
5696 ! maxtors = 4*maxatm
5698 ! zero out the surface area for the sphere of interest
5701 ! write (2,*) "ir",ir," radius",radius(ir)
5702 if (radius(ir) .eq. 0.0d0) return
5704 ! set the overlap significance and connectivity shift
5708 delta2 = delta * delta
5713 ! store coordinates and radius of the sphere of interest
5721 ! initialize values of some counters and summations
5730 ! test each sphere to see if it overlaps the sphere of interest
5733 if (i.eq.ir .or. radius(i).eq.0.0d0) goto 30
5734 rplus = rr + radius(i)
5736 if (abs(tx) .ge. rplus) goto 30
5738 if (abs(ty) .ge. rplus) goto 30
5740 if (abs(tz) .ge. rplus) goto 30
5742 ! check for sphere overlap by testing distance against radii
5744 xysq = tx*tx + ty*ty
5745 if (xysq .lt. delta2) then
5752 if (rplus-cc .le. delta) goto 30
5753 rminus = rr - radius(i)
5755 ! check to see if sphere of interest is completely buried
5757 if (cc-abs(rminus) .le. delta) then
5758 if (rminus .le. 0.0d0) goto 170
5762 ! check for too many overlaps with sphere of interest
5764 if (io .ge. maxarc) then
5766 20 format (/,' SURFATOM -- Increase the Value of MAXARC')
5770 ! get overlap between current sphere and sphere of interest
5779 gr(io) = (ccsq+rplus*rminus) / (2.0d0*rr*b1(io))
5785 ! case where no other spheres overlap the sphere of interest
5788 area = 4.0d0 * pi * rrsq
5792 ! case where only one sphere overlaps the sphere of interest
5795 area = pix2 * (1.0d0 + gr(1))
5796 area = mod(area,4.0d0*pi) * rrsq
5800 ! case where many spheres intersect the sphere of interest;
5801 ! sort the intersecting spheres by their degree of overlap
5803 call sort2 (io,gr,key)
5806 intag(i) = intag1(k)
5815 ! get radius of each overlap circle on surface of the sphere
5820 risq(i) = rrsq - gi*gi
5821 ri(i) = sqrt(risq(i))
5822 ther(i) = 0.5d0*pi - asin(min(1.0d0,max(-1.0d0,gr(i))))
5825 ! find boundary of inaccessible area on sphere of interest
5828 if (.not. omit(k)) then
5835 ! check to see if J circle is intersecting K circle;
5836 ! get distance between circle centers and sum of radii
5839 if (omit(j)) goto 60
5840 cc = (txk*xc(j)+tyk*yc(j)+tzk*zc(j))/(bk*b(j))
5841 cc = acos(min(1.0d0,max(-1.0d0,cc)))
5842 td = therk + ther(j)
5844 ! check to see if circles enclose separate regions
5846 if (cc .ge. td) goto 60
5848 ! check for circle J completely inside circle K
5850 if (cc+ther(j) .lt. therk) goto 40
5852 ! check for circles that are essentially parallel
5854 if (cc .gt. delta) goto 50
5859 ! check to see if sphere of interest is completely buried
5862 if (pix2-cc .le. td) goto 170
5868 ! find T value of circle intersections
5871 if (omit(k)) goto 110
5886 ! rotation matrix elements
5898 if (.not. omit(j)) then
5903 ! rotate spheres so K vector colinear with z-axis
5905 uxj = txj*axx + tyj*axy - tzj*axz
5906 uyj = tyj*ayy - txj*ayx
5907 uzj = txj*azx + tyj*azy + tzj*azz
5908 cosine = min(1.0d0,max(-1.0d0,uzj/b(j)))
5909 if (acos(cosine) .lt. therk+ther(j)) then
5910 dsqj = uxj*uxj + uyj*uyj
5915 tr2 = risqk*dsqj - tb*tb
5921 ! get T values of intersection for K circle
5924 tb = min(1.0d0,max(-1.0d0,tb))
5926 if (tyb-txr .lt. 0.0d0) tk1 = pix2 - tk1
5928 tb = min(1.0d0,max(-1.0d0,tb))
5930 if (tyb+txr .lt. 0.0d0) tk2 = pix2 - tk2
5931 thec = (rrsq*uzj-gk*bg(j)) / (rik*ri(j)*b(j))
5932 if (abs(thec) .lt. 1.0d0) then
5934 else if (thec .ge. 1.0d0) then
5936 else if (thec .le. -1.0d0) then
5940 ! see if "tk1" is entry or exit point; check t=0 point;
5941 ! "ti" is exit point, "tf" is entry point
5943 cosine = min(1.0d0,max(-1.0d0, &
5944 (uzj*gk-uxj*rik)/(b(j)*rr)))
5945 if ((acos(cosine)-ther(j))*(tk2-tk1) .le. 0.0d0) then
5953 if (narc .ge. maxarc) then
5955 70 format (/,' SURFATOM -- Increase the Value',&
5959 if (tf .le. ti) then
5980 ! special case; K circle without intersections
5982 if (narc .le. 0) goto 90
5984 ! general case; sum up arclength and set connectivity code
5986 call sort2 (narc,arci,key)
5991 if (narc .gt. 1) then
5994 if (t .lt. arci(j)) then
5995 arcsum = arcsum + arci(j) - t
5996 exang = exang + ex(ni)
5998 if (jb .ge. maxarc) then
6000 80 format (/,' SURFATOM -- Increase the Value',&
6005 kent(jb) = maxarc*i + k
6007 kout(jb) = maxarc*k + i
6016 arcsum = arcsum + pix2 - t
6018 exang = exang + ex(ni)
6021 kent(jb) = maxarc*i + k
6023 kout(jb) = maxarc*k + i
6030 arclen = arclen + gr(k)*arcsum
6033 if (arclen .eq. 0.0d0) goto 170
6034 if (jb .eq. 0) goto 150
6036 ! find number of independent boundaries and check connectivity
6040 if (kout(k) .ne. 0) then
6047 if (m .eq. kent(ii)) then
6050 if (j .eq. jb) goto 150
6062 ! attempt to fix connectivity error by moving atom slightly
6066 140 format (/,' SURFATOM -- Connectivity Error at Atom',i6)
6075 ! compute the exposed surface area for the sphere of interest
6078 area = ib*pix2 + exang + arclen
6079 area = mod(area,4.0d0*pi) * rrsq
6081 ! attempt to fix negative area by moving atom slightly
6083 if (area .lt. 0.0d0) then
6086 160 format (/,' SURFATOM -- Negative Area at Atom',i6)
6097 end subroutine surfatom
6098 !----------------------------------------------------------------
6099 !----------------------------------------------------------------
6100 subroutine alloc_MD_arrays
6101 !EL Allocation of arrays used by MD module
6103 integer :: nres2,nres6
6106 !----------------------
6110 allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
6111 allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
6112 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
6113 allocate(fricvec(nres2,nres2))
6114 allocate(pfric_mat(nres2,nres2),vfric_mat(nres2,nres2))
6115 allocate(afric_mat(nres2,nres2),prand_mat(nres2,nres2))
6116 allocate(vrand_mat1(nres2,nres2),vrand_mat2(nres2,nres2)) !(MAXRES2,MAXRES2)
6117 allocate(pfric0_mat(nres2,nres2,0:maxflag_stoch))
6118 allocate(afric0_mat(nres2,nres2,0:maxflag_stoch))
6119 allocate(vfric0_mat(nres2,nres2,0:maxflag_stoch))
6120 allocate(prand0_mat(nres2,nres2,0:maxflag_stoch))
6121 allocate(vrand0_mat1(nres2,nres2,0:maxflag_stoch))
6122 allocate(vrand0_mat2(nres2,nres2,0:maxflag_stoch)) !(MAXRES2,MAXRES2,0:maxflag_stoch)
6123 allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
6125 allocate(mt1(nres2,nres2),mt2(nres2,nres2),mt3(nres2,nres2)) !(maxres2,maxres2)
6126 !----------------------
6128 ! commom.langevin.lang0
6130 allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
6131 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
6132 allocate(fricvec(nres2,nres2)) !(MAXRES2,MAXRES2)
6133 allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
6134 allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
6137 if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2))
6138 !----------------------
6139 ! commom.hairpin in CSA module
6140 !----------------------
6141 ! common.mce in MCM_MD module
6142 !----------------------
6144 ! common /mdgrad/ in module.energy
6145 ! common /back_constr/ in module.energy
6146 ! common /qmeas/ in module.energy
6149 allocate(potEcomp(0:n_ene+4)) !(0:n_ene+4)
6151 allocate(d_t(3,0:nres2),d_a(3,0:nres2),d_t_old(3,0:nres2)) !(3,0:MAXRES2)
6152 allocate(d_a_work(nres6)) !(6*MAXRES)
6154 allocate(DM(nres2),DU1(nres2),DU2(nres2))
6155 allocate(DMorig(nres2),DU1orig(nres2),DU2orig(nres2))
6157 allocate(Gmat(nres2,nres2),A(nres2,nres2))
6158 if(.not.allocated(Ginv)) allocate(Ginv(nres2,nres2)) !in control: ergastulum
6159 allocate(Gsqrp(nres2,nres2),Gsqrm(nres2,nres2),Gvec(nres2,nres2)) !(maxres2,maxres2)
6161 allocate(Geigen(nres2)) !(maxres2)
6162 if(.not.allocated(vtot)) allocate(vtot(nres2)) !(maxres2)
6163 ! common /inertia/ in io_conf: parmread
6164 ! real(kind=8),dimension(:),allocatable :: ISC,msc !(ntyp+1)
6165 ! common /langevin/in io read_MDpar
6166 ! real(kind=8),dimension(:),allocatable :: gamsc !(ntyp1)
6167 ! real(kind=8),dimension(:),allocatable :: stdfsc !(ntyp)
6168 ! in io_conf: parmread
6169 ! real(kind=8),dimension(:),allocatable :: restok !(ntyp+1)
6170 ! common /mdpmpi/ in control: ergastulum
6171 if(.not.allocated(ng_start)) allocate(ng_start(0:nfgtasks-1))
6172 if(.not.allocated(ng_counts)) allocate(ng_counts(0:nfgtasks-1))
6173 if(.not.allocated(nginv_counts)) allocate(nginv_counts(0:nfgtasks-1)) !(0:MaxProcs-1)
6174 if(.not.allocated(nginv_start)) allocate(nginv_start(0:nfgtasks)) !(0:MaxProcs)
6175 !----------------------
6176 ! common.muca in read_muca
6177 ! common /double_muca/
6178 ! real(kind=8) :: elow,ehigh,factor,hbin,factor_min
6179 ! real(kind=8),dimension(:),allocatable :: emuca,nemuca,&
6180 ! nemuca2,hist !(4*maxres)
6181 ! common /integer_muca/
6182 ! integer :: nmuca,imtime,muca_smooth
6184 ! real(kind=8),dimension(:),allocatable :: elowi,ehighi !(maxprocs)
6185 !----------------------
6187 ! common /mdgrad/ in module.energy
6188 ! common /back_constr/ in module.energy
6189 ! common /qmeas/ in module.energy
6193 allocate(d_t_work(nres6),d_t_work_new(nres6),d_af_work(nres6))
6194 allocate(d_as_work(nres6),kinetic_force(nres6)) !(MAXRES6)
6195 allocate(d_t_new(3,0:nres2),d_a_old(3,0:nres2),d_a_short(3,0:nres2)) !,d_a !(3,0:MAXRES2)
6196 allocate(stdforcp(nres),stdforcsc(nres)) !(MAXRES)
6197 !----------------------
6199 allocate(D_ban(nres6)) !(MAXRES6) maxres6=6*maxres
6200 ! common /stochcalc/ stochforcvec
6201 allocate(stochforcvec(nres6)) !(MAXRES6) maxres6=6*maxres
6202 !----------------------
6204 end subroutine alloc_MD_arrays
6205 !-----------------------------------------------------------------------------
6206 !-----------------------------------------------------------------------------