2 !-----------------------------------------------------------------------------
15 !-----------------------------------------------------------------------------
17 ! common /mdgrad/ in module.energy
18 ! common /back_constr/ in module.energy
19 ! common /qmeas/ in module.energy
23 real(kind=8),dimension(:),allocatable :: d_t_work,&
24 d_t_work_new,d_af_work,d_as_work,kinetic_force !(MAXRES6)
25 real(kind=8),dimension(:,:),allocatable :: d_t_new,&
26 d_a_old,d_a_short!,d_a !(3,0:MAXRES2)
27 ! real(kind=8),dimension(:),allocatable :: d_a_work !(6*MAXRES)
28 ! real(kind=8),dimension(:,:),allocatable :: Gmat,Ginv,A,&
29 ! Gsqrp,Gsqrm,Gvec !(maxres2,maxres2)
30 ! real(kind=8),dimension(:),allocatable :: Geigen !(maxres2)
31 ! integer :: dimen,dimen1,dimen3
32 ! integer :: lang,count_reset_moment,count_reset_vel
33 ! logical :: reset_moment,reset_vel,rattle,RESPA
36 ! real(kind=8) :: rwat,etawat,stdfp,cPoise
37 ! real(kind=8),dimension(:),allocatable :: gamsc !(ntyp1)
38 ! real(kind=8),dimension(:),allocatable :: stdfsc !(ntyp)
39 real(kind=8),dimension(:),allocatable :: stdforcp,stdforcsc !(MAXRES)
40 !-----------------------------------------------------------------------------
44 ! ###################################################
45 ! ## COPYRIGHT (C) 1992 by Jay William Ponder ##
46 ! ## All Rights Reserved ##
47 ! ###################################################
49 ! #############################################################
51 ! ## sizes.i -- parameter values to set array dimensions ##
53 ! #############################################################
56 ! "sizes.i" sets values for critical array dimensions used
57 ! throughout the software; these parameters will fix the size
58 ! of the largest systems that can be handled; values too large
59 ! for the computer's memory and/or swap space to accomodate
60 ! will result in poor performance or outright failure
62 ! parameter: maximum allowed number of:
64 ! maxatm atoms in the molecular system
65 ! maxval atoms directly bonded to an atom
66 ! maxgrp ! user-defined groups of atoms
67 ! maxtyp force field atom type definitions
68 ! maxclass force field atom class definitions
69 ! maxkey lines in the keyword file
70 ! maxrot bonds for torsional rotation
71 ! maxvar optimization variables (vector storage)
72 ! maxopt optimization variables (matrix storage)
73 ! maxhess off-diagonal Hessian elements
74 ! maxlight sites for method of lights neighbors
75 ! maxvib vibrational frequencies
76 ! maxgeo distance geometry points
77 ! maxcell unit cells in replicated crystal
78 ! maxring 3-, 4-, or 5-membered rings
79 ! maxfix geometric restraints
80 ! maxbio biopolymer atom definitions
81 ! maxres residues in the macromolecule
82 ! maxamino amino acid residue types
83 ! maxnuc nucleic acid residue types
84 ! maxbnd covalent bonds in molecular system
85 ! maxang bond angles in molecular system
86 ! maxtors torsional angles in molecular system
87 ! maxpi atoms in conjugated pisystem
88 ! maxpib covalent bonds involving pisystem
89 ! maxpit torsional angles involving pisystem
92 !el integer maxatm,maxval,maxgrp
93 !el integer maxtyp,maxclass,maxkey
94 !el integer maxrot,maxopt
95 !el integer maxhess,maxlight,maxvib
96 !el integer maxgeo,maxcell,maxring
97 !el integer maxfix,maxbio
98 !el integer maxamino,maxnuc,maxbnd
99 !el integer maxang,maxtors,maxpi
100 !el integer maxpib,maxpit
101 ! integer :: maxatm !=2*nres !maxres2 maxres2=2*maxres
102 ! integer,parameter :: maxval=8
103 ! integer,parameter :: maxgrp=1000
104 ! integer,parameter :: maxtyp=3000
105 ! integer,parameter :: maxclass=500
106 ! integer,parameter :: maxkey=10000
107 ! integer,parameter :: maxrot=1000
108 ! integer,parameter :: maxopt=1000
109 ! integer,parameter :: maxhess=1000000
110 ! integer :: maxlight !=8*maxatm
111 ! integer,parameter :: maxvib=1000
112 ! integer,parameter :: maxgeo=1000
113 ! integer,parameter :: maxcell=10000
114 ! integer,parameter :: maxring=10000
115 ! integer,parameter :: maxfix=10000
116 ! integer,parameter :: maxbio=10000
117 ! integer,parameter :: maxamino=31
118 ! integer,parameter :: maxnuc=12
119 ! integer :: maxbnd !=2*maxatm
120 ! integer :: maxang !=3*maxatm
121 ! integer :: maxtors !=4*maxatm
122 ! integer,parameter :: maxpi=100
123 ! integer,parameter :: maxpib=2*maxpi
124 ! integer,parameter :: maxpit=4*maxpi
125 !-----------------------------------------------------------------------------
126 ! Maximum number of seed
127 ! integer,parameter :: max_seed=1
128 !-----------------------------------------------------------------------------
129 real(kind=8),dimension(:),allocatable :: stochforcvec !(MAXRES6) maxres6=6*maxres
130 ! common /stochcalc/ stochforcvec
131 !-----------------------------------------------------------------------------
132 ! common /przechowalnia/ subroutines: rattle1,rattle2,rattle_brown
133 real(kind=8),dimension(:,:),allocatable :: GGinv !(2*nres,2*nres) maxres2=2*maxres
134 real(kind=8),dimension(:,:,:),allocatable :: gdc !(3,2*nres,2*nres) maxres2=2*maxres
135 real(kind=8),dimension(:,:),allocatable :: Cmat !(2*nres,2*nres) maxres2=2*maxres
136 !-----------------------------------------------------------------------------
137 ! common /syfek/ subroutines: friction_force,setup_fricmat
138 !el real(kind=8),dimension(:),allocatable :: gamvec !(MAXRES6) or (MAXRES2)
139 !-----------------------------------------------------------------------------
140 ! common /przechowalnia/ subroutines: friction_force,setup_fricmat
141 real(kind=8),dimension(:,:),allocatable :: ginvfric !(2*nres,2*nres) !maxres2=2*maxres
142 !-----------------------------------------------------------------------------
143 ! common /przechowalnia/ subroutine: setup_fricmat
144 real(kind=8),dimension(:,:),allocatable :: fcopy !(2*nres,2*nres)
145 !-----------------------------------------------------------------------------
148 !-----------------------------------------------------------------------------
150 !-----------------------------------------------------------------------------
152 !-----------------------------------------------------------------------------
153 subroutine brown_step(itime)
154 !------------------------------------------------
155 ! Perform a single Euler integration step of Brownian dynamics
156 !------------------------------------------------
157 ! implicit real*8 (a-h,o-z)
159 use control, only: tcpu
162 ! use io_conf, only:cartprint
163 ! include 'DIMENSIONS'
167 ! include 'COMMON.CONTROL'
168 ! include 'COMMON.VAR'
169 ! include 'COMMON.MD'
171 ! include 'COMMON.LANGEVIN'
173 ! include 'COMMON.LANGEVIN.lang0'
175 ! include 'COMMON.CHAIN'
176 ! include 'COMMON.DERIV'
177 ! include 'COMMON.GEO'
178 ! include 'COMMON.LOCAL'
179 ! include 'COMMON.INTERACT'
180 ! include 'COMMON.IOUNITS'
181 ! include 'COMMON.NAMES'
182 ! include 'COMMON.TIME1'
183 real(kind=8),dimension(6*nres) :: zapas !(MAXRES6) maxres6=6*maxres
184 integer :: rstcount !ilen,
186 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
187 real(kind=8),dimension(6*nres,2*nres) :: Bmat,GBmat,Tmat !(MAXRES6,MAXRES2) (maxres2=2*maxres,maxres6=6*maxres)
188 real(kind=8),dimension(2*nres,2*nres) :: Cmat_,Cinv !(maxres2,maxres2) maxres2=2*maxres
189 real(kind=8),dimension(6*nres,6*nres) :: Pmat !(maxres6,maxres6) maxres6=6*maxres
190 ! real(kind=8),dimension(:,:),allocatable :: Bmat,GBmat,Tmat !(MAXRES6,MAXRES2) (maxres2=2*maxres,maxres6=6*maxres)
191 ! real(kind=8),dimension(:,:),allocatable :: Cmat_,Cinv !(maxres2,maxres2) maxres2=2*maxres
192 ! real(kind=8),dimension(:,:),allocatable :: Pmat !(maxres6,maxres6) maxres6=6*maxres
193 real(kind=8),dimension(6*nres) :: Td !(maxres6) maxres6=6*maxres
194 real(kind=8),dimension(2*nres) :: ppvec !(maxres2) maxres2=2*maxres
195 !el common /stochcalc/ stochforcvec
196 !el real(kind=8),dimension(3) :: cm !el
197 !el common /gucio/ cm
199 logical :: lprn = .false.,lprn1 = .false.
200 integer :: maxiter = 5
201 real(kind=8) :: difftol = 1.0d-5
202 real(kind=8) :: xx,diffmax,blen2,diffbond,tt0
203 integer :: i,j,nbond,k,ind,ind1,iter
204 integer :: nres2,nres6
208 ! if (.not.allocated(Bmat)) allocate(Bmat(nres6,nres2))
209 ! if (.not.allocated(GBmat)) allocate (GBmat(nres6,nres2))
210 ! if (.not.allocated(Tmat)) allocate (Tmat(nres6,nres2))
211 ! if (.not.allocated(Cmat_)) allocate(Cmat_(nres2,nres2))
212 ! if (.not.allocated(Cinv)) allocate (Cinv(nres2,nres2))
213 ! if (.not.allocated(Pmat)) allocate(Pmat(6*nres,6*nres))
215 if (.not.allocated(stochforcvec)) allocate(stochforcvec(nres6)) !(MAXRES6) maxres6=6*maxres
219 if (itype(i,1).ne.10) nbond=nbond+1
223 write (iout,*) "Generalized inverse of fricmat"
224 call matout(dimen,dimen,nres6,nres6,fricmat)
236 Bmat(ind+j,ind1)=dC_norm(j,i)
241 if (itype(i,1).ne.10) then
244 Bmat(ind+j,ind1)=dC_norm(j,i+nres)
250 write (iout,*) "Matrix Bmat"
251 call MATOUT(nbond,dimen,nres6,nres6,Bmat)
257 GBmat(i,j)=GBmat(i,j)+fricmat(i,k)*Bmat(k,j)
262 write (iout,*) "Matrix GBmat"
263 call MATOUT(nbond,dimen,nres6,nres2,Gbmat)
269 Cmat_(i,j)=Cmat_(i,j)+Bmat(k,i)*GBmat(k,j)
274 write (iout,*) "Matrix Cmat"
275 call MATOUT(nbond,nbond,nres2,nres2,Cmat_)
277 call matinvert(nbond,nres2,Cmat_,Cinv,osob)
279 write (iout,*) "Matrix Cinv"
280 call MATOUT(nbond,nbond,nres2,nres2,Cinv)
286 Tmat(i,j)=Tmat(i,j)+GBmat(i,k)*Cinv(k,j)
291 write (iout,*) "Matrix Tmat"
292 call MATOUT(nbond,dimen,nres6,nres2,Tmat)
302 Pmat(i,j)=Pmat(i,j)-Tmat(i,k)*Bmat(j,k)
307 write (iout,*) "Matrix Pmat"
308 call MATOUT(dimen,dimen,nres6,nres6,Pmat)
315 Td(i)=Td(i)+vbl*Tmat(i,ind)
318 if (itype(k,1).ne.10) then
320 Td(i)=Td(i)+vbldsc0(1,itype(k,1))*Tmat(i,ind)
325 write (iout,*) "Vector Td"
327 write (iout,'(i5,f10.5)') i,Td(i)
330 call stochastic_force(stochforcvec)
332 write (iout,*) "stochforcvec"
334 write (iout,*) i,stochforcvec(i)
338 zapas(j)=-gcart(j,0)+stochforcvec(j)
340 dC_work(j)=dC_old(j,0)
346 zapas(ind)=-gcart(j,i)+stochforcvec(ind)
347 dC_work(ind)=dC_old(j,i)
351 if (itype(i,1).ne.10) then
354 zapas(ind)=-gxcart(j,i)+stochforcvec(ind)
355 dC_work(ind)=dC_old(j,i+nres)
361 write (iout,*) "Initial d_t_work"
363 write (iout,*) i,d_t_work(i)
370 d_t_work(i)=d_t_work(i)+fricmat(i,j)*zapas(j)
377 zapas(i)=zapas(i)+Pmat(i,j)*(dC_work(j)+d_t_work(j)*d_time)
381 write (iout,*) "Final d_t_work and zapas"
383 write (iout,*) i,d_t_work(i),zapas(i)
397 dc_work(ind+j)=dc(j,i)
403 d_t(j,i+nres)=d_t_work(ind+j)
404 dc(j,i+nres)=zapas(ind+j)
405 dc_work(ind+j)=dc(j,i+nres)
411 write (iout,*) "Before correction for rotational lengthening"
412 write (iout,*) "New coordinates",&
413 " and differences between actual and standard bond lengths"
418 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
422 if (itype(i,1).ne.10) then
424 xx=vbld(i+nres)-vbldsc0(1,itype(i,1))
425 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
426 i,(dC(j,i+nres),j=1,3),xx
430 ! Second correction (rotational lengthening)
436 blen2 = scalar(dc(1,i),dc(1,i))
437 ppvec(ind)=2*vbl**2-blen2
438 diffbond=dabs(vbl-dsqrt(blen2))
439 if (diffbond.gt.diffmax) diffmax=diffbond
440 if (ppvec(ind).gt.0.0d0) then
441 ppvec(ind)=dsqrt(ppvec(ind))
446 write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind)
450 if (itype(i,1).ne.10) then
452 blen2 = scalar(dc(1,i+nres),dc(1,i+nres))
453 ppvec(ind)=2*vbldsc0(1,itype(i,1))**2-blen2
454 diffbond=dabs(vbldsc0(1,itype(i,1))-dsqrt(blen2))
455 if (diffbond.gt.diffmax) diffmax=diffbond
456 if (ppvec(ind).gt.0.0d0) then
457 ppvec(ind)=dsqrt(ppvec(ind))
462 write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind)
466 if (lprn) write (iout,*) "iter",iter," diffmax",diffmax
467 if (diffmax.lt.difftol) goto 10
471 Td(i)=Td(i)+ppvec(j)*Tmat(i,j)
477 zapas(i)=zapas(i)+Pmat(i,j)*dc_work(j)
488 dc_work(ind+j)=zapas(ind+j)
493 if (itype(i,1).ne.10) then
495 dc(j,i+nres)=zapas(ind+j)
496 dc_work(ind+j)=zapas(ind+j)
501 ! Building the chain from the newly calculated coordinates
504 if (large.and. mod(itime,ntwe).eq.0) then
505 write (iout,*) "Cartesian and internal coordinates: step 1"
508 write (iout,'(a)') "Potential forces"
510 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(-gcart(j,i),j=1,3),&
513 write (iout,'(a)') "Stochastic forces"
515 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(stochforc(j,i),j=1,3),&
516 (stochforc(j,i+nres),j=1,3)
518 write (iout,'(a)') "Velocities"
520 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
521 (d_t(j,i+nres),j=1,3)
526 write (iout,*) "After correction for rotational lengthening"
527 write (iout,*) "New coordinates",&
528 " and differences between actual and standard bond lengths"
533 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
537 if (itype(i,1).ne.10) then
539 xx=vbld(i+nres)-vbldsc0(1,itype(i,1))
540 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
541 i,(dC(j,i+nres),j=1,3),xx
546 ! write (iout,*) "Too many attempts at correcting the bonds"
554 ! Calculate energy and forces
556 call etotal(potEcomp)
557 potE=potEcomp(0)-potEcomp(20)
561 ! Calculate the kinetic and total energy and the kinetic temperature
564 t_enegrad=t_enegrad+MPI_Wtime()-tt0
566 t_enegrad=t_enegrad+tcpu()-tt0
569 kinetic_T=2.0d0/(dimen*Rb)*EK
571 end subroutine brown_step
572 !-----------------------------------------------------------------------------
574 !-----------------------------------------------------------------------------
575 subroutine gauss(RO,AP,MT,M,N,*)
577 ! CALCULATES (RO**(-1))*AP BY GAUSS ELIMINATION
578 ! RO IS A SQUARE MATRIX
579 ! THE CALCULATED PRODUCT IS STORED IN AP
580 ! ABNORMAL EXIT IF RO IS SINGULAR
582 integer :: MT, M, N, M1,I,J,IM,&
584 real(kind=8) :: RO(MT,M),AP(MT,N),X,RM,PR,Y
590 if(dabs(X).le.1.0D-13) return 1
602 if(DABS(RO(J,I)).LE.RM) goto 2
616 if(dabs(X).le.1.0E-13) return 1
625 8 AP(J,K)=AP(J,K)-Y*AP(I,K)
627 9 RO(J,K)=RO(J,K)-Y*RO(I,K)
631 if(dabs(X).le.1.0E-13) return 1
641 15 X=X-AP(K,J)*RO(MI,K)
646 !-----------------------------------------------------------------------------
648 !-----------------------------------------------------------------------------
649 subroutine kinetic(KE_total)
650 !----------------------------------------------------------------
651 ! This subroutine calculates the total kinetic energy of the chain
652 !-----------------------------------------------------------------
654 ! implicit real*8 (a-h,o-z)
655 ! include 'DIMENSIONS'
656 ! include 'COMMON.VAR'
657 ! include 'COMMON.CHAIN'
658 ! include 'COMMON.DERIV'
659 ! include 'COMMON.GEO'
660 ! include 'COMMON.LOCAL'
661 ! include 'COMMON.INTERACT'
662 ! include 'COMMON.MD'
663 ! include 'COMMON.IOUNITS'
664 real(kind=8) :: KE_total,mscab
666 integer :: i,j,k,iti,mnum
667 real(kind=8) :: KEt_p,KEt_sc,KEr_p,KEr_sc,incr(3),&
670 write (iout,*) "Velocities, kietic"
672 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
673 (d_t(j,i+nres),j=1,3)
678 ! write (iout,*) "ISC",(isc(itype(i,1)),i=1,nres)
679 ! The translational part for peptide virtual bonds
685 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
686 ! write (iout,*) "Kinetic trp:",i,(incr(j),j=1,3)
688 v(j)=incr(j)+0.5d0*d_t(j,i)
690 vtot(i)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
691 KEt_p=KEt_p+mp(mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
693 incr(j)=incr(j)+d_t(j,i)
696 ! write(iout,*) 'KEt_p', KEt_p
697 ! The translational part for the side chain virtual bond
698 ! Only now we can initialize incr with zeros. It must be equal
699 ! to the velocities of the first Calpha.
705 iti=iabs(itype(i,mnum))
711 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
712 if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
713 .or.(mnum.eq.5)) then
719 v(j)=incr(j)+d_t(j,nres+i)
722 ! write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3)
723 ! write (iout,*) "i",i," msc",msc(iti)," v",(v(j),j=1,3)
724 KEt_sc=KEt_sc+mscab*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
725 vtot(i+nres)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
727 incr(j)=incr(j)+d_t(j,i)
731 ! write(iout,*) 'KEt_sc', KEt_sc
732 ! The part due to stretching and rotation of the peptide groups
736 ! write (iout,*) "i",i
737 ! write (iout,*) "i",i," mag1",mag1," mag2",mag2
741 ! write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3)
742 KEr_p=KEr_p+Ip(mnum)*(incr(1)*incr(1)+incr(2)*incr(2) &
746 ! write(iout,*) 'KEr_p', KEr_p
747 ! The rotational part of the side chain virtual bond
751 iti=iabs(itype(i,mnum))
752 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
753 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
754 .and.(mnum.ne.5)) then
756 incr(j)=d_t(j,nres+i)
758 ! write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3)
759 KEr_sc=KEr_sc+Isc(iti,mnum)*(incr(1)*incr(1)+incr(2)*incr(2)+ &
763 ! The total kinetic energy
765 ! write(iout,*) 'KEr_sc', KEr_sc
766 KE_total=0.5d0*(KEt_p+KEt_sc+0.25d0*KEr_p+KEr_sc)
767 ! write (iout,*) "KE_total",KE_total
769 end subroutine kinetic
770 !-----------------------------------------------------------------------------
772 !-----------------------------------------------------------------------------
774 !------------------------------------------------
775 ! The driver for molecular dynamics subroutines
776 !------------------------------------------------
779 use control, only:tcpu,ovrtim
780 ! use io_comm, only:ilen
782 use compare, only:secondary2,hairpin
783 use io, only:cartout,statout
784 ! implicit real*8 (a-h,o-z)
785 ! include 'DIMENSIONS'
788 integer :: IERROR,ERRCODE
790 ! include 'COMMON.SETUP'
791 ! include 'COMMON.CONTROL'
792 ! include 'COMMON.VAR'
793 ! include 'COMMON.MD'
795 ! include 'COMMON.LANGEVIN'
797 ! include 'COMMON.LANGEVIN.lang0'
799 ! include 'COMMON.CHAIN'
800 ! include 'COMMON.DERIV'
801 ! include 'COMMON.GEO'
802 ! include 'COMMON.LOCAL'
803 ! include 'COMMON.INTERACT'
804 ! include 'COMMON.IOUNITS'
805 ! include 'COMMON.NAMES'
806 ! include 'COMMON.TIME1'
807 ! include 'COMMON.HAIRPIN'
808 real(kind=8),dimension(3) :: L,vcm
810 real(kind=8),dimension(6*nres) :: v_work,v_transf !(maxres6) maxres6=6*maxres
812 integer :: rstcount !ilen,
814 character(len=50) :: tytul
815 !el common /gucio/ cm
816 integer :: itime,i,j,nharp
817 integer,dimension(4,nres/3) :: iharp !(4,nres/3)(4,maxres/3)
819 real(kind=8) :: tt0,scalfac
825 print *,"MY tmpdir",tmpdir,ilen(tmpdir)
826 if (ilen(tmpdir).gt.0) &
827 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_" &
828 //liczba(:ilen(liczba))//'.rst')
830 if (ilen(tmpdir).gt.0) &
831 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"//'.rst')
838 write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
844 print *,"just befor setup matix",nres
845 ! Determine the inverse of the inertia matrix.
846 call setup_MD_matrices
848 print *,"AFTER SETUP MATRICES"
850 print *,"AFTER INIT MD"
853 t_MDsetup = MPI_Wtime()-tt0
855 t_MDsetup = tcpu()-tt0
858 ! Entering the MD loop
864 if (lang.eq.2 .or. lang.eq.3) then
868 call sd_verlet_p_setup
870 call sd_verlet_ciccotti_setup
874 pfric0_mat(i,j,0)=pfric_mat(i,j)
875 afric0_mat(i,j,0)=afric_mat(i,j)
876 vfric0_mat(i,j,0)=vfric_mat(i,j)
877 prand0_mat(i,j,0)=prand_mat(i,j)
878 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
879 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
884 flag_stoch(i)=.false.
888 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
890 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
894 else if (lang.eq.1 .or. lang.eq.4) then
895 print *,"before setup_fricmat"
897 print *,"after setup_fricmat"
900 t_langsetup=MPI_Wtime()-tt0
903 t_langsetup=tcpu()-tt0
906 do itime=1,n_timestep
908 if (large.and. mod(itime,ntwe).eq.0) &
909 write (iout,*) "itime",itime
911 if (lang.gt.0 .and. surfarea .and. &
912 mod(itime,reset_fricmat).eq.0) then
913 if (lang.eq.2 .or. lang.eq.3) then
917 call sd_verlet_p_setup
919 call sd_verlet_ciccotti_setup
923 pfric0_mat(i,j,0)=pfric_mat(i,j)
924 afric0_mat(i,j,0)=afric_mat(i,j)
925 vfric0_mat(i,j,0)=vfric_mat(i,j)
926 prand0_mat(i,j,0)=prand_mat(i,j)
927 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
928 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
933 flag_stoch(i)=.false.
936 else if (lang.eq.1 .or. lang.eq.4) then
939 write (iout,'(a,i10)') &
940 "Friction matrix reset based on surface area, itime",itime
942 if (reset_vel .and. tbf .and. lang.eq.0 &
943 .and. mod(itime,count_reset_vel).eq.0) then
945 write(iout,'(a,f20.2)') &
946 "Velocities reset to random values, time",totT
949 d_t_old(j,i)=d_t(j,i)
953 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
957 d_t(j,0)=d_t(j,0)-vcm(j)
960 kinetic_T=2.0d0/(dimen3*Rb)*EK
961 scalfac=dsqrt(T_bath/kinetic_T)
962 write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
965 d_t_old(j,i)=scalfac*d_t(j,i)
971 ! Time-reversible RESPA algorithm
972 ! (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
973 call RESPA_step(itime)
975 ! Variable time step algorithm.
976 call velverlet_step(itime)
980 call brown_step(itime)
982 print *,"Brown dynamics not here!"
984 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
990 if (mod(itime,ntwe).eq.0) then
993 ! call check_ecartint
1003 v_work(ind)=d_t(j,i)
1008 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1.and.mnum.ne.5) then
1011 v_work(ind)=d_t(j,i+nres)
1016 write (66,'(80f10.5)') &
1017 ((d_t(j,i),j=1,3),i=0,nres-1),((d_t(j,i+nres),j=1,3),i=1,nres)
1021 v_transf(i)=v_transf(i)+gvec(j,i)*v_work(j)
1023 v_transf(i)= v_transf(i)*dsqrt(geigen(i))
1025 write (67,'(80f10.5)') (v_transf(i),i=1,ind)
1028 if (mod(itime,ntwx).eq.0) then
1030 write (tytul,'("time",f8.2)') totT
1032 call hairpin(.true.,nharp,iharp)
1033 call secondary2(.true.)
1034 call pdbout(potE,tytul,ipdb)
1039 if (rstcount.eq.1000.or.itime.eq.n_timestep) then
1040 open(irest2,file=rest2name,status='unknown')
1041 write(irest2,*) totT,EK,potE,totE,t_bath
1043 ! AL 4/17/17: Now writing d_t(0,:) too
1045 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
1047 ! AL 4/17/17: Now writing d_c(0,:) too
1049 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
1057 t_MD=MPI_Wtime()-tt0
1061 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))') &
1063 'MD calculations setup:',t_MDsetup,&
1064 'Energy & gradient evaluation:',t_enegrad,&
1065 'Stochastic MD setup:',t_langsetup,&
1066 'Stochastic MD step setup:',t_sdsetup,&
1068 write (iout,'(/28(1h=),a25,27(1h=))') &
1069 ' End of MD calculation '
1071 write (iout,*) "time for etotal",t_etotal," elong",t_elong,&
1073 write (iout,*) "time_fric",time_fric," time_stoch",time_stoch,&
1074 " time_fricmatmult",time_fricmatmult," time_fsample ",&
1079 !-----------------------------------------------------------------------------
1080 subroutine velverlet_step(itime)
1081 !-------------------------------------------------------------------------------
1082 ! Perform a single velocity Verlet step; the time step can be rescaled if
1083 ! increments in accelerations exceed the threshold
1084 !-------------------------------------------------------------------------------
1085 ! implicit real*8 (a-h,o-z)
1086 ! include 'DIMENSIONS'
1088 use control, only:tcpu
1090 use minimm, only:minim_dc
1093 integer :: ierror,ierrcode
1094 real(kind=8) :: errcode
1096 ! include 'COMMON.SETUP'
1097 ! include 'COMMON.VAR'
1098 ! include 'COMMON.MD'
1100 ! include 'COMMON.LANGEVIN'
1102 ! include 'COMMON.LANGEVIN.lang0'
1104 ! include 'COMMON.CHAIN'
1105 ! include 'COMMON.DERIV'
1106 ! include 'COMMON.GEO'
1107 ! include 'COMMON.LOCAL'
1108 ! include 'COMMON.INTERACT'
1109 ! include 'COMMON.IOUNITS'
1110 ! include 'COMMON.NAMES'
1111 ! include 'COMMON.TIME1'
1112 ! include 'COMMON.MUCA'
1113 real(kind=8),dimension(3) :: vcm,incr
1114 real(kind=8),dimension(3) :: L
1115 integer :: count,rstcount !ilen,
1117 character(len=50) :: tytul
1118 integer :: maxcount_scale = 30
1119 !el common /gucio/ cm
1120 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1121 !el common /stochcalc/ stochforcvec
1122 integer :: itime,icount_scale,itime_scal,i,j,ifac_time,iretcode
1124 real(kind=8) :: epdrift,tt0,fac_time
1126 if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres)) !(MAXRES6) maxres6=6*maxres
1132 else if (lang.eq.2 .or. lang.eq.3) then
1134 call stochastic_force(stochforcvec)
1137 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
1139 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1146 icount_scale=icount_scale+1
1147 ! write(iout,*) "icount_scale",icount_scale,scalel
1148 if (icount_scale.gt.maxcount_scale) then
1150 "ERROR: too many attempts at scaling down the time step. ",&
1151 "amax=",amax,"epdrift=",epdrift,&
1152 "damax=",damax,"edriftmax=",edriftmax,&
1156 call MPI_Abort(MPI_COMM_WORLD,IERROR,IERRCODE)
1160 ! First step of the velocity Verlet algorithm
1165 else if (lang.eq.3) then
1167 call sd_verlet1_ciccotti
1169 else if (lang.eq.1) then
1174 ! Build the chain from the newly calculated coordinates
1175 call chainbuild_cart
1176 if (rattle) call rattle1
1178 if (large.and. mod(itime,ntwe).eq.0) then
1179 write (iout,*) "Cartesian and internal coordinates: step 1"
1184 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),&
1185 (dc(j,i+nres),j=1,3)
1187 write (iout,*) "Accelerations"
1189 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1190 (d_a(j,i+nres),j=1,3)
1192 write (iout,*) "Velocities, step 1"
1194 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1195 (d_t(j,i+nres),j=1,3)
1204 ! Calculate energy and forces
1206 call etotal(potEcomp)
1207 ! AL 4/17/17: Reduce the steps if NaNs occurred.
1208 if (potEcomp(0).gt.0.99e18 .or. isnan(potEcomp(0)).gt.0) then
1210 if (icount_scale.gt.15) then
1211 write (iout,*) "Tu jest problem",potEcomp(0),d_time
1212 ! call gen_rand_conf(1,*335)
1213 ! call minim_dc(potEcomp(0),iretcode,100)
1216 ! call etotal(potEcomp)
1217 ! write(iout,*) "needed to repara,",potEcomp
1220 ! 335 write(iout,*) "Failed genrand"
1224 if (large.and. mod(itime,ntwe).eq.0) &
1225 call enerprint(potEcomp)
1228 t_etotal=t_etotal+MPI_Wtime()-tt0
1230 t_etotal=t_etotal+tcpu()-tt0
1233 potE=potEcomp(0)-potEcomp(20)
1235 ! Get the new accelerations
1238 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1240 t_enegrad=t_enegrad+tcpu()-tt0
1242 ! Determine maximum acceleration and scale down the timestep if needed
1244 amax=amax/(itime_scal+1)**2
1245 call predict_edrift(epdrift)
1246 ! write(iout,*) "amax=",amax,damax,epdrift,edriftmax,amax/(itime_scal+1)
1248 ! write (iout,*) "before enter if",scalel,icount_scale
1249 if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then
1250 ! write(iout,*) "I enter if"
1251 ! Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
1253 ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax)) &
1255 itime_scal=itime_scal+ifac_time
1256 ! fac_time=dmin1(damax/amax,0.5d0)
1257 fac_time=0.5d0**ifac_time
1258 d_time=d_time*fac_time
1259 if (lang.eq.2 .or. lang.eq.3) then
1261 ! write (iout,*) "Calling sd_verlet_setup: 1"
1262 ! Rescale the stochastic forces and recalculate or restore
1263 ! the matrices of tinker integrator
1264 if (itime_scal.gt.maxflag_stoch) then
1265 if (large) write (iout,'(a,i5,a)') &
1266 "Calculate matrices for stochastic step;",&
1267 " itime_scal ",itime_scal
1269 call sd_verlet_p_setup
1271 call sd_verlet_ciccotti_setup
1273 write (iout,'(2a,i3,a,i3,1h.)') &
1274 "Warning: cannot store matrices for stochastic",&
1275 " integration because the index",itime_scal,&
1276 " is greater than",maxflag_stoch
1277 write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",&
1278 " integration Langevin algorithm for better efficiency."
1279 else if (flag_stoch(itime_scal)) then
1280 if (large) write (iout,'(a,i5,a,l1)') &
1281 "Restore matrices for stochastic step; itime_scal ",&
1282 itime_scal," flag ",flag_stoch(itime_scal)
1285 pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
1286 afric_mat(i,j)=afric0_mat(i,j,itime_scal)
1287 vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
1288 prand_mat(i,j)=prand0_mat(i,j,itime_scal)
1289 vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
1290 vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
1294 if (large) write (iout,'(2a,i5,a,l1)') &
1295 "Calculate & store matrices for stochastic step;",&
1296 " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
1298 call sd_verlet_p_setup
1300 call sd_verlet_ciccotti_setup
1302 flag_stoch(ifac_time)=.true.
1305 pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
1306 afric0_mat(i,j,itime_scal)=afric_mat(i,j)
1307 vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
1308 prand0_mat(i,j,itime_scal)=prand_mat(i,j)
1309 vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
1310 vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
1314 fac_time=1.0d0/dsqrt(fac_time)
1316 stochforcvec(i)=fac_time*stochforcvec(i)
1319 else if (lang.eq.1) then
1320 ! Rescale the accelerations due to stochastic forces
1321 fac_time=1.0d0/dsqrt(fac_time)
1323 d_as_work(i)=d_as_work(i)*fac_time
1326 if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)') &
1327 "itime",itime," Timestep scaled down to ",&
1328 d_time," ifac_time",ifac_time," itime_scal",itime_scal
1330 ! Second step of the velocity Verlet algorithm
1335 else if (lang.eq.3) then
1337 call sd_verlet2_ciccotti
1339 else if (lang.eq.1) then
1344 if (rattle) call rattle2
1347 if (d_time.ne.d_time0) then
1350 if (lang.eq.2 .or. lang.eq.3) then
1351 if (large) write (iout,'(a)') &
1352 "Restore original matrices for stochastic step"
1353 ! write (iout,*) "Calling sd_verlet_setup: 2"
1354 ! Restore the matrices of tinker integrator if the time step has been restored
1357 pfric_mat(i,j)=pfric0_mat(i,j,0)
1358 afric_mat(i,j)=afric0_mat(i,j,0)
1359 vfric_mat(i,j)=vfric0_mat(i,j,0)
1360 prand_mat(i,j)=prand0_mat(i,j,0)
1361 vrand_mat1(i,j)=vrand0_mat1(i,j,0)
1362 vrand_mat2(i,j)=vrand0_mat2(i,j,0)
1370 ! Calculate the kinetic and the total energy and the kinetic temperature
1374 ! call kinetic1(EK1)
1375 ! write (iout,*) "step",itime," EK",EK," EK1",EK1
1377 ! Couple the system to Berendsen bath if needed
1378 if (tbf .and. lang.eq.0) then
1381 kinetic_T=2.0d0/(dimen3*Rb)*EK
1382 ! Backup the coordinates, velocities, and accelerations
1386 d_t_old(j,i)=d_t(j,i)
1387 d_a_old(j,i)=d_a(j,i)
1391 if (mod(itime,ntwe).eq.0 .and. large) then
1392 write (iout,*) "Velocities, step 2"
1394 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1395 (d_t(j,i+nres),j=1,3)
1400 end subroutine velverlet_step
1401 !-----------------------------------------------------------------------------
1402 subroutine RESPA_step(itime)
1403 !-------------------------------------------------------------------------------
1404 ! Perform a single RESPA step.
1405 !-------------------------------------------------------------------------------
1406 ! implicit real*8 (a-h,o-z)
1407 ! include 'DIMENSIONS'
1411 use control, only:tcpu
1413 ! use io_conf, only:cartprint
1416 integer :: IERROR,ERRCODE
1418 ! include 'COMMON.SETUP'
1419 ! include 'COMMON.CONTROL'
1420 ! include 'COMMON.VAR'
1421 ! include 'COMMON.MD'
1423 ! include 'COMMON.LANGEVIN'
1425 ! include 'COMMON.LANGEVIN.lang0'
1427 ! include 'COMMON.CHAIN'
1428 ! include 'COMMON.DERIV'
1429 ! include 'COMMON.GEO'
1430 ! include 'COMMON.LOCAL'
1431 ! include 'COMMON.INTERACT'
1432 ! include 'COMMON.IOUNITS'
1433 ! include 'COMMON.NAMES'
1434 ! include 'COMMON.TIME1'
1435 real(kind=8),dimension(0:n_ene) :: energia_short,energia_long
1436 real(kind=8),dimension(3) :: L,vcm,incr
1437 real(kind=8),dimension(3,0:2*nres) :: dc_old0,d_t_old0,d_a_old0 !(3,0:maxres2) maxres2=2*maxres
1438 logical :: PRINT_AMTS_MSG = .false.
1439 integer :: count,rstcount !ilen,
1441 character(len=50) :: tytul
1442 integer :: maxcount_scale = 10
1443 !el common /gucio/ cm
1444 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1445 !el common /stochcalc/ stochforcvec
1446 integer :: itime,itt,i,j,itsplit
1448 !el common /cipiszcze/ itt
1450 real(kind=8) :: epdrift,tt0,epdriftmax
1453 if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres)) !(MAXRES6) maxres6=6*maxres
1457 if (large.and. mod(itime,ntwe).eq.0) then
1458 write (iout,*) "***************** RESPA itime",itime
1459 write (iout,*) "Cartesian and internal coordinates: step 0"
1461 call pdbout(0.0d0,"cipiszcze",iout)
1463 write (iout,*) "Accelerations from long-range forces"
1465 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1466 (d_a(j,i+nres),j=1,3)
1468 write (iout,*) "Velocities, step 0"
1470 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1471 (d_t(j,i+nres),j=1,3)
1476 ! Perform the initial RESPA step (increment velocities)
1477 ! write (iout,*) "*********************** RESPA ini"
1480 if (mod(itime,ntwe).eq.0 .and. large) then
1481 write (iout,*) "Velocities, end"
1483 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1484 (d_t(j,i+nres),j=1,3)
1488 ! Compute the short-range forces
1494 ! 7/2/2009 commented out
1496 ! call etotal_short(energia_short)
1499 ! 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
1502 d_a(j,i)=d_a_short(j,i)
1506 if (large.and. mod(itime,ntwe).eq.0) then
1507 write (iout,*) "energia_short",energia_short(0)
1508 write (iout,*) "Accelerations from short-range forces"
1510 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1511 (d_a(j,i+nres),j=1,3)
1516 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1518 t_enegrad=t_enegrad+tcpu()-tt0
1523 d_t_old(j,i)=d_t(j,i)
1524 d_a_old(j,i)=d_a(j,i)
1527 ! 6/30/08 A-MTS: attempt at increasing the split number
1530 dc_old0(j,i)=dc_old(j,i)
1531 d_t_old0(j,i)=d_t_old(j,i)
1532 d_a_old0(j,i)=d_a_old(j,i)
1535 if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
1536 if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
1543 ! write (iout,*) "itime",itime," ntime_split",ntime_split
1544 ! Split the time step
1545 d_time=d_time0/ntime_split
1546 ! Perform the short-range RESPA steps (velocity Verlet increments of
1547 ! positions and velocities using short-range forces)
1548 ! write (iout,*) "*********************** RESPA split"
1549 do itsplit=1,ntime_split
1552 else if (lang.eq.2 .or. lang.eq.3) then
1554 call stochastic_force(stochforcvec)
1557 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
1559 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1564 ! First step of the velocity Verlet algorithm
1569 else if (lang.eq.3) then
1571 call sd_verlet1_ciccotti
1573 else if (lang.eq.1) then
1578 ! Build the chain from the newly calculated coordinates
1579 call chainbuild_cart
1580 if (rattle) call rattle1
1582 if (large.and. mod(itime,ntwe).eq.0) then
1583 write (iout,*) "***** ITSPLIT",itsplit
1584 write (iout,*) "Cartesian and internal coordinates: step 1"
1585 call pdbout(0.0d0,"cipiszcze",iout)
1588 write (iout,*) "Velocities, step 1"
1590 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1591 (d_t(j,i+nres),j=1,3)
1600 ! Calculate energy and forces
1602 call etotal_short(energia_short)
1603 ! AL 4/17/17: Exit itime_split loop when energy goes infinite
1604 if (energia_short(0).gt.0.99e20 .or. isnan(energia_short(0)) ) then
1605 if (PRINT_AMTS_MSG) &
1606 write (iout,*) "Infinities/NaNs in energia_short",energia_short(0),"; increasing ntime_split to",ntime_split
1607 ntime_split=ntime_split*2
1608 if (ntime_split.gt.maxtime_split) then
1611 "Cannot rescue the run; aborting job. Retry with a smaller time step"
1613 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1616 "Cannot rescue the run; terminating. Retry with a smaller time step"
1622 if (large.and. mod(itime,ntwe).eq.0) &
1623 call enerprint(energia_short)
1626 t_eshort=t_eshort+MPI_Wtime()-tt0
1628 t_eshort=t_eshort+tcpu()-tt0
1632 ! Get the new accelerations
1634 ! 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
1637 d_a_short(j,i)=d_a(j,i)
1641 if (large.and. mod(itime,ntwe).eq.0) then
1642 write (iout,*)"energia_short",energia_short(0)
1643 write (iout,*) "Accelerations from short-range forces"
1645 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1646 (d_a(j,i+nres),j=1,3)
1651 ! Determine maximum acceleration and scale down the timestep if needed
1653 amax=amax/ntime_split**2
1654 call predict_edrift(epdrift)
1655 if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0) &
1656 write (iout,*) "amax",amax," damax",damax,&
1657 " epdrift",epdrift," epdriftmax",epdriftmax
1658 ! Exit loop and try with increased split number if the change of
1659 ! acceleration is too big
1660 if (amax.gt.damax .or. epdrift.gt.edriftmax) then
1661 if (ntime_split.lt.maxtime_split) then
1663 ntime_split=ntime_split*2
1664 ! AL 4/17/17: We should exit the itime_split loop when acceleration change is too big
1668 dc_old(j,i)=dc_old0(j,i)
1669 d_t_old(j,i)=d_t_old0(j,i)
1670 d_a_old(j,i)=d_a_old0(j,i)
1673 if (PRINT_AMTS_MSG) then
1674 write (iout,*) "acceleration/energy drift too large",amax,&
1675 epdrift," split increased to ",ntime_split," itime",itime,&
1681 "Uh-hu. Bumpy landscape. Maximum splitting number",&
1683 " already reached!!! Trying to carry on!"
1687 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1689 t_enegrad=t_enegrad+tcpu()-tt0
1691 ! Second step of the velocity Verlet algorithm
1696 else if (lang.eq.3) then
1698 call sd_verlet2_ciccotti
1700 else if (lang.eq.1) then
1705 if (rattle) call rattle2
1706 ! Backup the coordinates, velocities, and accelerations
1710 d_t_old(j,i)=d_t(j,i)
1711 d_a_old(j,i)=d_a(j,i)
1718 ! Restore the time step
1720 ! Compute long-range forces
1727 call etotal_long(energia_long)
1728 if (energia_long(0).gt.0.99e20 .or. isnan(energia_long(0))) then
1731 "Infinitied/NaNs in energia_long, Aborting MPI job."
1733 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1735 write (iout,*) "Infinitied/NaNs in energia_long, terminating."
1739 if (large.and. mod(itime,ntwe).eq.0) &
1740 call enerprint(energia_long)
1743 t_elong=t_elong+MPI_Wtime()-tt0
1745 t_elong=t_elong+tcpu()-tt0
1751 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1753 t_enegrad=t_enegrad+tcpu()-tt0
1755 ! Compute accelerations from long-range forces
1757 if (large.and. mod(itime,ntwe).eq.0) then
1758 write (iout,*) "energia_long",energia_long(0)
1759 write (iout,*) "Cartesian and internal coordinates: step 2"
1761 call pdbout(0.0d0,"cipiszcze",iout)
1763 write (iout,*) "Accelerations from long-range forces"
1765 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1766 (d_a(j,i+nres),j=1,3)
1768 write (iout,*) "Velocities, step 2"
1770 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1771 (d_t(j,i+nres),j=1,3)
1775 ! Compute the final RESPA step (increment velocities)
1776 ! write (iout,*) "*********************** RESPA fin"
1778 ! Compute the complete potential energy
1780 potEcomp(i)=energia_short(i)+energia_long(i)
1782 potE=potEcomp(0)-potEcomp(20)
1783 ! potE=energia_short(0)+energia_long(0)
1786 ! Calculate the kinetic and the total energy and the kinetic temperature
1789 ! Couple the system to Berendsen bath if needed
1790 if (tbf .and. lang.eq.0) then
1793 kinetic_T=2.0d0/(dimen3*Rb)*EK
1794 ! Backup the coordinates, velocities, and accelerations
1796 if (mod(itime,ntwe).eq.0 .and. large) then
1797 write (iout,*) "Velocities, end"
1799 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1800 (d_t(j,i+nres),j=1,3)
1805 end subroutine RESPA_step
1806 !-----------------------------------------------------------------------------
1807 subroutine RESPA_vel
1808 ! First and last RESPA step (incrementing velocities using long-range
1811 ! implicit real*8 (a-h,o-z)
1812 ! include 'DIMENSIONS'
1813 ! include 'COMMON.CONTROL'
1814 ! include 'COMMON.VAR'
1815 ! include 'COMMON.MD'
1816 ! include 'COMMON.CHAIN'
1817 ! include 'COMMON.DERIV'
1818 ! include 'COMMON.GEO'
1819 ! include 'COMMON.LOCAL'
1820 ! include 'COMMON.INTERACT'
1821 ! include 'COMMON.IOUNITS'
1822 ! include 'COMMON.NAMES'
1823 integer :: i,j,inres,mnum
1826 d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
1830 d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
1835 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
1836 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
1837 .and.(mnum.ne.5)) then
1840 d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
1845 end subroutine RESPA_vel
1846 !-----------------------------------------------------------------------------
1848 ! Applying velocity Verlet algorithm - step 1 to coordinates
1850 ! implicit real*8 (a-h,o-z)
1851 ! include 'DIMENSIONS'
1852 ! include 'COMMON.CONTROL'
1853 ! include 'COMMON.VAR'
1854 ! include 'COMMON.MD'
1855 ! include 'COMMON.CHAIN'
1856 ! include 'COMMON.DERIV'
1857 ! include 'COMMON.GEO'
1858 ! include 'COMMON.LOCAL'
1859 ! include 'COMMON.INTERACT'
1860 ! include 'COMMON.IOUNITS'
1861 ! include 'COMMON.NAMES'
1862 real(kind=8) :: adt,adt2
1863 integer :: i,j,inres,mnum
1866 write (iout,*) "VELVERLET1 START: DC"
1868 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),&
1869 (dc(j,i+nres),j=1,3)
1873 adt=d_a_old(j,0)*d_time
1875 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1876 d_t_new(j,0)=d_t_old(j,0)+adt2
1877 d_t(j,0)=d_t_old(j,0)+adt
1881 adt=d_a_old(j,i)*d_time
1883 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1884 d_t_new(j,i)=d_t_old(j,i)+adt2
1885 d_t(j,i)=d_t_old(j,i)+adt
1890 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
1891 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
1892 .and.(mnum.ne.5)) then
1895 adt=d_a_old(j,inres)*d_time
1897 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1898 d_t_new(j,inres)=d_t_old(j,inres)+adt2
1899 d_t(j,inres)=d_t_old(j,inres)+adt
1904 write (iout,*) "VELVERLET1 END: DC"
1906 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),&
1907 (dc(j,i+nres),j=1,3)
1911 end subroutine verlet1
1912 !-----------------------------------------------------------------------------
1914 ! Step 2 of the velocity Verlet algorithm: update velocities
1916 ! implicit real*8 (a-h,o-z)
1917 ! include 'DIMENSIONS'
1918 ! include 'COMMON.CONTROL'
1919 ! include 'COMMON.VAR'
1920 ! include 'COMMON.MD'
1921 ! include 'COMMON.CHAIN'
1922 ! include 'COMMON.DERIV'
1923 ! include 'COMMON.GEO'
1924 ! include 'COMMON.LOCAL'
1925 ! include 'COMMON.INTERACT'
1926 ! include 'COMMON.IOUNITS'
1927 ! include 'COMMON.NAMES'
1928 integer :: i,j,inres,mnum
1931 d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
1935 d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
1940 ! iti=iabs(itype(i,mnum))
1941 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
1942 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
1943 .and.(mnum.ne.5)) then
1946 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
1951 end subroutine verlet2
1952 !-----------------------------------------------------------------------------
1953 subroutine sddir_precalc
1954 ! Applying velocity Verlet algorithm - step 1 to coordinates
1955 ! implicit real*8 (a-h,o-z)
1956 ! include 'DIMENSIONS'
1962 ! include 'COMMON.CONTROL'
1963 ! include 'COMMON.VAR'
1964 ! include 'COMMON.MD'
1966 ! include 'COMMON.LANGEVIN'
1968 ! include 'COMMON.LANGEVIN.lang0'
1970 ! include 'COMMON.CHAIN'
1971 ! include 'COMMON.DERIV'
1972 ! include 'COMMON.GEO'
1973 ! include 'COMMON.LOCAL'
1974 ! include 'COMMON.INTERACT'
1975 ! include 'COMMON.IOUNITS'
1976 ! include 'COMMON.NAMES'
1977 ! include 'COMMON.TIME1'
1978 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1979 !el common /stochcalc/ stochforcvec
1980 real(kind=8) :: time00
1982 ! Compute friction and stochastic forces
1987 time_fric=time_fric+MPI_Wtime()-time00
1989 call stochastic_force(stochforcvec)
1990 time_stoch=time_stoch+MPI_Wtime()-time00
1993 ! Compute the acceleration due to friction forces (d_af_work) and stochastic
1994 ! forces (d_as_work)
1996 call ginv_mult(fric_work, d_af_work)
1997 call ginv_mult(stochforcvec, d_as_work)
1999 end subroutine sddir_precalc
2000 !-----------------------------------------------------------------------------
2001 subroutine sddir_verlet1
2002 ! Applying velocity Verlet algorithm - step 1 to velocities
2005 ! implicit real*8 (a-h,o-z)
2006 ! include 'DIMENSIONS'
2007 ! include 'COMMON.CONTROL'
2008 ! include 'COMMON.VAR'
2009 ! include 'COMMON.MD'
2011 ! include 'COMMON.LANGEVIN'
2013 ! include 'COMMON.LANGEVIN.lang0'
2015 ! include 'COMMON.CHAIN'
2016 ! include 'COMMON.DERIV'
2017 ! include 'COMMON.GEO'
2018 ! include 'COMMON.LOCAL'
2019 ! include 'COMMON.INTERACT'
2020 ! include 'COMMON.IOUNITS'
2021 ! include 'COMMON.NAMES'
2022 ! Revised 3/31/05 AL: correlation between random contributions to
2023 ! position and velocity increments included.
2024 real(kind=8) :: sqrt13 = 0.57735026918962576451d0 ! 1/sqrt(3)
2025 real(kind=8) :: adt,adt2
2026 integer :: i,j,ind,inres,mnum
2028 ! Add the contribution from BOTH friction and stochastic force to the
2029 ! coordinates, but ONLY the contribution from the friction forces to velocities
2032 adt=(d_a_old(j,0)+d_af_work(j))*d_time
2033 adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
2034 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
2035 d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
2036 d_t(j,0)=d_t_old(j,0)+adt
2041 adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
2042 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
2043 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
2044 d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
2045 d_t(j,i)=d_t_old(j,i)+adt
2051 ! iti=iabs(itype(i,mnum))
2052 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2053 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2054 .and.(mnum.ne.5)) then
2057 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
2058 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
2059 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
2060 d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
2061 d_t(j,inres)=d_t_old(j,inres)+adt
2067 end subroutine sddir_verlet1
2068 !-----------------------------------------------------------------------------
2069 subroutine sddir_verlet2
2070 ! Calculating the adjusted velocities for accelerations
2073 ! implicit real*8 (a-h,o-z)
2074 ! include 'DIMENSIONS'
2075 ! include 'COMMON.CONTROL'
2076 ! include 'COMMON.VAR'
2077 ! include 'COMMON.MD'
2079 ! include 'COMMON.LANGEVIN'
2081 ! include 'COMMON.LANGEVIN.lang0'
2083 ! include 'COMMON.CHAIN'
2084 ! include 'COMMON.DERIV'
2085 ! include 'COMMON.GEO'
2086 ! include 'COMMON.LOCAL'
2087 ! include 'COMMON.INTERACT'
2088 ! include 'COMMON.IOUNITS'
2089 ! include 'COMMON.NAMES'
2090 real(kind=8),dimension(6*nres) :: stochforcvec,d_as_work1 !(MAXRES6) maxres6=6*maxres
2091 real(kind=8) :: cos60 = 0.5d0, sin60 = 0.86602540378443864676d0
2092 integer :: i,j,ind,inres,mnum
2093 ! Revised 3/31/05 AL: correlation between random contributions to
2094 ! position and velocity increments included.
2095 ! The correlation coefficients are calculated at low-friction limit.
2096 ! Also, friction forces are now not calculated with new velocities.
2098 ! call friction_force
2099 call stochastic_force(stochforcvec)
2101 ! Compute the acceleration due to friction forces (d_af_work) and stochastic
2102 ! forces (d_as_work)
2104 call ginv_mult(stochforcvec, d_as_work1)
2110 d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j)) &
2111 +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
2116 d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j)) &
2117 +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
2123 ! iti=iabs(itype(i,mnum))
2124 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2125 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2126 .and.(mnum.ne.5)) then
2129 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres) &
2130 +d_af_work(ind+j))+sin60*d_as_work(ind+j) &
2131 +cos60*d_as_work1(ind+j))*d_time
2137 end subroutine sddir_verlet2
2138 !-----------------------------------------------------------------------------
2139 subroutine max_accel
2141 ! Find the maximum difference in the accelerations of the the sites
2142 ! at the beginning and the end of the time step.
2145 ! implicit real*8 (a-h,o-z)
2146 ! include 'DIMENSIONS'
2147 ! include 'COMMON.CONTROL'
2148 ! include 'COMMON.VAR'
2149 ! include 'COMMON.MD'
2150 ! include 'COMMON.CHAIN'
2151 ! include 'COMMON.DERIV'
2152 ! include 'COMMON.GEO'
2153 ! include 'COMMON.LOCAL'
2154 ! include 'COMMON.INTERACT'
2155 ! include 'COMMON.IOUNITS'
2156 real(kind=8),dimension(3) :: aux,accel,accel_old
2157 real(kind=8) :: dacc
2161 ! aux(j)=d_a(j,0)-d_a_old(j,0)
2162 accel_old(j)=d_a_old(j,0)
2169 ! 7/3/08 changed to asymmetric difference
2171 ! accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
2172 accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
2173 accel(j)=accel(j)+0.5d0*d_a(j,i)
2174 ! if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
2175 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
2176 dacc=dabs(accel(j)-accel_old(j))
2177 ! write (iout,*) i,dacc
2178 if (dacc.gt.amax) amax=dacc
2186 accel_old(j)=d_a_old(j,0)
2191 accel_old(j)=accel_old(j)+d_a_old(j,1)
2192 accel(j)=accel(j)+d_a(j,1)
2197 ! iti=iabs(itype(i,mnum))
2198 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2199 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2200 .and.(mnum.ne.5)) then
2202 ! accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
2203 accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
2204 accel(j)=accel(j)+d_a(j,i+nres)
2208 ! if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
2209 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
2210 dacc=dabs(accel(j)-accel_old(j))
2211 ! write (iout,*) "side-chain",i,dacc
2212 if (dacc.gt.amax) amax=dacc
2216 accel_old(j)=accel_old(j)+d_a_old(j,i)
2217 accel(j)=accel(j)+d_a(j,i)
2218 ! aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
2222 end subroutine max_accel
2223 !-----------------------------------------------------------------------------
2224 subroutine predict_edrift(epdrift)
2226 ! Predict the drift of the potential energy
2229 use control_data, only: lmuca
2230 ! implicit real*8 (a-h,o-z)
2231 ! include 'DIMENSIONS'
2232 ! include 'COMMON.CONTROL'
2233 ! include 'COMMON.VAR'
2234 ! include 'COMMON.MD'
2235 ! include 'COMMON.CHAIN'
2236 ! include 'COMMON.DERIV'
2237 ! include 'COMMON.GEO'
2238 ! include 'COMMON.LOCAL'
2239 ! include 'COMMON.INTERACT'
2240 ! include 'COMMON.IOUNITS'
2241 ! include 'COMMON.MUCA'
2242 real(kind=8) :: epdrift,epdriftij
2244 ! Drift of the potential energy
2250 epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
2251 if (lmuca) epdriftij=epdriftij*factor
2252 ! write (iout,*) "back",i,j,epdriftij
2253 if (epdriftij.gt.epdrift) epdrift=epdriftij
2257 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1.and.&
2258 molnum(i).ne.5) then
2261 dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
2262 if (lmuca) epdriftij=epdriftij*factor
2263 ! write (iout,*) "side",i,j,epdriftij
2264 if (epdriftij.gt.epdrift) epdrift=epdriftij
2268 epdrift=0.5d0*epdrift*d_time*d_time
2269 ! write (iout,*) "epdrift",epdrift
2271 end subroutine predict_edrift
2272 !-----------------------------------------------------------------------------
2273 subroutine verlet_bath
2275 ! Coupling to the thermostat by using the Berendsen algorithm
2278 ! implicit real*8 (a-h,o-z)
2279 ! include 'DIMENSIONS'
2280 ! include 'COMMON.CONTROL'
2281 ! include 'COMMON.VAR'
2282 ! include 'COMMON.MD'
2283 ! include 'COMMON.CHAIN'
2284 ! include 'COMMON.DERIV'
2285 ! include 'COMMON.GEO'
2286 ! include 'COMMON.LOCAL'
2287 ! include 'COMMON.INTERACT'
2288 ! include 'COMMON.IOUNITS'
2289 ! include 'COMMON.NAMES'
2290 real(kind=8) :: T_half,fact
2291 integer :: i,j,inres,mnum
2293 T_half=2.0d0/(dimen3*Rb)*EK
2294 fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
2295 ! write(iout,*) "T_half", T_half
2296 ! write(iout,*) "EK", EK
2297 ! write(iout,*) "fact", fact
2299 d_t(j,0)=fact*d_t(j,0)
2303 d_t(j,i)=fact*d_t(j,i)
2308 ! iti=iabs(itype(i,mnum))
2309 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2310 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2311 .and.(mnum.ne.5)) then
2314 d_t(j,inres)=fact*d_t(j,inres)
2319 end subroutine verlet_bath
2320 !-----------------------------------------------------------------------------
2322 ! Set up the initial conditions of a MD simulation
2325 use control, only:tcpu
2326 !el use io_basic, only:ilen
2329 use minimm, only:minim_dc,minimize,sc_move
2330 use io_config, only:readrst
2331 use io, only:statout
2332 ! implicit real*8 (a-h,o-z)
2333 ! include 'DIMENSIONS'
2336 character(len=16) :: form
2337 integer :: IERROR,ERRCODE
2339 integer :: iranmin,itrial,itmp
2340 ! include 'COMMON.SETUP'
2341 ! include 'COMMON.CONTROL'
2342 ! include 'COMMON.VAR'
2343 ! include 'COMMON.MD'
2345 ! include 'COMMON.LANGEVIN'
2347 ! include 'COMMON.LANGEVIN.lang0'
2349 ! include 'COMMON.CHAIN'
2350 ! include 'COMMON.DERIV'
2351 ! include 'COMMON.GEO'
2352 ! include 'COMMON.LOCAL'
2353 ! include 'COMMON.INTERACT'
2354 ! include 'COMMON.IOUNITS'
2355 ! include 'COMMON.NAMES'
2356 ! include 'COMMON.REMD'
2357 real(kind=8),dimension(0:n_ene) :: energia_long,energia_short
2358 real(kind=8),dimension(3) :: vcm,incr,L
2359 real(kind=8) :: xv,sigv,lowb,highb
2360 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
2361 character(len=256) :: qstr
2364 character(len=50) :: tytul
2365 logical :: file_exist
2366 !el common /gucio/ cm
2367 integer :: i,j,ipos,iq,iw,nft_sc,iretcode,nfun,itime,ierr,mnum
2368 real(kind=8) :: etot,tt0
2372 ! write(iout,*) "d_time", d_time
2373 ! Compute the standard deviations of stochastic forces for Langevin dynamics
2374 ! if the friction coefficients do not depend on surface area
2375 if (lang.gt.0 .and. .not.surfarea) then
2378 stdforcp(i)=stdfp(mnum)*dsqrt(gamp(mnum))
2382 stdforcsc(i)=stdfsc(iabs(itype(i,mnum)),mnum) &
2383 *dsqrt(gamsc(iabs(itype(i,mnum)),mnum))
2386 ! Open the pdb file for snapshotshots
2389 if (ilen(tmpdir).gt.0) &
2390 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2391 liczba(:ilen(liczba))//".pdb")
2393 file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2397 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) &
2398 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2399 liczba(:ilen(liczba))//".x")
2400 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2403 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) &
2404 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2405 liczba(:ilen(liczba))//".cx")
2406 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2412 if (ilen(tmpdir).gt.0) &
2413 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
2414 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
2416 if (ilen(tmpdir).gt.0) &
2417 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
2418 cartname=prefix(:ilen(prefix))//"_MD.cx"
2422 write (qstr,'(256(1h ))')
2425 iq = qinfrag(i,iset)*10
2426 iw = wfrag(i,iset)/100
2428 if(me.eq.king.or..not.out1file) &
2429 write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
2430 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
2435 iq = qinpair(i,iset)*10
2436 iw = wpair(i,iset)/100
2438 if(me.eq.king.or..not.out1file) &
2439 write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
2440 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
2444 ! pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
2446 ! cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
2448 ! cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
2450 ! statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
2454 if (restart1file) then
2456 inquire(file=mremd_rst_name,exist=file_exist)
2457 write (*,*) me," Before broadcast: file_exist",file_exist
2459 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,&
2462 write (*,*) me," After broadcast: file_exist",file_exist
2463 ! inquire(file=mremd_rst_name,exist=file_exist)
2464 if(me.eq.king.or..not.out1file) &
2465 write(iout,*) "Initial state read by master and distributed"
2467 if (ilen(tmpdir).gt.0) &
2468 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_' &
2469 //liczba(:ilen(liczba))//'.rst')
2470 inquire(file=rest2name,exist=file_exist)
2473 if(.not.restart1file) then
2474 if(me.eq.king.or..not.out1file) &
2475 write(iout,*) "Initial state will be read from file ",&
2476 rest2name(:ilen(rest2name))
2479 call rescale_weights(t_bath)
2481 if(me.eq.king.or..not.out1file)then
2482 if (restart1file) then
2483 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),&
2486 write(iout,*) "File ",rest2name(:ilen(rest2name)),&
2489 write(iout,*) "Initial velocities randomly generated"
2496 ! Generate initial velocities
2497 if(me.eq.king.or..not.out1file) &
2498 write(iout,*) "Initial velocities randomly generated"
2503 ! rest2name = prefix(:ilen(prefix))//'.rst'
2504 if(me.eq.king.or..not.out1file)then
2505 write (iout,*) "Initial velocities"
2507 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2508 (d_t(j,i+nres),j=1,3)
2510 ! Zeroing the total angular momentum of the system
2511 write(iout,*) "Calling the zero-angular momentum subroutine"
2514 ! Getting the potential energy and forces and velocities and accelerations
2516 ! write (iout,*) "velocity of the center of the mass:"
2517 ! write (iout,*) (vcm(j),j=1,3)
2519 d_t(j,0)=d_t(j,0)-vcm(j)
2521 ! Removing the velocity of the center of mass
2523 if(me.eq.king.or..not.out1file)then
2524 write (iout,*) "vcm right after adjustment:"
2525 write (iout,*) (vcm(j),j=1,3)
2529 if(iranconf.ne.0 .or.indpdb.gt.0.and..not.unres_pdb .or.preminim) then
2531 print *, 'Calling OVERLAP_SC'
2532 call overlap_sc(fail)
2533 print *,'after OVERLAP'
2536 print *,'call SC_MOVE'
2537 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
2538 print *,'SC_move',nft_sc,etot
2539 if(me.eq.king.or..not.out1file) &
2540 write(iout,*) 'SC_move',nft_sc,etot
2544 print *, 'Calling MINIM_DC'
2545 call minim_dc(etot,iretcode,nfun)
2547 call geom_to_var(nvar,varia)
2548 print *,'Calling MINIMIZE.'
2549 call minimize(etot,varia,iretcode,nfun)
2550 call var_to_geom(nvar,varia)
2552 if(me.eq.king.or..not.out1file) &
2553 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2555 if(iranconf.ne.0) then
2556 !c 8/22/17 AL Loop to produce a low-energy random conformation
2559 if(me.eq.king.or..not.out1file) &
2560 write (iout,*) 'Calling OVERLAP_SC'
2561 call overlap_sc(fail)
2562 endif !endif overlap
2565 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
2566 print *,'SC_move',nft_sc,etot
2567 if(me.eq.king.or..not.out1file) &
2568 write(iout,*) 'SC_move',nft_sc,etot
2572 print *, 'Calling MINIM_DC'
2573 call minim_dc(etot,iretcode,nfun)
2574 call int_from_cart1(.false.)
2576 call geom_to_var(nvar,varia)
2577 print *,'Calling MINIMIZE.'
2578 call minimize(etot,varia,iretcode,nfun)
2579 call var_to_geom(nvar,varia)
2581 if(me.eq.king.or..not.out1file) &
2582 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2584 if (isnan(etot) .or. etot.gt.4.0d6) then
2585 write (iout,*) "Energy too large",etot, &
2586 " trying another random conformation"
2589 call gen_rand_conf(itmp,*30)
2591 30 write (iout,*) 'Failed to generate random conformation', &
2593 write (*,*) 'Processor:',me, &
2594 ' Failed to generate random conformation',&
2603 write (iout,'(a,i3,a)') 'Processor:',me, &
2604 ' error in generating random conformation.'
2605 write (*,'(a,i3,a)') 'Processor:',me, &
2606 ' error in generating random conformation.'
2609 ! call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
2610 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
2620 write (iout,'(a,i3,a)') 'Processor:',me, &
2621 ' failed to generate a low-energy random conformation.'
2622 write (*,'(a,i3,a,f10.3)') 'Processor:',me, &
2623 ' failed to generate a low-energy random conformation.',etot
2627 ! call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
2628 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
2635 call chainbuild_cart
2640 kinetic_T=2.0d0/(dimen3*Rb)*EK
2641 if(me.eq.king.or..not.out1file)then
2651 write(iout,*) "before ETOTAL"
2652 call etotal(potEcomp)
2653 if (large) call enerprint(potEcomp)
2656 t_etotal=t_etotal+MPI_Wtime()-tt0
2658 t_etotal=t_etotal+tcpu()-tt0
2665 if (amax*d_time .gt. dvmax) then
2666 d_time=d_time*dvmax/amax
2667 if(me.eq.king.or..not.out1file) write (iout,*) &
2668 "Time step reduced to",d_time,&
2669 " because of too large initial acceleration."
2671 if(me.eq.king.or..not.out1file)then
2672 write(iout,*) "Potential energy and its components"
2673 call enerprint(potEcomp)
2674 ! write(iout,*) (potEcomp(i),i=0,n_ene)
2676 potE=potEcomp(0)-potEcomp(20)
2679 if (ntwe.ne.0) call statout(itime)
2680 if(me.eq.king.or..not.out1file) &
2681 write (iout,'(/a/3(a25,1pe14.5/))') "Initial:", &
2682 " Kinetic energy",EK," Potential energy",potE, &
2683 " Total energy",totE," Maximum acceleration ", &
2686 write (iout,*) "Initial coordinates"
2688 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),&
2691 write (iout,*) "Initial dC"
2693 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),&
2694 (dc(j,i+nres),j=1,3)
2696 write (iout,*) "Initial velocities"
2697 write (iout,"(13x,' backbone ',23x,' side chain')")
2699 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2700 (d_t(j,i+nres),j=1,3)
2702 write (iout,*) "Initial accelerations"
2704 ! write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2705 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),&
2706 (d_a(j,i+nres),j=1,3)
2712 d_t_old(j,i)=d_t(j,i)
2713 d_a_old(j,i)=d_a(j,i)
2715 ! write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
2724 call etotal_short(energia_short)
2725 if (large) call enerprint(potEcomp)
2728 t_eshort=t_eshort+MPI_Wtime()-tt0
2730 t_eshort=t_eshort+tcpu()-tt0
2735 if(.not.out1file .and. large) then
2736 write (iout,*) "energia_long",energia_long(0),&
2737 " energia_short",energia_short(0),&
2738 " total",energia_long(0)+energia_short(0)
2739 write (iout,*) "Initial fast-force accelerations"
2741 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
2742 (d_a(j,i+nres),j=1,3)
2745 ! 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
2748 d_a_short(j,i)=d_a(j,i)
2757 call etotal_long(energia_long)
2758 if (large) call enerprint(potEcomp)
2761 t_elong=t_elong+MPI_Wtime()-tt0
2763 t_elong=t_elong+tcpu()-tt0
2768 if(.not.out1file .and. large) then
2769 write (iout,*) "energia_long",energia_long(0)
2770 write (iout,*) "Initial slow-force accelerations"
2772 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
2773 (d_a(j,i+nres),j=1,3)
2777 t_enegrad=t_enegrad+MPI_Wtime()-tt0
2779 t_enegrad=t_enegrad+tcpu()-tt0
2783 end subroutine init_MD
2784 !-----------------------------------------------------------------------------
2785 subroutine random_vel
2787 ! implicit real*8 (a-h,o-z)
2789 use random, only:anorm_distr
2791 ! include 'DIMENSIONS'
2792 ! include 'COMMON.CONTROL'
2793 ! include 'COMMON.VAR'
2794 ! include 'COMMON.MD'
2796 ! include 'COMMON.LANGEVIN'
2798 ! include 'COMMON.LANGEVIN.lang0'
2800 ! include 'COMMON.CHAIN'
2801 ! include 'COMMON.DERIV'
2802 ! include 'COMMON.GEO'
2803 ! include 'COMMON.LOCAL'
2804 ! include 'COMMON.INTERACT'
2805 ! include 'COMMON.IOUNITS'
2806 ! include 'COMMON.NAMES'
2807 ! include 'COMMON.TIME1'
2808 real(kind=8) :: xv,sigv,lowb,highb ,Ek1
2811 real(kind=8) ,allocatable, dimension(:) :: DDU1,DDU2,DL2,DL1,xsolv,DML,rs
2812 real(kind=8) :: sumx
2814 real(kind=8) ,allocatable, dimension(:) :: rsold
2815 real (kind=8),allocatable,dimension(:,:) :: matold
2819 integer :: i,j,ii,k,ind,mark,imark,mnum
2820 ! Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
2821 ! First generate velocities in the eigenspace of the G matrix
2822 ! write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
2829 sigv=dsqrt((Rb*t_bath)/geigen(i))
2832 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
2834 write (iout,*) "i",i," ii",ii," geigen",geigen(i),&
2835 " d_t_work_new",d_t_work_new(ii)
2846 Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
2849 write (iout,*) "Ek from eigenvectors",Ek1
2850 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
2854 ! Transform velocities to UNRES coordinate space
2855 allocate (DL1(2*nres))
2856 allocate (DDU1(2*nres))
2857 allocate (DL2(2*nres))
2858 allocate (DDU2(2*nres))
2859 allocate (xsolv(2*nres))
2860 allocate (DML(2*nres))
2861 allocate (rs(2*nres))
2863 allocate (rsold(2*nres))
2864 allocate (matold(2*nres,2*nres))
2866 matold(1,1)=DMorig(1)
2867 matold(1,2)=DU1orig(1)
2868 matold(1,3)=DU2orig(1)
2869 write (*,*) DMorig(1),DU1orig(1),DU2orig(1)
2874 matold(i,j)=DMorig(i)
2875 matold(i,j-1)=DU1orig(i-1)
2877 matold(i,j-2)=DU2orig(i-2)
2885 matold(i,j+1)=DU1orig(i)
2891 matold(i,j+2)=DU2orig(i)
2895 matold(dimen,dimen)=DMorig(dimen)
2896 matold(dimen,dimen-1)=DU1orig(dimen-1)
2897 matold(dimen,dimen-2)=DU2orig(dimen-2)
2898 write(iout,*) "old gmatrix"
2899 call matout(dimen,dimen,2*nres,2*nres,matold)
2903 ! Find the ith eigenvector of the pentadiagonal inertiq matrix
2907 DML(j)=DMorig(j)-geigen(i)
2910 DML(j-1)=DMorig(j)-geigen(i)
2915 DDU1(imark-1)=DU2orig(imark-1)
2916 do j=imark+1,dimen-1
2917 DDU1(j-1)=DU1orig(j)
2925 DDU2(j)=DU2orig(j+1)
2934 write (iout,*) "DL2,DL1,DML,DDU1,DDU2"
2935 write(iout,'(10f10.5)') (DL2(k),k=3,dimen-1)
2936 write(iout,'(10f10.5)') (DL1(k),k=2,dimen-1)
2937 write(iout,'(10f10.5)') (DML(k),k=1,dimen-1)
2938 write(iout,'(10f10.5)') (DDU1(k),k=1,dimen-2)
2939 write(iout,'(10f10.5)') (DDU2(k),k=1,dimen-3)
2942 if (imark.gt.2) rs(imark-2)=-DU2orig(imark-2)
2943 if (imark.gt.1) rs(imark-1)=-DU1orig(imark-1)
2944 if (imark.lt.dimen) rs(imark)=-DU1orig(imark)
2945 if (imark.lt.dimen-1) rs(imark+1)=-DU2orig(imark)
2949 ! write (iout,*) "Vector rs"
2951 ! write (iout,*) j,rs(j)
2954 call FDIAG(dimen-1,DL2,DL1,DML,DDU1,DDU2,rs,xsolv,mark)
2961 sumx=-geigen(i)*xsolv(j)
2963 sumx=sumx+matold(j,k)*xsolv(k)
2966 sumx=sumx+matold(j,k)*xsolv(k-1)
2968 write(iout,'(i5,3f10.5)') j,sumx,rsold(j),sumx-rsold(j)
2971 sumx=-geigen(i)*xsolv(j-1)
2973 sumx=sumx+matold(j,k)*xsolv(k)
2976 sumx=sumx+matold(j,k)*xsolv(k-1)
2978 write(iout,'(i5,3f10.5)') j-1,sumx,rsold(j-1),sumx-rsold(j-1)
2982 "Solution of equations system",i," for eigenvalue",geigen(i)
2984 write(iout,'(i5,f10.5)') j,xsolv(j)
2987 do j=dimen-1,imark,-1
2992 write (iout,*) "Un-normalized eigenvector",i," for eigenvalue",geigen(i)
2994 write(iout,'(i5,f10.5)') j,xsolv(j)
2997 ! Normalize ith eigenvector
3000 sumx=sumx+xsolv(j)**2
3004 xsolv(j)=xsolv(j)/sumx
3007 write (iout,*) "Eigenvector",i," for eigenvalue",geigen(i)
3009 write(iout,'(i5,3f10.5)') j,xsolv(j)
3012 ! All done at this point for eigenvector i, exit loop
3020 write (iout,*) "Unable to find eigenvector",i
3023 ! write (iout,*) "i=",i
3025 ! write (iout,*) "k=",k
3028 ! write(iout,*) "ind",ind," ind1",3*(i-1)+k
3029 d_t_work(ind)=d_t_work(ind) &
3030 +xsolv(j)*d_t_work_new(3*(i-1)+k)
3033 enddo ! i (loop over the eigenvectors)
3036 write (iout,*) "d_t_work"
3038 write (iout,"(i5,f10.5)") i,d_t_work(i)
3043 ! if (itype(i,1).eq.10) then
3045 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
3046 iti=iabs(itype(i,mnum))
3047 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3048 if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
3049 .or.(mnum.eq.5)) then
3056 ! write (iout,*) "k",k," ii+k",ii+k," ii+j+k",ii+j+k,"EK1",Ek1
3057 Ek1=Ek1+0.5d0*mp(mnum)*((d_t_work(ii+j+k)+d_t_work_new(ii+k))/2)**2+&
3058 0.5d0*0.25d0*IP(mnum)*(d_t_work(ii+j+k)-d_t_work_new(ii+k))**2
3061 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3062 .and.(mnum.ne.5)) ii=ii+3
3063 write (iout,*) "i",i," itype",itype(i,mnum)," mass",msc(itype(i,mnum),mnum)
3064 write (iout,*) "ii",ii
3067 write (iout,*) "k",k," ii",ii,"EK1",EK1
3068 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3070 Ek1=Ek1+0.5d0*Isc(iabs(itype(i,mnum)),mnum)*(d_t_work(ii)-d_t_work(ii-3))**2
3071 Ek1=Ek1+0.5d0*msc(iabs(itype(i,mnum)),mnum)*d_t_work(ii)**2
3073 write (iout,*) "i",i," ii",ii
3075 write (iout,*) "Ek from d_t_work",Ek1
3076 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
3078 if(allocated(DDU1)) deallocate(DDU1)
3079 if(allocated(DDU2)) deallocate(DDU2)
3080 if(allocated(DL2)) deallocate(DL2)
3081 if(allocated(DL1)) deallocate(DL1)
3082 if(allocated(xsolv)) deallocate(xsolv)
3083 if(allocated(DML)) deallocate(DML)
3084 if(allocated(rs)) deallocate(rs)
3086 if(allocated(matold)) deallocate(matold)
3087 if(allocated(rsold)) deallocate(rsold)
3092 d_t(k,j)=d_t_work(ind)
3096 if (itype(j,1).ne.10 .and. itype(j,mnum).ne.ntyp1_molec(mnum)&
3097 .and.(mnum.ne.5)) then
3099 d_t(k,j+nres)=d_t_work(ind)
3105 write (iout,*) "Random velocities in the Calpha,SC space"
3107 write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3)
3110 write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3)
3117 ! if (itype(i,1).eq.10) then
3119 if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
3120 .or.(mnum.eq.5)) then
3122 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
3126 d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
3127 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
3132 write (iout,*)"Random velocities in the virtual-bond-vector space"
3134 write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3)
3137 write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3)
3140 write (iout,*) "Ek from d_t_work",Ek1
3141 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
3149 d_t_work(ind)=d_t_work(ind) &
3150 +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
3152 ! write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
3156 ! Transfer to the d_t vector
3158 d_t(j,0)=d_t_work(j)
3164 d_t(j,i)=d_t_work(ind)
3169 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3170 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3171 .and.(mnum.ne.5)) then
3174 d_t(j,i+nres)=d_t_work(ind)
3180 ! write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",&
3181 ! 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
3183 ! write(iout,*) "end init MD"
3185 end subroutine random_vel
3186 !-----------------------------------------------------------------------------
3188 subroutine sd_verlet_p_setup
3189 ! Sets up the parameters of stochastic Verlet algorithm
3190 ! implicit real*8 (a-h,o-z)
3191 ! include 'DIMENSIONS'
3192 use control, only: tcpu
3197 ! include 'COMMON.CONTROL'
3198 ! include 'COMMON.VAR'
3199 ! include 'COMMON.MD'
3201 ! include 'COMMON.LANGEVIN'
3203 ! include 'COMMON.LANGEVIN.lang0'
3205 ! include 'COMMON.CHAIN'
3206 ! include 'COMMON.DERIV'
3207 ! include 'COMMON.GEO'
3208 ! include 'COMMON.LOCAL'
3209 ! include 'COMMON.INTERACT'
3210 ! include 'COMMON.IOUNITS'
3211 ! include 'COMMON.NAMES'
3212 ! include 'COMMON.TIME1'
3213 real(kind=8),dimension(6*nres) :: emgdt !(MAXRES6) maxres6=6*maxres
3214 real(kind=8) :: pterm,vterm,rho,rhoc,vsig
3215 real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
3216 prand_vec,vrand_vec1,vrand_vec2 !(MAXRES6) maxres6=6*maxres
3217 logical :: lprn = .false.
3218 real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
3219 real(kind=8) :: ktm,gdt,egdt,gdt2,gdt3,gdt4,gdt5,gdt6,gdt7,gdt8,&
3221 integer :: i,maxres2
3228 ! AL 8/17/04 Code adapted from tinker
3230 ! Get the frictional and random terms for stochastic dynamics in the
3231 ! eigenspace of mass-scaled UNRES friction matrix
3235 gdt = fricgam(i) * d_time
3237 ! Stochastic dynamics reduces to simple MD for zero friction
3239 if (gdt .le. zero) then
3240 pfric_vec(i) = 1.0d0
3241 vfric_vec(i) = d_time
3242 afric_vec(i) = 0.5d0 * d_time * d_time
3243 prand_vec(i) = 0.0d0
3244 vrand_vec1(i) = 0.0d0
3245 vrand_vec2(i) = 0.0d0
3247 ! Analytical expressions when friction coefficient is large
3250 if (gdt .ge. gdt_radius) then
3253 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
3254 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
3255 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
3256 vterm = 1.0d0 - egdt**2
3257 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
3259 ! Use series expansions when friction coefficient is small
3270 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0 &
3271 - gdt5/120.0d0 + gdt6/720.0d0 &
3272 - gdt7/5040.0d0 + gdt8/40320.0d0 &
3273 - gdt9/362880.0d0) / fricgam(i)**2
3274 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
3275 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
3276 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0 &
3277 + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0 &
3278 + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0 &
3279 + 127.0d0*gdt9/90720.0d0
3280 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0 &
3281 - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0 &
3282 - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0 &
3283 - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
3284 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0 &
3285 - 17.0d0*gdt2/1280.0d0 &
3286 + 17.0d0*gdt3/6144.0d0 &
3287 + 40967.0d0*gdt4/34406400.0d0 &
3288 - 57203.0d0*gdt5/275251200.0d0 &
3289 - 1429487.0d0*gdt6/13212057600.0d0)
3292 ! Compute the scaling factors of random terms for the nonzero friction case
3294 ktm = 0.5d0*d_time/fricgam(i)
3295 psig = dsqrt(ktm*pterm) / fricgam(i)
3296 vsig = dsqrt(ktm*vterm)
3297 rhoc = dsqrt(1.0d0 - rho*rho)
3299 vrand_vec1(i) = vsig * rho
3300 vrand_vec2(i) = vsig * rhoc
3305 "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",&
3308 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),&
3309 afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
3313 ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
3316 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
3317 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
3318 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
3319 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
3320 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
3321 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
3324 t_sdsetup=t_sdsetup+MPI_Wtime()
3326 t_sdsetup=t_sdsetup+tcpu()-tt0
3329 end subroutine sd_verlet_p_setup
3330 !-----------------------------------------------------------------------------
3331 subroutine eigtransf1(n,ndim,ab,d,c)
3335 real(kind=8) :: ab(ndim,ndim,n),c(ndim,n),d(ndim)
3341 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
3346 end subroutine eigtransf1
3347 !-----------------------------------------------------------------------------
3348 subroutine eigtransf(n,ndim,a,b,d,c)
3352 real(kind=8) :: a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
3358 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
3363 end subroutine eigtransf
3364 !-----------------------------------------------------------------------------
3365 subroutine sd_verlet1
3367 ! Applying stochastic velocity Verlet algorithm - step 1 to velocities
3369 ! implicit real*8 (a-h,o-z)
3370 ! include 'DIMENSIONS'
3371 ! include 'COMMON.CONTROL'
3372 ! include 'COMMON.VAR'
3373 ! include 'COMMON.MD'
3375 ! include 'COMMON.LANGEVIN'
3377 ! include 'COMMON.LANGEVIN.lang0'
3379 ! include 'COMMON.CHAIN'
3380 ! include 'COMMON.DERIV'
3381 ! include 'COMMON.GEO'
3382 ! include 'COMMON.LOCAL'
3383 ! include 'COMMON.INTERACT'
3384 ! include 'COMMON.IOUNITS'
3385 ! include 'COMMON.NAMES'
3386 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
3387 !el common /stochcalc/ stochforcvec
3388 logical :: lprn = .false.
3389 real(kind=8) :: ddt1,ddt2
3390 integer :: i,j,ind,inres
3392 ! write (iout,*) "dc_old"
3394 ! write (iout,'(i5,3f10.5,5x,3f10.5)')
3395 ! & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3398 dc_work(j)=dc_old(j,0)
3399 d_t_work(j)=d_t_old(j,0)
3400 d_a_work(j)=d_a_old(j,0)
3405 dc_work(ind+j)=dc_old(j,i)
3406 d_t_work(ind+j)=d_t_old(j,i)
3407 d_a_work(ind+j)=d_a_old(j,i)
3413 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3414 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3415 .and.(mnum.ne.5)) then
3417 dc_work(ind+j)=dc_old(j,i+nres)
3418 d_t_work(ind+j)=d_t_old(j,i+nres)
3419 d_a_work(ind+j)=d_a_old(j,i+nres)
3427 "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",&
3431 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
3432 vfric_mat(i,j),afric_mat(i,j),&
3433 prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3441 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) &
3442 +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3443 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3444 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3446 d_t_work_new(i)=ddt1+0.5d0*ddt2
3447 d_t_work(i)=ddt1+ddt2
3452 d_t(j,0)=d_t_work(j)
3457 dc(j,i)=dc_work(ind+j)
3458 d_t(j,i)=d_t_work(ind+j)
3464 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3465 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3466 .and.(mnum.ne.5)) then
3469 dc(j,inres)=dc_work(ind+j)
3470 d_t(j,inres)=d_t_work(ind+j)
3476 end subroutine sd_verlet1
3477 !-----------------------------------------------------------------------------
3478 subroutine sd_verlet2
3480 ! Calculating the adjusted velocities for accelerations
3482 ! implicit real*8 (a-h,o-z)
3483 ! include 'DIMENSIONS'
3484 ! include 'COMMON.CONTROL'
3485 ! include 'COMMON.VAR'
3486 ! include 'COMMON.MD'
3488 ! include 'COMMON.LANGEVIN'
3490 ! include 'COMMON.LANGEVIN.lang0'
3492 ! include 'COMMON.CHAIN'
3493 ! include 'COMMON.DERIV'
3494 ! include 'COMMON.GEO'
3495 ! include 'COMMON.LOCAL'
3496 ! include 'COMMON.INTERACT'
3497 ! include 'COMMON.IOUNITS'
3498 ! include 'COMMON.NAMES'
3499 !el real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV !(MAXRES6) maxres6=6*maxres
3500 real(kind=8),dimension(6*nres) :: stochforcvecV !(MAXRES6) maxres6=6*maxres
3501 !el common /stochcalc/ stochforcvec
3503 real(kind=8) :: ddt1,ddt2
3504 integer :: i,j,ind,inres
3505 ! Compute the stochastic forces which contribute to velocity change
3507 call stochastic_force(stochforcvecV)
3514 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3515 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+ &
3516 vrand_mat2(i,j)*stochforcvecV(j)
3518 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
3522 d_t(j,0)=d_t_work(j)
3527 d_t(j,i)=d_t_work(ind+j)
3533 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3534 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3535 .and.(mnum.ne.5)) then
3538 d_t(j,inres)=d_t_work(ind+j)
3544 end subroutine sd_verlet2
3545 !-----------------------------------------------------------------------------
3546 subroutine sd_verlet_ciccotti_setup
3548 ! Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
3550 ! implicit real*8 (a-h,o-z)
3551 ! include 'DIMENSIONS'
3552 use control, only: tcpu
3557 ! include 'COMMON.CONTROL'
3558 ! include 'COMMON.VAR'
3559 ! include 'COMMON.MD'
3561 ! include 'COMMON.LANGEVIN'
3563 ! include 'COMMON.LANGEVIN.lang0'
3565 ! include 'COMMON.CHAIN'
3566 ! include 'COMMON.DERIV'
3567 ! include 'COMMON.GEO'
3568 ! include 'COMMON.LOCAL'
3569 ! include 'COMMON.INTERACT'
3570 ! include 'COMMON.IOUNITS'
3571 ! include 'COMMON.NAMES'
3572 ! include 'COMMON.TIME1'
3573 real(kind=8),dimension(6*nres) :: emgdt !(MAXRES6) maxres6=6*maxres
3574 real(kind=8) :: pterm,vterm,rho,rhoc,vsig
3575 real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
3576 prand_vec,vrand_vec1,vrand_vec2 !(MAXRES6) maxres6=6*maxres
3577 logical :: lprn = .false.
3578 real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
3579 real(kind=8) :: ktm,gdt,egdt,tt0
3580 integer :: i,maxres2
3587 ! AL 8/17/04 Code adapted from tinker
3589 ! Get the frictional and random terms for stochastic dynamics in the
3590 ! eigenspace of mass-scaled UNRES friction matrix
3594 write (iout,*) "i",i," fricgam",fricgam(i)
3595 gdt = fricgam(i) * d_time
3597 ! Stochastic dynamics reduces to simple MD for zero friction
3599 if (gdt .le. zero) then
3600 pfric_vec(i) = 1.0d0
3601 vfric_vec(i) = d_time
3602 afric_vec(i) = 0.5d0*d_time*d_time
3603 prand_vec(i) = afric_vec(i)
3604 vrand_vec2(i) = vfric_vec(i)
3606 ! Analytical expressions when friction coefficient is large
3611 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
3612 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
3613 prand_vec(i) = afric_vec(i)
3614 vrand_vec2(i) = vfric_vec(i)
3616 ! Compute the scaling factors of random terms for the nonzero friction case
3618 ! ktm = 0.5d0*d_time/fricgam(i)
3619 ! psig = dsqrt(ktm*pterm) / fricgam(i)
3620 ! vsig = dsqrt(ktm*vterm)
3621 ! prand_vec(i) = psig*afric_vec(i)
3622 ! vrand_vec2(i) = vsig*vfric_vec(i)
3627 "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",&
3630 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),&
3631 afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
3635 ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
3637 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
3638 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
3639 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
3640 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
3641 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
3643 t_sdsetup=t_sdsetup+MPI_Wtime()
3645 t_sdsetup=t_sdsetup+tcpu()-tt0
3648 end subroutine sd_verlet_ciccotti_setup
3649 !-----------------------------------------------------------------------------
3650 subroutine sd_verlet1_ciccotti
3652 ! Applying stochastic velocity Verlet algorithm - step 1 to velocities
3653 ! implicit real*8 (a-h,o-z)
3655 ! include 'DIMENSIONS'
3659 ! include 'COMMON.CONTROL'
3660 ! include 'COMMON.VAR'
3661 ! include 'COMMON.MD'
3663 ! include 'COMMON.LANGEVIN'
3665 ! include 'COMMON.LANGEVIN.lang0'
3667 ! include 'COMMON.CHAIN'
3668 ! include 'COMMON.DERIV'
3669 ! include 'COMMON.GEO'
3670 ! include 'COMMON.LOCAL'
3671 ! include 'COMMON.INTERACT'
3672 ! include 'COMMON.IOUNITS'
3673 ! include 'COMMON.NAMES'
3674 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
3675 !el common /stochcalc/ stochforcvec
3676 logical :: lprn = .false.
3677 real(kind=8) :: ddt1,ddt2
3678 integer :: i,j,ind,inres
3679 ! write (iout,*) "dc_old"
3681 ! write (iout,'(i5,3f10.5,5x,3f10.5)')
3682 ! & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3685 dc_work(j)=dc_old(j,0)
3686 d_t_work(j)=d_t_old(j,0)
3687 d_a_work(j)=d_a_old(j,0)
3692 dc_work(ind+j)=dc_old(j,i)
3693 d_t_work(ind+j)=d_t_old(j,i)
3694 d_a_work(ind+j)=d_a_old(j,i)
3699 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3701 dc_work(ind+j)=dc_old(j,i+nres)
3702 d_t_work(ind+j)=d_t_old(j,i+nres)
3703 d_a_work(ind+j)=d_a_old(j,i+nres)
3712 "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",&
3716 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
3717 vfric_mat(i,j),afric_mat(i,j),&
3718 prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3726 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) &
3727 +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3728 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3729 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3731 d_t_work_new(i)=ddt1+0.5d0*ddt2
3732 d_t_work(i)=ddt1+ddt2
3737 d_t(j,0)=d_t_work(j)
3742 dc(j,i)=dc_work(ind+j)
3743 d_t(j,i)=d_t_work(ind+j)
3749 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3750 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3751 .and.(mnum.ne.5)) then
3754 dc(j,inres)=dc_work(ind+j)
3755 d_t(j,inres)=d_t_work(ind+j)
3761 end subroutine sd_verlet1_ciccotti
3762 !-----------------------------------------------------------------------------
3763 subroutine sd_verlet2_ciccotti
3765 ! Calculating the adjusted velocities for accelerations
3767 ! implicit real*8 (a-h,o-z)
3768 ! include 'DIMENSIONS'
3769 ! include 'COMMON.CONTROL'
3770 ! include 'COMMON.VAR'
3771 ! include 'COMMON.MD'
3773 ! include 'COMMON.LANGEVIN'
3775 ! include 'COMMON.LANGEVIN.lang0'
3777 ! include 'COMMON.CHAIN'
3778 ! include 'COMMON.DERIV'
3779 ! include 'COMMON.GEO'
3780 ! include 'COMMON.LOCAL'
3781 ! include 'COMMON.INTERACT'
3782 ! include 'COMMON.IOUNITS'
3783 ! include 'COMMON.NAMES'
3784 !el real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV !(MAXRES6) maxres6=6*maxres
3785 real(kind=8),dimension(6*nres) :: stochforcvecV !(MAXRES6) maxres6=6*maxres
3786 !el common /stochcalc/ stochforcvec
3787 real(kind=8) :: ddt1,ddt2
3788 integer :: i,j,ind,inres
3790 ! Compute the stochastic forces which contribute to velocity change
3792 call stochastic_force(stochforcvecV)
3799 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3800 ! ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
3801 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
3803 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
3807 d_t(j,0)=d_t_work(j)
3812 d_t(j,i)=d_t_work(ind+j)
3818 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3820 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3823 d_t(j,inres)=d_t_work(ind+j)
3829 end subroutine sd_verlet2_ciccotti
3831 !-----------------------------------------------------------------------------
3833 !-----------------------------------------------------------------------------
3834 subroutine inertia_tensor
3836 ! Calculating the intertia tensor for the entire protein in order to
3837 ! remove the perpendicular components of velocity matrix which cause
3838 ! the molecule to rotate.
3841 ! implicit real*8 (a-h,o-z)
3842 ! include 'DIMENSIONS'
3843 ! include 'COMMON.CONTROL'
3844 ! include 'COMMON.VAR'
3845 ! include 'COMMON.MD'
3846 ! include 'COMMON.CHAIN'
3847 ! include 'COMMON.DERIV'
3848 ! include 'COMMON.GEO'
3849 ! include 'COMMON.LOCAL'
3850 ! include 'COMMON.INTERACT'
3851 ! include 'COMMON.IOUNITS'
3852 ! include 'COMMON.NAMES'
3854 real(kind=8),dimension(3,3) :: Im,Imcp,eigvec,Id
3855 real(kind=8),dimension(3) :: pr,eigval,L,vp,vrot
3856 real(kind=8) :: M_SC,mag,mag2,M_PEP
3857 real(kind=8),dimension(3,0:nres) :: vpp !(3,0:MAXRES)
3858 real(kind=8),dimension(3) :: vs_p,pp,incr,v
3859 real(kind=8),dimension(3,3) :: pr1,pr2
3861 !el common /gucio/ cm
3862 integer :: iti,inres,i,j,k,mnum
3873 ! calculating the center of the mass of the protein
3877 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
3878 if (itype(i,mnum).eq.ntyp1_molec(mnum)) cycle
3879 M_PEP=M_PEP+mp(mnum)
3881 cm(j)=cm(j)+(c(j,i)+0.5d0*dc(j,i))*mp(mnum)
3890 if (mnum.eq.5) cycle
3891 iti=iabs(itype(i,mnum))
3892 M_SC=M_SC+msc(iabs(iti),mnum)
3895 cm(j)=cm(j)+msc(iabs(iti),mnum)*c(j,inres)
3899 cm(j)=cm(j)/(M_SC+M_PEP)
3904 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
3906 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
3908 Im(1,1)=Im(1,1)+mp(mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
3909 Im(1,2)=Im(1,2)-mp(mnum)*pr(1)*pr(2)
3910 Im(1,3)=Im(1,3)-mp(mnum)*pr(1)*pr(3)
3911 Im(2,3)=Im(2,3)-mp(mnum)*pr(2)*pr(3)
3912 Im(2,2)=Im(2,2)+mp(mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
3913 Im(3,3)=Im(3,3)+mp(mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
3918 iti=iabs(itype(i,mnum))
3919 if (mnum.eq.5) cycle
3922 pr(j)=c(j,inres)-cm(j)
3924 Im(1,1)=Im(1,1)+msc(iabs(iti),mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
3925 Im(1,2)=Im(1,2)-msc(iabs(iti),mnum)*pr(1)*pr(2)
3926 Im(1,3)=Im(1,3)-msc(iabs(iti),mnum)*pr(1)*pr(3)
3927 Im(2,3)=Im(2,3)-msc(iabs(iti),mnum)*pr(2)*pr(3)
3928 Im(2,2)=Im(2,2)+msc(iabs(iti),mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
3929 Im(3,3)=Im(3,3)+msc(iabs(iti),mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
3934 Im(1,1)=Im(1,1)+Ip(mnum)*(1-dc_norm(1,i)*dc_norm(1,i))* &
3935 vbld(i+1)*vbld(i+1)*0.25d0
3936 Im(1,2)=Im(1,2)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(2,i))* &
3937 vbld(i+1)*vbld(i+1)*0.25d0
3938 Im(1,3)=Im(1,3)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(3,i))* &
3939 vbld(i+1)*vbld(i+1)*0.25d0
3940 Im(2,3)=Im(2,3)+Ip(mnum)*(-dc_norm(2,i)*dc_norm(3,i))* &
3941 vbld(i+1)*vbld(i+1)*0.25d0
3942 Im(2,2)=Im(2,2)+Ip(mnum)*(1-dc_norm(2,i)*dc_norm(2,i))* &
3943 vbld(i+1)*vbld(i+1)*0.25d0
3944 Im(3,3)=Im(3,3)+Ip(mnum)*(1-dc_norm(3,i)*dc_norm(3,i))* &
3945 vbld(i+1)*vbld(i+1)*0.25d0
3951 ! if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
3952 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3953 .and.(mnum.ne.5)) then
3954 iti=iabs(itype(i,mnum))
3956 Im(1,1)=Im(1,1)+Isc(iti,mnum)*(1-dc_norm(1,inres)* &
3957 dc_norm(1,inres))*vbld(inres)*vbld(inres)
3958 Im(1,2)=Im(1,2)-Isc(iti,mnum)*(dc_norm(1,inres)* &
3959 dc_norm(2,inres))*vbld(inres)*vbld(inres)
3960 Im(1,3)=Im(1,3)-Isc(iti,mnum)*(dc_norm(1,inres)* &
3961 dc_norm(3,inres))*vbld(inres)*vbld(inres)
3962 Im(2,3)=Im(2,3)-Isc(iti,mnum)*(dc_norm(2,inres)* &
3963 dc_norm(3,inres))*vbld(inres)*vbld(inres)
3964 Im(2,2)=Im(2,2)+Isc(iti,mnum)*(1-dc_norm(2,inres)* &
3965 dc_norm(2,inres))*vbld(inres)*vbld(inres)
3966 Im(3,3)=Im(3,3)+Isc(iti,mnum)*(1-dc_norm(3,inres)* &
3967 dc_norm(3,inres))*vbld(inres)*vbld(inres)
3972 ! write(iout,*) "The angular momentum before adjustment:"
3973 ! write(iout,*) (L(j),j=1,3)
3979 ! Copying the Im matrix for the djacob subroutine
3987 ! Finding the eigenvectors and eignvalues of the inertia tensor
3988 call djacob(3,3,10000,1.0d-10,Imcp,eigvec,eigval)
3989 ! write (iout,*) "Eigenvalues & Eigenvectors"
3990 ! write (iout,'(5x,3f10.5)') (eigval(i),i=1,3)
3993 ! write (iout,'(i5,3f10.5)') i,(eigvec(i,j),j=1,3)
3995 ! Constructing the diagonalized matrix
3997 if (dabs(eigval(i)).gt.1.0d-15) then
3998 Id(i,i)=1.0d0/eigval(i)
4005 Imcp(i,j)=eigvec(j,i)
4011 pr1(i,j)=pr1(i,j)+Id(i,k)*Imcp(k,j)
4018 pr2(i,j)=pr2(i,j)+eigvec(i,k)*pr1(k,j)
4022 ! Calculating the total rotational velocity of the molecule
4025 vrot(i)=vrot(i)+pr2(i,j)*L(j)
4028 ! Resetting the velocities
4030 call vecpr(vrot(1),dc(1,i),vp)
4032 d_t(j,i)=d_t(j,i)-vp(j)
4037 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4038 .and.(mnum.ne.5)) then
4039 ! if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
4040 ! if(itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
4042 call vecpr(vrot(1),dc(1,inres),vp)
4044 d_t(j,inres)=d_t(j,inres)-vp(j)
4049 ! write(iout,*) "The angular momentum after adjustment:"
4050 ! write(iout,*) (L(j),j=1,3)
4053 end subroutine inertia_tensor
4054 !-----------------------------------------------------------------------------
4055 subroutine angmom(cm,L)
4058 ! implicit real*8 (a-h,o-z)
4059 ! include 'DIMENSIONS'
4060 ! include 'COMMON.CONTROL'
4061 ! include 'COMMON.VAR'
4062 ! include 'COMMON.MD'
4063 ! include 'COMMON.CHAIN'
4064 ! include 'COMMON.DERIV'
4065 ! include 'COMMON.GEO'
4066 ! include 'COMMON.LOCAL'
4067 ! include 'COMMON.INTERACT'
4068 ! include 'COMMON.IOUNITS'
4069 ! include 'COMMON.NAMES'
4070 real(kind=8) :: mscab
4071 real(kind=8),dimension(3) :: L,cm,pr,vp,vrot,incr,v,pp
4072 integer :: iti,inres,i,j,mnum
4073 ! Calculate the angular momentum
4082 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
4084 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
4087 v(j)=incr(j)+0.5d0*d_t(j,i)
4090 incr(j)=incr(j)+d_t(j,i)
4092 call vecpr(pr(1),v(1),vp)
4094 L(j)=L(j)+mp(mnum)*vp(j)
4098 pp(j)=0.5d0*d_t(j,i)
4100 call vecpr(pr(1),pp(1),vp)
4102 L(j)=L(j)+Ip(mnum)*vp(j)
4110 iti=iabs(itype(i,mnum))
4118 pr(j)=c(j,inres)-cm(j)
4120 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4121 .and.(mnum.ne.5)) then
4123 v(j)=incr(j)+d_t(j,inres)
4130 call vecpr(pr(1),v(1),vp)
4131 ! write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),&
4132 ! " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
4134 L(j)=L(j)+mscab*vp(j)
4136 ! write (iout,*) "L",(l(j),j=1,3)
4137 if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4138 .and.(mnum.ne.5)) then
4140 v(j)=incr(j)+d_t(j,inres)
4142 call vecpr(dc(1,inres),d_t(1,inres),vp)
4144 L(j)=L(j)+Isc(iti,mnum)*vp(j)
4148 incr(j)=incr(j)+d_t(j,i)
4152 end subroutine angmom
4153 !-----------------------------------------------------------------------------
4154 subroutine vcm_vel(vcm)
4157 ! implicit real*8 (a-h,o-z)
4158 ! include 'DIMENSIONS'
4159 ! include 'COMMON.VAR'
4160 ! include 'COMMON.MD'
4161 ! include 'COMMON.CHAIN'
4162 ! include 'COMMON.DERIV'
4163 ! include 'COMMON.GEO'
4164 ! include 'COMMON.LOCAL'
4165 ! include 'COMMON.INTERACT'
4166 ! include 'COMMON.IOUNITS'
4167 real(kind=8),dimension(3) :: vcm,vv
4168 real(kind=8) :: summas,amas
4178 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
4180 summas=summas+mp(mnum)
4182 vcm(j)=vcm(j)+mp(mnum)*(vv(j)+0.5d0*d_t(j,i))
4186 amas=msc(iabs(itype(i,mnum)),mnum)
4191 if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4192 .and.(mnum.ne.5)) then
4194 vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres))
4198 vcm(j)=vcm(j)+amas*vv(j)
4202 vv(j)=vv(j)+d_t(j,i)
4205 ! write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas
4207 vcm(j)=vcm(j)/summas
4210 end subroutine vcm_vel
4211 !-----------------------------------------------------------------------------
4213 !-----------------------------------------------------------------------------
4215 ! RATTLE algorithm for velocity Verlet - step 1, UNRES
4219 ! implicit real*8 (a-h,o-z)
4220 ! include 'DIMENSIONS'
4222 ! include 'COMMON.CONTROL'
4223 ! include 'COMMON.VAR'
4224 ! include 'COMMON.MD'
4226 ! include 'COMMON.LANGEVIN'
4228 ! include 'COMMON.LANGEVIN.lang0'
4230 ! include 'COMMON.CHAIN'
4231 ! include 'COMMON.DERIV'
4232 ! include 'COMMON.GEO'
4233 ! include 'COMMON.LOCAL'
4234 ! include 'COMMON.INTERACT'
4235 ! include 'COMMON.IOUNITS'
4236 ! include 'COMMON.NAMES'
4237 ! include 'COMMON.TIME1'
4238 !el real(kind=8) :: gginv(2*nres,2*nres),&
4239 !el gdc(3,2*nres,2*nres)
4240 real(kind=8) :: dC_uncor(3,2*nres) !,&
4241 !el real(kind=8) :: Cmat(2*nres,2*nres)
4242 real(kind=8) :: x(2*nres),xcorr(3,2*nres) !maxres2=2*maxres
4243 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
4244 !el common /przechowalnia/ nbond
4245 integer :: max_rattle = 5
4246 logical :: lprn = .false., lprn1 = .false., not_done
4247 real(kind=8) :: tol_rattle = 1.0d-5
4249 integer :: ii,i,j,jj,l,ind,ind1,nres2
4252 !el /common/ przechowalnia
4254 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4255 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4256 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4258 if (lprn) write (iout,*) "RATTLE1"
4262 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4263 .and.(mnum.ne.5)) nbond=nbond+1
4265 ! Make a folded form of the Ginv-matrix
4278 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4283 if (itype(k,1).ne.10 .and. itype(k,mnum).ne.ntyp1_molec(mnum)&
4284 .and.(mnum.ne.5)) then
4288 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4296 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4307 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4311 if (itype(k,1).ne.10) then
4315 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4323 write (iout,*) "Matrix GGinv"
4324 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
4330 if (iter.gt.max_rattle) then
4331 write (iout,*) "Error - too many iterations in RATTLE."
4334 ! Calculate the matrix C = GG**(-1) dC_old o dC
4339 dC_uncor(j,ind1)=dC(j,i)
4343 if (itype(i,1).ne.10) then
4346 dC_uncor(j,ind1)=dC(j,i+nres)
4355 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
4359 if (itype(k,1).ne.10) then
4362 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
4367 ! Calculate deviations from standard virtual-bond lengths
4371 x(ind)=vbld(i+1)**2-vbl**2
4374 if (itype(i,1).ne.10) then
4376 x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4380 write (iout,*) "Coordinates and violations"
4382 write(iout,'(i5,3f10.5,5x,e15.5)') &
4383 i,(dC_uncor(j,i),j=1,3),x(i)
4385 write (iout,*) "Velocities and violations"
4389 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4390 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4394 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4395 .and.(mnum.ne.5)) then
4398 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4399 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4400 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4403 ! write (iout,*) "gdc"
4405 ! write (iout,*) "i",i
4407 ! write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4413 if (dabs(x(i)).gt.xmax) then
4417 if (xmax.lt.tol_rattle) then
4421 ! Calculate the matrix of the system of equations
4426 Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
4431 write (iout,*) "Matrix Cmat"
4432 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4434 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4435 ! Add constraint term to positions
4442 xx = xx+x(ii)*gdc(j,ind,ii)
4446 d_t_new(j,i)=d_t_new(j,i)-xx/d_time
4450 if (itype(i,1).ne.10) then
4455 xx = xx+x(ii)*gdc(j,ind,ii)
4458 dC(j,i+nres)=dC(j,i+nres)-xx
4459 d_t_new(j,i+nres)=d_t_new(j,i+nres)-xx/d_time
4463 ! Rebuild the chain using the new coordinates
4464 call chainbuild_cart
4466 write (iout,*) "New coordinates, Lagrange multipliers,",&
4467 " and differences between actual and standard bond lengths"
4471 xx=vbld(i+1)**2-vbl**2
4472 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4473 i,(dC(j,i),j=1,3),x(ind),xx
4477 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4480 xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4481 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4482 i,(dC(j,i+nres),j=1,3),x(ind),xx
4485 write (iout,*) "Velocities and violations"
4489 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4490 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4493 if (itype(i,1).ne.10) then
4495 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4496 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4497 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4504 10 write (iout,*) "Error - singularity in solving the system",&
4505 " of equations for Lagrange multipliers."
4509 "RATTLE inactive; use -DRATTLE switch at compile time."
4512 end subroutine rattle1
4513 !-----------------------------------------------------------------------------
4515 ! RATTLE algorithm for velocity Verlet - step 2, UNRES
4519 ! implicit real*8 (a-h,o-z)
4520 ! include 'DIMENSIONS'
4522 ! include 'COMMON.CONTROL'
4523 ! include 'COMMON.VAR'
4524 ! include 'COMMON.MD'
4526 ! include 'COMMON.LANGEVIN'
4528 ! include 'COMMON.LANGEVIN.lang0'
4530 ! include 'COMMON.CHAIN'
4531 ! include 'COMMON.DERIV'
4532 ! include 'COMMON.GEO'
4533 ! include 'COMMON.LOCAL'
4534 ! include 'COMMON.INTERACT'
4535 ! include 'COMMON.IOUNITS'
4536 ! include 'COMMON.NAMES'
4537 ! include 'COMMON.TIME1'
4538 !el real(kind=8) :: gginv(2*nres,2*nres),&
4539 !el gdc(3,2*nres,2*nres)
4540 real(kind=8) :: dC_uncor(3,2*nres) !,&
4541 !el Cmat(2*nres,2*nres)
4542 real(kind=8) :: x(2*nres) !maxres2=2*maxres
4543 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
4544 !el common /przechowalnia/ nbond
4545 integer :: max_rattle = 5
4546 logical :: lprn = .false., lprn1 = .false., not_done
4547 real(kind=8) :: tol_rattle = 1.0d-5
4551 !el /common/ przechowalnia
4552 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4553 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4554 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4556 if (lprn) write (iout,*) "RATTLE2"
4557 if (lprn) write (iout,*) "Velocity correction"
4558 ! Calculate the matrix G dC
4564 gdc(j,i,ind)=GGinv(i,ind)*dC(j,k)
4569 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4570 .and.(mnum.ne.5)) then
4573 gdc(j,i,ind)=GGinv(i,ind)*dC(j,k+nres)
4579 ! write (iout,*) "gdc"
4581 ! write (iout,*) "i",i
4583 ! write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4587 ! Calculate the matrix of the system of equations
4594 Cmat(ind,j)=Cmat(ind,j)+dC(k,i)*gdc(k,ind,j)
4600 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4601 .and.(mnum.ne.5)) then
4606 Cmat(ind,j)=Cmat(ind,j)+dC(k,i+nres)*gdc(k,ind,j)
4611 ! Calculate the scalar product dC o d_t_new
4615 x(ind)=scalar(d_t(1,i),dC(1,i))
4619 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4620 .and.(mnum.ne.5)) then
4622 x(ind)=scalar(d_t(1,i+nres),dC(1,i+nres))
4626 write (iout,*) "Velocities and violations"
4630 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4631 i,ind,(d_t(j,i),j=1,3),x(ind)
4635 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4636 .and.(mnum.ne.5)) then
4638 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4639 i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind)
4645 if (dabs(x(i)).gt.xmax) then
4649 if (xmax.lt.tol_rattle) then
4654 write (iout,*) "Matrix Cmat"
4655 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4657 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4658 ! Add constraint term to velocities
4665 xx = xx+x(ii)*gdc(j,ind,ii)
4667 d_t(j,i)=d_t(j,i)-xx
4672 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4673 .and.(mnum.ne.5)) then
4678 xx = xx+x(ii)*gdc(j,ind,ii)
4680 d_t(j,i+nres)=d_t(j,i+nres)-xx
4686 "New velocities, Lagrange multipliers violations"
4690 if (lprn) write (iout,'(2i5,3f10.5,5x,2e15.5)') &
4691 i,ind,(d_t(j,i),j=1,3),x(ind),scalar(d_t(1,i),dC(1,i))
4695 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4698 write (iout,'(2i5,3f10.5,5x,2e15.5)') &
4699 i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind),&
4700 scalar(d_t(1,i+nres),dC(1,i+nres))
4706 10 write (iout,*) "Error - singularity in solving the system",&
4707 " of equations for Lagrange multipliers."
4711 "RATTLE inactive; use -DRATTLE option at compile time."
4714 end subroutine rattle2
4715 !-----------------------------------------------------------------------------
4716 subroutine rattle_brown
4717 ! RATTLE/LINCS algorithm for Brownian dynamics, UNRES
4721 ! implicit real*8 (a-h,o-z)
4722 ! include 'DIMENSIONS'
4724 ! include 'COMMON.CONTROL'
4725 ! include 'COMMON.VAR'
4726 ! include 'COMMON.MD'
4728 ! include 'COMMON.LANGEVIN'
4730 ! include 'COMMON.LANGEVIN.lang0'
4732 ! include 'COMMON.CHAIN'
4733 ! include 'COMMON.DERIV'
4734 ! include 'COMMON.GEO'
4735 ! include 'COMMON.LOCAL'
4736 ! include 'COMMON.INTERACT'
4737 ! include 'COMMON.IOUNITS'
4738 ! include 'COMMON.NAMES'
4739 ! include 'COMMON.TIME1'
4740 !el real(kind=8) :: gginv(2*nres,2*nres),&
4741 !el gdc(3,2*nres,2*nres)
4742 real(kind=8) :: dC_uncor(3,2*nres) !,&
4743 !el real(kind=8) :: Cmat(2*nres,2*nres)
4744 real(kind=8) :: x(2*nres) !maxres2=2*maxres
4745 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
4746 !el common /przechowalnia/ nbond
4747 integer :: max_rattle = 5
4748 logical :: lprn = .false., lprn1 = .false., not_done
4749 real(kind=8) :: tol_rattle = 1.0d-5
4753 !el /common/ przechowalnia
4754 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4755 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4756 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4759 if (lprn) write (iout,*) "RATTLE_BROWN"
4762 if (itype(i,1).ne.10) nbond=nbond+1
4764 ! Make a folded form of the Ginv-matrix
4777 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4781 if (itype(k,1).ne.10) then
4785 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4792 if (itype(i,1).ne.10) then
4802 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4806 if (itype(k,1).ne.10) then
4810 if (j.eq.1 .and. l.eq.1)GGinv(ii,jj)=fricmat(ind,ind1)
4818 write (iout,*) "Matrix GGinv"
4819 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
4825 if (iter.gt.max_rattle) then
4826 write (iout,*) "Error - too many iterations in RATTLE."
4829 ! Calculate the matrix C = GG**(-1) dC_old o dC
4834 dC_uncor(j,ind1)=dC(j,i)
4838 if (itype(i,1).ne.10) then
4841 dC_uncor(j,ind1)=dC(j,i+nres)
4850 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
4854 if (itype(k,1).ne.10) then
4857 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
4862 ! Calculate deviations from standard virtual-bond lengths
4866 x(ind)=vbld(i+1)**2-vbl**2
4869 if (itype(i,1).ne.10) then
4871 x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4875 write (iout,*) "Coordinates and violations"
4877 write(iout,'(i5,3f10.5,5x,e15.5)') &
4878 i,(dC_uncor(j,i),j=1,3),x(i)
4880 write (iout,*) "Velocities and violations"
4884 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4885 i,ind,(d_t(j,i),j=1,3),scalar(d_t(1,i),dC_old(1,i))
4888 if (itype(i,1).ne.10) then
4890 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4891 i+nres,ind,(d_t(j,i+nres),j=1,3),&
4892 scalar(d_t(1,i+nres),dC_old(1,i+nres))
4895 write (iout,*) "gdc"
4897 write (iout,*) "i",i
4899 write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4905 if (dabs(x(i)).gt.xmax) then
4909 if (xmax.lt.tol_rattle) then
4913 ! Calculate the matrix of the system of equations
4918 Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
4923 write (iout,*) "Matrix Cmat"
4924 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4926 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4927 ! Add constraint term to positions
4934 xx = xx+x(ii)*gdc(j,ind,ii)
4937 d_t(j,i)=d_t(j,i)+xx/d_time
4942 if (itype(i,1).ne.10) then
4947 xx = xx+x(ii)*gdc(j,ind,ii)
4950 d_t(j,i+nres)=d_t(j,i+nres)+xx/d_time
4951 dC(j,i+nres)=dC(j,i+nres)+xx
4955 ! Rebuild the chain using the new coordinates
4956 call chainbuild_cart
4958 write (iout,*) "New coordinates, Lagrange multipliers,",&
4959 " and differences between actual and standard bond lengths"
4963 xx=vbld(i+1)**2-vbl**2
4964 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4965 i,(dC(j,i),j=1,3),x(ind),xx
4968 if (itype(i,1).ne.10) then
4970 xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4971 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4972 i,(dC(j,i+nres),j=1,3),x(ind),xx
4975 write (iout,*) "Velocities and violations"
4979 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4980 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4983 if (itype(i,1).ne.10) then
4985 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4986 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4987 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4994 10 write (iout,*) "Error - singularity in solving the system",&
4995 " of equations for Lagrange multipliers."
4999 "RATTLE inactive; use -DRATTLE option at compile time"
5002 end subroutine rattle_brown
5003 !-----------------------------------------------------------------------------
5005 !-----------------------------------------------------------------------------
5006 subroutine friction_force
5011 ! implicit real*8 (a-h,o-z)
5012 ! include 'DIMENSIONS'
5013 ! include 'COMMON.VAR'
5014 ! include 'COMMON.CHAIN'
5015 ! include 'COMMON.DERIV'
5016 ! include 'COMMON.GEO'
5017 ! include 'COMMON.LOCAL'
5018 ! include 'COMMON.INTERACT'
5019 ! include 'COMMON.MD'
5021 ! include 'COMMON.LANGEVIN'
5023 ! include 'COMMON.LANGEVIN.lang0'
5025 ! include 'COMMON.IOUNITS'
5026 !el real(kind=8),dimension(6*nres) :: gamvec !(MAXRES6) maxres6=6*maxres
5027 !el common /syfek/ gamvec
5028 real(kind=8) :: vv(3),vvtot(3,nres),v_work(6*nres) !,&
5029 !el ginvfric(2*nres,2*nres) !maxres2=2*maxres
5030 !el common /przechowalnia/ ginvfric
5032 logical :: lprn = .false., checkmode = .false.
5033 integer :: i,j,ind,k,nres2,nres6,mnum
5037 if(.not.allocated(gamvec)) allocate(gamvec(nres6)) !(MAXRES6)
5038 if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
5046 d_t_work(j)=d_t(j,0)
5051 d_t_work(ind+j)=d_t(j,i)
5057 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5058 .and.(mnum.ne.5)) then
5060 d_t_work(ind+j)=d_t(j,i+nres)
5066 call fricmat_mult(d_t_work,fric_work)
5068 if (.not.checkmode) return
5071 write (iout,*) "d_t_work and fric_work"
5073 write (iout,'(i3,2e15.5)') i,d_t_work(i),fric_work(i)
5077 friction(j,0)=fric_work(j)
5082 friction(j,i)=fric_work(ind+j)
5088 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5089 .and.(mnum.ne.5)) then
5090 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
5092 friction(j,i+nres)=fric_work(ind+j)
5098 write(iout,*) "Friction backbone"
5100 write(iout,'(i5,3e15.5,5x,3e15.5)') &
5101 i,(friction(j,i),j=1,3),(d_t(j,i),j=1,3)
5103 write(iout,*) "Friction side chain"
5105 write(iout,'(i5,3e15.5,5x,3e15.5)') &
5106 i,(friction(j,i+nres),j=1,3),(d_t(j,i+nres),j=1,3)
5115 vvtot(j,i)=vv(j)+0.5d0*d_t(j,i)
5116 vvtot(j,i+nres)=vv(j)+d_t(j,i+nres)
5117 vv(j)=vv(j)+d_t(j,i)
5120 write (iout,*) "vvtot backbone and sidechain"
5122 write (iout,'(i5,3e15.5,5x,3e15.5)') i,(vvtot(j,i),j=1,3),&
5123 (vvtot(j,i+nres),j=1,3)
5128 v_work(ind+j)=vvtot(j,i)
5134 v_work(ind+j)=vvtot(j,i+nres)
5138 write (iout,*) "v_work gamvec and site-based friction forces"
5140 write (iout,'(i5,3e15.5)') i,v_work(i),gamvec(i),&
5144 ! fric_work1(i)=0.0d0
5146 ! fric_work1(i)=fric_work1(i)-A(j,i)*gamvec(j)*v_work(j)
5149 ! write (iout,*) "fric_work and fric_work1"
5151 ! write (iout,'(i5,2e15.5)') i,fric_work(i),fric_work1(i)
5157 ginvfric(i,j)=ginvfric(i,j)+ginv(i,k)*fricmat(k,j)
5161 write (iout,*) "ginvfric"
5163 write (iout,'(i5,100f8.3)') i,(ginvfric(i,j),j=1,dimen)
5165 write (iout,*) "symmetry check"
5168 write (iout,*) i,j,ginvfric(i,j)-ginvfric(j,i)
5173 end subroutine friction_force
5174 !-----------------------------------------------------------------------------
5175 subroutine setup_fricmat
5179 use control_data, only:time_Bcast
5180 use control, only:tcpu
5182 ! implicit real*8 (a-h,o-z)
5186 real(kind=8) :: time00
5188 ! include 'DIMENSIONS'
5189 ! include 'COMMON.VAR'
5190 ! include 'COMMON.CHAIN'
5191 ! include 'COMMON.DERIV'
5192 ! include 'COMMON.GEO'
5193 ! include 'COMMON.LOCAL'
5194 ! include 'COMMON.INTERACT'
5195 ! include 'COMMON.MD'
5196 ! include 'COMMON.SETUP'
5197 ! include 'COMMON.TIME1'
5198 ! integer licznik /0/
5201 ! include 'COMMON.LANGEVIN'
5203 ! include 'COMMON.LANGEVIN.lang0'
5205 ! include 'COMMON.IOUNITS'
5207 integer :: i,j,ind,ind1,m
5208 logical :: lprn = .false.
5209 real(kind=8) :: dtdi !el ,gamvec(2*nres)
5210 !el real(kind=8),dimension(2*nres,2*nres) :: ginvfric,fcopy
5211 ! real(kind=8),allocatable,dimension(:,:) :: fcopy
5212 !el real(kind=8),dimension(2*nres*(2*nres+1)/2) :: Ghalf !(mmaxres2) (mmaxres2=(maxres2*(maxres2+1)/2))
5213 !el common /syfek/ gamvec
5214 real(kind=8) :: work(8*2*nres)
5215 integer :: iwork(2*nres)
5216 !el common /przechowalnia/ ginvfric,Ghalf,fcopy
5217 integer :: ii,iti,k,l,nzero,nres2,nres6,ierr,mnum
5221 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
5222 if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
5223 if (fg_rank.ne.king) goto 10
5228 if(.not.allocated(gamvec)) allocate(gamvec(nres2)) !(MAXRES2)
5229 if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
5230 if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
5231 !el allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
5232 if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !maxres2=2*maxres
5234 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
5235 ! Zeroing out fricmat
5241 ! Load the friction coefficients corresponding to peptide groups
5246 gamvec(ind1)=gamp(mnum)
5248 ! Load the friction coefficients corresponding to side chains
5252 gamsc(ntyp1_molec(j),j)=1.0d0
5259 gamvec(ii)=gamsc(iabs(iti),mnum)
5261 if (surfarea) call sdarea(gamvec)
5263 ! write (iout,*) "Matrix A and vector gamma"
5265 ! write (iout,'(i2,$)') i
5267 ! write (iout,'(f4.1,$)') A(i,j)
5269 ! write (iout,'(f8.3)') gamvec(i)
5273 write (iout,*) "Vector gamvec"
5275 write (iout,'(i5,f10.5)') i, gamvec(i)
5279 ! The friction matrix
5284 dtdi=dtdi+A(j,k)*A(j,i)*gamvec(j)
5291 write (iout,'(//a)') "Matrix fricmat"
5292 call matout2(dimen,dimen,nres2,nres2,fricmat)
5294 if (lang.eq.2 .or. lang.eq.3) then
5295 ! Mass-scale the friction matrix if non-direct integration will be performed
5301 Ginvfric(i,j)=Ginvfric(i,j)+ &
5302 Gsqrm(i,k)*Gsqrm(l,j)*fricmat(k,l)
5307 ! Diagonalize the friction matrix
5312 Ghalf(ind)=Ginvfric(i,j)
5315 call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
5318 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
5319 " mass-scaled friction matrix"
5320 call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
5322 ! Precompute matrices for tinker stochastic integrator
5329 mt1(i,j)=mt1(i,j)+fricvec(k,i)*gsqrm(k,j)
5330 mt2(i,j)=mt2(i,j)+fricvec(k,i)*gsqrp(k,j)
5336 else if (lang.eq.4) then
5337 ! Diagonalize the friction matrix
5342 Ghalf(ind)=fricmat(i,j)
5345 call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
5348 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
5350 call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
5352 ! Determine the number of zero eigenvalues of the friction matrix
5353 nzero=max0(dimen-dimen1,0)
5354 ! do while (fricgam(nzero+1).le.1.0d-5 .and. nzero.lt.dimen)
5357 write (iout,*) "Number of zero eigenvalues:",nzero
5362 fricmat(i,j)=fricmat(i,j) &
5363 +fricvec(i,k)*fricvec(j,k)/fricgam(k)
5368 write (iout,'(//a)') "Generalized inverse of fricmat"
5369 call matout(dimen,dimen,nres6,nres6,fricmat)
5374 if (nfgtasks.gt.1) then
5375 if (fg_rank.eq.0) then
5376 ! The matching BROADCAST for fg processors is called in ERGASTULUM
5382 call MPI_Bcast(10,1,MPI_INTEGER,king,FG_COMM,IERROR)
5384 time_Bcast=time_Bcast+MPI_Wtime()-time00
5386 time_Bcast=time_Bcast+tcpu()-time00
5388 ! print *,"Processor",myrank,
5389 ! & " BROADCAST iorder in SETUP_FRICMAT"
5392 write (iout,*) "setup_fricmat licznik"!,licznik !sp
5398 ! Scatter the friction matrix
5399 call MPI_Scatterv(fricmat(1,1),nginv_counts(0),&
5400 nginv_start(0),MPI_DOUBLE_PRECISION,fcopy(1,1),&
5401 myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
5404 time_scatter=time_scatter+MPI_Wtime()-time00
5405 time_scatter_fmat=time_scatter_fmat+MPI_Wtime()-time00
5407 time_scatter=time_scatter+tcpu()-time00
5408 time_scatter_fmat=time_scatter_fmat+tcpu()-time00
5412 do j=1,2*my_ng_count
5413 fricmat(j,i)=fcopy(i,j)
5416 ! write (iout,*) "My chunk of fricmat"
5417 ! call MATOUT2(my_ng_count,dimen,maxres2,maxres2,fcopy)
5421 end subroutine setup_fricmat
5422 !-----------------------------------------------------------------------------
5423 subroutine stochastic_force(stochforcvec)
5426 use random, only:anorm_distr
5427 ! implicit real*8 (a-h,o-z)
5428 ! include 'DIMENSIONS'
5429 use control, only: tcpu
5434 ! include 'COMMON.VAR'
5435 ! include 'COMMON.CHAIN'
5436 ! include 'COMMON.DERIV'
5437 ! include 'COMMON.GEO'
5438 ! include 'COMMON.LOCAL'
5439 ! include 'COMMON.INTERACT'
5440 ! include 'COMMON.MD'
5441 ! include 'COMMON.TIME1'
5443 ! include 'COMMON.LANGEVIN'
5445 ! include 'COMMON.LANGEVIN.lang0'
5447 ! include 'COMMON.IOUNITS'
5449 real(kind=8) :: x,sig,lowb,highb
5450 real(kind=8) :: ff(3),force(3,0:2*nres),zeta2,lowb2
5451 real(kind=8) :: highb2,sig2,forcvec(6*nres),stochforcvec(6*nres)
5452 real(kind=8) :: time00
5453 logical :: lprn = .false.
5454 integer :: i,j,ind,mnum
5458 stochforc(j,i)=0.0d0
5468 ! Compute the stochastic forces acting on bodies. Store in force.
5474 force(j,i)=anorm_distr(x,sig,lowb,highb)
5482 force(j,i+nres)=anorm_distr(x,sig2,lowb2,highb2)
5486 time_fsample=time_fsample+MPI_Wtime()-time00
5488 time_fsample=time_fsample+tcpu()-time00
5490 ! Compute the stochastic forces acting on virtual-bond vectors.
5496 stochforc(j,i)=ff(j)+0.5d0*force(j,i)
5499 ff(j)=ff(j)+force(j,i)
5501 ! if (itype(i+1,1).ne.ntyp1) then
5503 if (itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
5505 stochforc(j,i)=stochforc(j,i)+force(j,i+nres+1)
5506 ff(j)=ff(j)+force(j,i+nres+1)
5511 stochforc(j,0)=ff(j)+force(j,nnt+nres)
5515 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5516 .and.(mnum.ne.5)) then
5517 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
5519 stochforc(j,i+nres)=force(j,i+nres)
5525 stochforcvec(j)=stochforc(j,0)
5530 stochforcvec(ind+j)=stochforc(j,i)
5536 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5537 .and.(mnum.ne.5)) then
5538 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
5540 stochforcvec(ind+j)=stochforc(j,i+nres)
5546 write (iout,*) "stochforcvec"
5548 write(iout,'(i5,e15.5)') i,stochforcvec(i)
5550 write(iout,*) "Stochastic forces backbone"
5552 write(iout,'(i5,3e15.5)') i,(stochforc(j,i),j=1,3)
5554 write(iout,*) "Stochastic forces side chain"
5556 write(iout,'(i5,3e15.5)') &
5557 i,(stochforc(j,i+nres),j=1,3)
5565 write (iout,*) i,ind
5567 forcvec(ind+j)=force(j,i)
5572 write (iout,*) i,ind
5574 forcvec(j+ind)=force(j,i+nres)
5579 write (iout,*) "forcvec"
5583 write (iout,'(2i3,2f10.5)') i,j,force(j,i),&
5590 write (iout,'(2i3,2f10.5)') i,j,force(j,i+nres),&
5599 end subroutine stochastic_force
5600 !-----------------------------------------------------------------------------
5601 subroutine sdarea(gamvec)
5603 ! Scale the friction coefficients according to solvent accessible surface areas
5604 ! Code adapted from TINKER
5608 ! implicit real*8 (a-h,o-z)
5609 ! include 'DIMENSIONS'
5610 ! include 'COMMON.CONTROL'
5611 ! include 'COMMON.VAR'
5612 ! include 'COMMON.MD'
5614 ! include 'COMMON.LANGEVIN'
5616 ! include 'COMMON.LANGEVIN.lang0'
5618 ! include 'COMMON.CHAIN'
5619 ! include 'COMMON.DERIV'
5620 ! include 'COMMON.GEO'
5621 ! include 'COMMON.LOCAL'
5622 ! include 'COMMON.INTERACT'
5623 ! include 'COMMON.IOUNITS'
5624 ! include 'COMMON.NAMES'
5625 real(kind=8),dimension(2*nres) :: radius,gamvec !(maxres2)
5626 real(kind=8),parameter :: twosix = 1.122462048309372981d0
5627 logical :: lprn = .false.
5628 real(kind=8) :: probe,area,ratio
5629 integer :: i,j,ind,iti,mnum
5631 ! determine new friction coefficients every few SD steps
5633 ! set the atomic radii to estimates of sigma values
5635 ! print *,"Entered sdarea"
5641 ! Load peptide group radii
5644 radius(i)=pstok(mnum)
5646 ! Load side chain radii
5650 radius(i+nres)=restok(iti,mnum)
5653 ! write (iout,*) "i",i," radius",radius(i)
5656 radius(i) = radius(i) / twosix
5657 if (radius(i) .ne. 0.0d0) radius(i) = radius(i) + probe
5660 ! scale atomic friction coefficients by accessible area
5662 if (lprn) write (iout,*) &
5663 "Original gammas, surface areas, scaling factors, new gammas, ",&
5664 "std's of stochastic forces"
5667 if (radius(i).gt.0.0d0) then
5668 call surfatom (i,area,radius)
5669 ratio = dmax1(area/(4.0d0*pi*radius(i)**2),1.0d-1)
5670 if (lprn) write (iout,'(i5,3f10.5,$)') &
5671 i,gamvec(ind+1),area,ratio
5674 gamvec(ind) = ratio * gamvec(ind)
5677 stdforcp(i)=stdfp(mnum)*dsqrt(gamvec(ind))
5678 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcp(i)
5682 if (radius(i+nres).gt.0.0d0) then
5683 call surfatom (i+nres,area,radius)
5684 ratio = dmax1(area/(4.0d0*pi*radius(i+nres)**2),1.0d-1)
5685 if (lprn) write (iout,'(i5,3f10.5,$)') &
5686 i,gamvec(ind+1),area,ratio
5689 gamvec(ind) = ratio * gamvec(ind)
5692 stdforcsc(i)=stdfsc(itype(i,mnum),mnum)*dsqrt(gamvec(ind))
5693 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcsc(i)
5698 end subroutine sdarea
5699 !-----------------------------------------------------------------------------
5701 !-----------------------------------------------------------------------------
5704 ! ###################################################
5705 ! ## COPYRIGHT (C) 1996 by Jay William Ponder ##
5706 ! ## All Rights Reserved ##
5707 ! ###################################################
5709 ! ################################################################
5711 ! ## subroutine surfatom -- exposed surface area of an atom ##
5713 ! ################################################################
5716 ! "surfatom" performs an analytical computation of the surface
5717 ! area of a specified atom; a simplified version of "surface"
5719 ! literature references:
5721 ! T. J. Richmond, "Solvent Accessible Surface Area and
5722 ! Excluded Volume in Proteins", Journal of Molecular Biology,
5725 ! L. Wesson and D. Eisenberg, "Atomic Solvation Parameters
5726 ! Applied to Molecular Dynamics of Proteins in Solution",
5727 ! Protein Science, 1, 227-235 (1992)
5729 ! variables and parameters:
5731 ! ir number of atom for which area is desired
5732 ! area accessible surface area of the atom
5733 ! radius radii of each of the individual atoms
5736 subroutine surfatom(ir,area,radius)
5738 ! implicit real*8 (a-h,o-z)
5739 ! include 'DIMENSIONS'
5741 ! include 'COMMON.GEO'
5742 ! include 'COMMON.IOUNITS'
5744 integer :: nsup,nstart_sup
5745 ! double precision c,dc,dc_old,d_c_work,xloc,xrot,dc_norm
5746 ! common /chain/ c(3,maxres2+2),dc(3,0:maxres2),dc_old(3,0:maxres2),
5747 ! & xloc(3,maxres),xrot(3,maxres),dc_norm(3,0:maxres2),
5748 ! & dc_work(MAXRES6),nres,nres0
5749 integer,parameter :: maxarc=300
5753 integer :: mi,ni,narc
5754 integer :: key(maxarc)
5755 integer :: intag(maxarc)
5756 integer :: intag1(maxarc)
5757 real(kind=8) :: area,arcsum
5758 real(kind=8) :: arclen,exang
5759 real(kind=8) :: delta,delta2
5760 real(kind=8) :: eps,rmove
5761 real(kind=8) :: xr,yr,zr
5762 real(kind=8) :: rr,rrsq
5763 real(kind=8) :: rplus,rminus
5764 real(kind=8) :: axx,axy,axz
5765 real(kind=8) :: ayx,ayy
5766 real(kind=8) :: azx,azy,azz
5767 real(kind=8) :: uxj,uyj,uzj
5768 real(kind=8) :: tx,ty,tz
5769 real(kind=8) :: txb,tyb,td
5770 real(kind=8) :: tr2,tr,txr,tyr
5771 real(kind=8) :: tk1,tk2
5772 real(kind=8) :: thec,the,t,tb
5773 real(kind=8) :: txk,tyk,tzk
5774 real(kind=8) :: t1,ti,tf,tt
5775 real(kind=8) :: txj,tyj,tzj
5776 real(kind=8) :: ccsq,cc,xysq
5777 real(kind=8) :: bsqk,bk,cosine
5778 real(kind=8) :: dsqj,gi,pix2
5779 real(kind=8) :: therk,dk,gk
5780 real(kind=8) :: risqk,rik
5781 real(kind=8) :: radius(2*nres) !(maxatm) (maxatm=maxres2)
5782 real(kind=8) :: ri(maxarc),risq(maxarc)
5783 real(kind=8) :: ux(maxarc),uy(maxarc),uz(maxarc)
5784 real(kind=8) :: xc(maxarc),yc(maxarc),zc(maxarc)
5785 real(kind=8) :: xc1(maxarc),yc1(maxarc),zc1(maxarc)
5786 real(kind=8) :: dsq(maxarc),bsq(maxarc)
5787 real(kind=8) :: dsq1(maxarc),bsq1(maxarc)
5788 real(kind=8) :: arci(maxarc),arcf(maxarc)
5789 real(kind=8) :: ex(maxarc),lt(maxarc),gr(maxarc)
5790 real(kind=8) :: b(maxarc),b1(maxarc),bg(maxarc)
5791 real(kind=8) :: kent(maxarc),kout(maxarc)
5792 real(kind=8) :: ther(maxarc)
5793 logical :: moved,top
5794 logical :: omit(maxarc)
5797 ! maxatm = 2*nres !maxres2 maxres2=2*maxres
5798 ! maxlight = 8*maxatm
5801 ! maxtors = 4*maxatm
5803 ! zero out the surface area for the sphere of interest
5806 ! write (2,*) "ir",ir," radius",radius(ir)
5807 if (radius(ir) .eq. 0.0d0) return
5809 ! set the overlap significance and connectivity shift
5813 delta2 = delta * delta
5818 ! store coordinates and radius of the sphere of interest
5826 ! initialize values of some counters and summations
5835 ! test each sphere to see if it overlaps the sphere of interest
5838 if (i.eq.ir .or. radius(i).eq.0.0d0) goto 30
5839 rplus = rr + radius(i)
5841 if (abs(tx) .ge. rplus) goto 30
5843 if (abs(ty) .ge. rplus) goto 30
5845 if (abs(tz) .ge. rplus) goto 30
5847 ! check for sphere overlap by testing distance against radii
5849 xysq = tx*tx + ty*ty
5850 if (xysq .lt. delta2) then
5857 if (rplus-cc .le. delta) goto 30
5858 rminus = rr - radius(i)
5860 ! check to see if sphere of interest is completely buried
5862 if (cc-abs(rminus) .le. delta) then
5863 if (rminus .le. 0.0d0) goto 170
5867 ! check for too many overlaps with sphere of interest
5869 if (io .ge. maxarc) then
5871 20 format (/,' SURFATOM -- Increase the Value of MAXARC')
5875 ! get overlap between current sphere and sphere of interest
5884 gr(io) = (ccsq+rplus*rminus) / (2.0d0*rr*b1(io))
5890 ! case where no other spheres overlap the sphere of interest
5893 area = 4.0d0 * pi * rrsq
5897 ! case where only one sphere overlaps the sphere of interest
5900 area = pix2 * (1.0d0 + gr(1))
5901 area = mod(area,4.0d0*pi) * rrsq
5905 ! case where many spheres intersect the sphere of interest;
5906 ! sort the intersecting spheres by their degree of overlap
5908 call sort2 (io,gr,key)
5911 intag(i) = intag1(k)
5920 ! get radius of each overlap circle on surface of the sphere
5925 risq(i) = rrsq - gi*gi
5926 ri(i) = sqrt(risq(i))
5927 ther(i) = 0.5d0*pi - asin(min(1.0d0,max(-1.0d0,gr(i))))
5930 ! find boundary of inaccessible area on sphere of interest
5933 if (.not. omit(k)) then
5940 ! check to see if J circle is intersecting K circle;
5941 ! get distance between circle centers and sum of radii
5944 if (omit(j)) goto 60
5945 cc = (txk*xc(j)+tyk*yc(j)+tzk*zc(j))/(bk*b(j))
5946 cc = acos(min(1.0d0,max(-1.0d0,cc)))
5947 td = therk + ther(j)
5949 ! check to see if circles enclose separate regions
5951 if (cc .ge. td) goto 60
5953 ! check for circle J completely inside circle K
5955 if (cc+ther(j) .lt. therk) goto 40
5957 ! check for circles that are essentially parallel
5959 if (cc .gt. delta) goto 50
5964 ! check to see if sphere of interest is completely buried
5967 if (pix2-cc .le. td) goto 170
5973 ! find T value of circle intersections
5976 if (omit(k)) goto 110
5991 ! rotation matrix elements
6003 if (.not. omit(j)) then
6008 ! rotate spheres so K vector colinear with z-axis
6010 uxj = txj*axx + tyj*axy - tzj*axz
6011 uyj = tyj*ayy - txj*ayx
6012 uzj = txj*azx + tyj*azy + tzj*azz
6013 cosine = min(1.0d0,max(-1.0d0,uzj/b(j)))
6014 if (acos(cosine) .lt. therk+ther(j)) then
6015 dsqj = uxj*uxj + uyj*uyj
6020 tr2 = risqk*dsqj - tb*tb
6026 ! get T values of intersection for K circle
6029 tb = min(1.0d0,max(-1.0d0,tb))
6031 if (tyb-txr .lt. 0.0d0) tk1 = pix2 - tk1
6033 tb = min(1.0d0,max(-1.0d0,tb))
6035 if (tyb+txr .lt. 0.0d0) tk2 = pix2 - tk2
6036 thec = (rrsq*uzj-gk*bg(j)) / (rik*ri(j)*b(j))
6037 if (abs(thec) .lt. 1.0d0) then
6039 else if (thec .ge. 1.0d0) then
6041 else if (thec .le. -1.0d0) then
6045 ! see if "tk1" is entry or exit point; check t=0 point;
6046 ! "ti" is exit point, "tf" is entry point
6048 cosine = min(1.0d0,max(-1.0d0, &
6049 (uzj*gk-uxj*rik)/(b(j)*rr)))
6050 if ((acos(cosine)-ther(j))*(tk2-tk1) .le. 0.0d0) then
6058 if (narc .ge. maxarc) then
6060 70 format (/,' SURFATOM -- Increase the Value',&
6064 if (tf .le. ti) then
6085 ! special case; K circle without intersections
6087 if (narc .le. 0) goto 90
6089 ! general case; sum up arclength and set connectivity code
6091 call sort2 (narc,arci,key)
6096 if (narc .gt. 1) then
6099 if (t .lt. arci(j)) then
6100 arcsum = arcsum + arci(j) - t
6101 exang = exang + ex(ni)
6103 if (jb .ge. maxarc) then
6105 80 format (/,' SURFATOM -- Increase the Value',&
6110 kent(jb) = maxarc*i + k
6112 kout(jb) = maxarc*k + i
6121 arcsum = arcsum + pix2 - t
6123 exang = exang + ex(ni)
6126 kent(jb) = maxarc*i + k
6128 kout(jb) = maxarc*k + i
6135 arclen = arclen + gr(k)*arcsum
6138 if (arclen .eq. 0.0d0) goto 170
6139 if (jb .eq. 0) goto 150
6141 ! find number of independent boundaries and check connectivity
6145 if (kout(k) .ne. 0) then
6152 if (m .eq. kent(ii)) then
6155 if (j .eq. jb) goto 150
6167 ! attempt to fix connectivity error by moving atom slightly
6171 140 format (/,' SURFATOM -- Connectivity Error at Atom',i6)
6180 ! compute the exposed surface area for the sphere of interest
6183 area = ib*pix2 + exang + arclen
6184 area = mod(area,4.0d0*pi) * rrsq
6186 ! attempt to fix negative area by moving atom slightly
6188 if (area .lt. 0.0d0) then
6191 160 format (/,' SURFATOM -- Negative Area at Atom',i6)
6202 end subroutine surfatom
6203 !----------------------------------------------------------------
6204 !----------------------------------------------------------------
6205 subroutine alloc_MD_arrays
6206 !EL Allocation of arrays used by MD module
6208 integer :: nres2,nres6
6211 !----------------------
6215 allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
6216 allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
6217 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
6218 allocate(fricvec(nres2,nres2))
6219 allocate(pfric_mat(nres2,nres2),vfric_mat(nres2,nres2))
6220 allocate(afric_mat(nres2,nres2),prand_mat(nres2,nres2))
6221 allocate(vrand_mat1(nres2,nres2),vrand_mat2(nres2,nres2)) !(MAXRES2,MAXRES2)
6222 allocate(pfric0_mat(nres2,nres2,0:maxflag_stoch))
6223 allocate(afric0_mat(nres2,nres2,0:maxflag_stoch))
6224 allocate(vfric0_mat(nres2,nres2,0:maxflag_stoch))
6225 allocate(prand0_mat(nres2,nres2,0:maxflag_stoch))
6226 allocate(vrand0_mat1(nres2,nres2,0:maxflag_stoch))
6227 allocate(vrand0_mat2(nres2,nres2,0:maxflag_stoch)) !(MAXRES2,MAXRES2,0:maxflag_stoch)
6228 allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
6230 allocate(mt1(nres2,nres2),mt2(nres2,nres2),mt3(nres2,nres2)) !(maxres2,maxres2)
6231 !----------------------
6233 ! commom.langevin.lang0
6235 allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
6236 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
6237 allocate(fricvec(nres2,nres2)) !(MAXRES2,MAXRES2)
6238 allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
6239 allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
6242 if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2))
6243 !----------------------
6244 ! commom.hairpin in CSA module
6245 !----------------------
6246 ! common.mce in MCM_MD module
6247 !----------------------
6249 ! common /mdgrad/ in module.energy
6250 ! common /back_constr/ in module.energy
6251 ! common /qmeas/ in module.energy
6254 allocate(potEcomp(0:n_ene+4)) !(0:n_ene+4)
6256 allocate(d_t(3,0:nres2),d_a(3,0:nres2),d_t_old(3,0:nres2)) !(3,0:MAXRES2)
6257 allocate(d_a_work(nres6)) !(6*MAXRES)
6259 allocate(DM(nres2),DU1(nres2),DU2(nres2))
6260 allocate(DMorig(nres2),DU1orig(nres2),DU2orig(nres2))
6262 write (iout,*) "Before A Allocation",nfgtasks-1
6264 allocate(Gmat(nres2,nres2),A(nres2,nres2))
6265 if(.not.allocated(Ginv)) allocate(Ginv(nres2,nres2)) !in control: ergastulum
6266 allocate(Gsqrp(nres2,nres2),Gsqrm(nres2,nres2),Gvec(nres2,nres2)) !(maxres2,maxres2)
6268 allocate(Geigen(nres2)) !(maxres2)
6269 if(.not.allocated(vtot)) allocate(vtot(nres2)) !(maxres2)
6270 ! common /inertia/ in io_conf: parmread
6271 ! real(kind=8),dimension(:),allocatable :: ISC,msc !(ntyp+1)
6272 ! common /langevin/in io read_MDpar
6273 ! real(kind=8),dimension(:),allocatable :: gamsc !(ntyp1)
6274 ! real(kind=8),dimension(:),allocatable :: stdfsc !(ntyp)
6275 ! in io_conf: parmread
6276 ! real(kind=8),dimension(:),allocatable :: restok !(ntyp+1)
6277 ! common /mdpmpi/ in control: ergastulum
6278 if(.not.allocated(ng_start)) allocate(ng_start(0:nfgtasks-1))
6279 if(.not.allocated(ng_counts)) allocate(ng_counts(0:nfgtasks-1))
6280 if(.not.allocated(nginv_counts)) allocate(nginv_counts(0:nfgtasks-1)) !(0:MaxProcs-1)
6281 if(.not.allocated(nginv_start)) allocate(nginv_start(0:nfgtasks)) !(0:MaxProcs)
6282 !----------------------
6283 ! common.muca in read_muca
6284 ! common /double_muca/
6285 ! real(kind=8) :: elow,ehigh,factor,hbin,factor_min
6286 ! real(kind=8),dimension(:),allocatable :: emuca,nemuca,&
6287 ! nemuca2,hist !(4*maxres)
6288 ! common /integer_muca/
6289 ! integer :: nmuca,imtime,muca_smooth
6291 ! real(kind=8),dimension(:),allocatable :: elowi,ehighi !(maxprocs)
6292 !----------------------
6294 ! common /mdgrad/ in module.energy
6295 ! common /back_constr/ in module.energy
6296 ! common /qmeas/ in module.energy
6300 allocate(d_t_work(nres6),d_t_work_new(nres6),d_af_work(nres6))
6301 allocate(d_as_work(nres6),kinetic_force(nres6)) !(MAXRES6)
6302 allocate(d_t_new(3,0:nres2),d_a_old(3,0:nres2),d_a_short(3,0:nres2)) !,d_a !(3,0:MAXRES2)
6303 allocate(stdforcp(nres),stdforcsc(nres)) !(MAXRES)
6304 !----------------------
6306 allocate(D_ban(nres6)) !(MAXRES6) maxres6=6*maxres
6307 ! common /stochcalc/ stochforcvec
6308 allocate(stochforcvec(nres6)) !(MAXRES6) maxres6=6*maxres
6309 !----------------------
6311 end subroutine alloc_MD_arrays
6312 !-----------------------------------------------------------------------------
6313 !-----------------------------------------------------------------------------