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
817 integer,dimension(4,nres/3) :: iharp !(4,nres/3)(4,maxres/3)
819 real(kind=8) :: tt0,scalfac
820 integer :: nres2,itime
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)
991 if (mod(itime,ntwe).eq.0) then
994 ! call check_ecartint
1004 v_work(ind)=d_t(j,i)
1009 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1.and.mnum.ne.5) then
1012 v_work(ind)=d_t(j,i+nres)
1017 write (66,'(80f10.5)') &
1018 ((d_t(j,i),j=1,3),i=0,nres-1),((d_t(j,i+nres),j=1,3),i=1,nres)
1022 v_transf(i)=v_transf(i)+gvec(j,i)*v_work(j)
1024 v_transf(i)= v_transf(i)*dsqrt(geigen(i))
1026 write (67,'(80f10.5)') (v_transf(i),i=1,ind)
1029 if (mod(itime,ntwx).eq.0) then
1031 write (tytul,'("time",f8.2)') totT
1033 call hairpin(.true.,nharp,iharp)
1034 call secondary2(.true.)
1035 call pdbout(potE,tytul,ipdb)
1040 if (rstcount.eq.1000.or.itime.eq.n_timestep) then
1041 open(irest2,file=rest2name,status='unknown')
1042 write(irest2,*) totT,EK,potE,totE,t_bath
1044 ! AL 4/17/17: Now writing d_t(0,:) too
1046 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
1048 ! AL 4/17/17: Now writing d_c(0,:) too
1050 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
1058 t_MD=MPI_Wtime()-tt0
1062 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))') &
1064 'MD calculations setup:',t_MDsetup,&
1065 'Energy & gradient evaluation:',t_enegrad,&
1066 'Stochastic MD setup:',t_langsetup,&
1067 'Stochastic MD step setup:',t_sdsetup,&
1069 write (iout,'(/28(1h=),a25,27(1h=))') &
1070 ' End of MD calculation '
1072 write (iout,*) "time for etotal",t_etotal," elong",t_elong,&
1074 write (iout,*) "time_fric",time_fric," time_stoch",time_stoch,&
1075 " time_fricmatmult",time_fricmatmult," time_fsample ",&
1080 !-----------------------------------------------------------------------------
1081 subroutine velverlet_step(itime)
1082 !-------------------------------------------------------------------------------
1083 ! Perform a single velocity Verlet step; the time step can be rescaled if
1084 ! increments in accelerations exceed the threshold
1085 !-------------------------------------------------------------------------------
1086 ! implicit real*8 (a-h,o-z)
1087 ! include 'DIMENSIONS'
1089 use control, only:tcpu
1091 use minimm, only:minim_dc
1094 integer :: ierror,ierrcode
1095 real(kind=8) :: errcode
1097 ! include 'COMMON.SETUP'
1098 ! include 'COMMON.VAR'
1099 ! include 'COMMON.MD'
1101 ! include 'COMMON.LANGEVIN'
1103 ! include 'COMMON.LANGEVIN.lang0'
1105 ! include 'COMMON.CHAIN'
1106 ! include 'COMMON.DERIV'
1107 ! include 'COMMON.GEO'
1108 ! include 'COMMON.LOCAL'
1109 ! include 'COMMON.INTERACT'
1110 ! include 'COMMON.IOUNITS'
1111 ! include 'COMMON.NAMES'
1112 ! include 'COMMON.TIME1'
1113 ! include 'COMMON.MUCA'
1114 real(kind=8),dimension(3) :: vcm,incr
1115 real(kind=8),dimension(3) :: L
1116 integer :: count,rstcount !ilen,
1118 character(len=50) :: tytul
1119 integer :: maxcount_scale = 30
1120 !el common /gucio/ cm
1121 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1122 !el common /stochcalc/ stochforcvec
1123 integer :: icount_scale,itime_scal,i,j,ifac_time,iretcode,itime
1125 real(kind=8) :: epdrift,tt0,fac_time
1127 if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres)) !(MAXRES6) maxres6=6*maxres
1133 else if (lang.eq.2 .or. lang.eq.3) then
1135 call stochastic_force(stochforcvec)
1138 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
1140 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1147 icount_scale=icount_scale+1
1148 ! write(iout,*) "icount_scale",icount_scale,scalel
1149 if (icount_scale.gt.maxcount_scale) then
1151 "ERROR: too many attempts at scaling down the time step. ",&
1152 "amax=",amax,"epdrift=",epdrift,&
1153 "damax=",damax,"edriftmax=",edriftmax,&
1157 call MPI_Abort(MPI_COMM_WORLD,IERROR,IERRCODE)
1161 ! First step of the velocity Verlet algorithm
1166 else if (lang.eq.3) then
1168 call sd_verlet1_ciccotti
1170 else if (lang.eq.1) then
1175 ! Build the chain from the newly calculated coordinates
1176 call chainbuild_cart
1177 if (rattle) call rattle1
1179 if (large.and. mod(itime,ntwe).eq.0) then
1180 write (iout,*) "Cartesian and internal coordinates: step 1"
1185 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),&
1186 (dc(j,i+nres),j=1,3)
1188 write (iout,*) "Accelerations"
1190 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1191 (d_a(j,i+nres),j=1,3)
1193 write (iout,*) "Velocities, step 1"
1195 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1196 (d_t(j,i+nres),j=1,3)
1205 ! Calculate energy and forces
1207 call etotal(potEcomp)
1208 ! AL 4/17/17: Reduce the steps if NaNs occurred.
1209 if (potEcomp(0).gt.0.99e18 .or. isnan(potEcomp(0)).gt.0) then
1210 call enerprint(potEcomp)
1212 if (icount_scale.gt.15) then
1213 write (iout,*) "Tu jest problem",potEcomp(0),d_time
1214 ! call gen_rand_conf(1,*335)
1215 ! call minim_dc(potEcomp(0),iretcode,100)
1218 ! call etotal(potEcomp)
1219 ! write(iout,*) "needed to repara,",potEcomp
1222 ! 335 write(iout,*) "Failed genrand"
1226 if (large.and. mod(itime,ntwe).eq.0) &
1227 call enerprint(potEcomp)
1230 t_etotal=t_etotal+MPI_Wtime()-tt0
1232 t_etotal=t_etotal+tcpu()-tt0
1235 potE=potEcomp(0)-potEcomp(20)
1237 ! Get the new accelerations
1240 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1242 t_enegrad=t_enegrad+tcpu()-tt0
1244 ! Determine maximum acceleration and scale down the timestep if needed
1246 amax=amax/(itime_scal+1)**2
1247 call predict_edrift(epdrift)
1248 ! write(iout,*) "amax=",amax,damax,epdrift,edriftmax,amax/(itime_scal+1)
1250 ! write (iout,*) "before enter if",scalel,icount_scale
1251 if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then
1252 ! write(iout,*) "I enter if"
1253 ! Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
1255 ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax)) &
1257 itime_scal=itime_scal+ifac_time
1258 ! fac_time=dmin1(damax/amax,0.5d0)
1259 fac_time=0.5d0**ifac_time
1260 d_time=d_time*fac_time
1261 if (lang.eq.2 .or. lang.eq.3) then
1263 ! write (iout,*) "Calling sd_verlet_setup: 1"
1264 ! Rescale the stochastic forces and recalculate or restore
1265 ! the matrices of tinker integrator
1266 if (itime_scal.gt.maxflag_stoch) then
1267 if (large) write (iout,'(a,i5,a)') &
1268 "Calculate matrices for stochastic step;",&
1269 " itime_scal ",itime_scal
1271 call sd_verlet_p_setup
1273 call sd_verlet_ciccotti_setup
1275 write (iout,'(2a,i3,a,i3,1h.)') &
1276 "Warning: cannot store matrices for stochastic",&
1277 " integration because the index",itime_scal,&
1278 " is greater than",maxflag_stoch
1279 write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",&
1280 " integration Langevin algorithm for better efficiency."
1281 else if (flag_stoch(itime_scal)) then
1282 if (large) write (iout,'(a,i5,a,l1)') &
1283 "Restore matrices for stochastic step; itime_scal ",&
1284 itime_scal," flag ",flag_stoch(itime_scal)
1287 pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
1288 afric_mat(i,j)=afric0_mat(i,j,itime_scal)
1289 vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
1290 prand_mat(i,j)=prand0_mat(i,j,itime_scal)
1291 vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
1292 vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
1296 if (large) write (iout,'(2a,i5,a,l1)') &
1297 "Calculate & store matrices for stochastic step;",&
1298 " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
1300 call sd_verlet_p_setup
1302 call sd_verlet_ciccotti_setup
1304 flag_stoch(ifac_time)=.true.
1307 pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
1308 afric0_mat(i,j,itime_scal)=afric_mat(i,j)
1309 vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
1310 prand0_mat(i,j,itime_scal)=prand_mat(i,j)
1311 vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
1312 vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
1316 fac_time=1.0d0/dsqrt(fac_time)
1318 stochforcvec(i)=fac_time*stochforcvec(i)
1321 else if (lang.eq.1) then
1322 ! Rescale the accelerations due to stochastic forces
1323 fac_time=1.0d0/dsqrt(fac_time)
1325 d_as_work(i)=d_as_work(i)*fac_time
1328 if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)') &
1329 "itime",itime," Timestep scaled down to ",&
1330 d_time," ifac_time",ifac_time," itime_scal",itime_scal
1332 ! Second step of the velocity Verlet algorithm
1337 else if (lang.eq.3) then
1339 call sd_verlet2_ciccotti
1341 else if (lang.eq.1) then
1346 if (rattle) call rattle2
1349 if (d_time.ne.d_time0) then
1352 if (lang.eq.2 .or. lang.eq.3) then
1353 if (large) write (iout,'(a)') &
1354 "Restore original matrices for stochastic step"
1355 ! write (iout,*) "Calling sd_verlet_setup: 2"
1356 ! Restore the matrices of tinker integrator if the time step has been restored
1359 pfric_mat(i,j)=pfric0_mat(i,j,0)
1360 afric_mat(i,j)=afric0_mat(i,j,0)
1361 vfric_mat(i,j)=vfric0_mat(i,j,0)
1362 prand_mat(i,j)=prand0_mat(i,j,0)
1363 vrand_mat1(i,j)=vrand0_mat1(i,j,0)
1364 vrand_mat2(i,j)=vrand0_mat2(i,j,0)
1372 ! Calculate the kinetic and the total energy and the kinetic temperature
1376 ! call kinetic1(EK1)
1377 ! write (iout,*) "step",itime," EK",EK," EK1",EK1
1379 ! Couple the system to Berendsen bath if needed
1380 if (tbf .and. lang.eq.0) then
1383 kinetic_T=2.0d0/(dimen3*Rb)*EK
1384 ! Backup the coordinates, velocities, and accelerations
1388 d_t_old(j,i)=d_t(j,i)
1389 d_a_old(j,i)=d_a(j,i)
1393 if (mod(itime,ntwe).eq.0 .and. large) then
1394 write (iout,*) "Velocities, step 2"
1396 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1397 (d_t(j,i+nres),j=1,3)
1402 end subroutine velverlet_step
1403 !-----------------------------------------------------------------------------
1404 subroutine RESPA_step(itime)
1405 !-------------------------------------------------------------------------------
1406 ! Perform a single RESPA step.
1407 !-------------------------------------------------------------------------------
1408 ! implicit real*8 (a-h,o-z)
1409 ! include 'DIMENSIONS'
1413 use control, only:tcpu
1415 ! use io_conf, only:cartprint
1418 integer :: IERROR,ERRCODE
1420 ! include 'COMMON.SETUP'
1421 ! include 'COMMON.CONTROL'
1422 ! include 'COMMON.VAR'
1423 ! include 'COMMON.MD'
1425 ! include 'COMMON.LANGEVIN'
1427 ! include 'COMMON.LANGEVIN.lang0'
1429 ! include 'COMMON.CHAIN'
1430 ! include 'COMMON.DERIV'
1431 ! include 'COMMON.GEO'
1432 ! include 'COMMON.LOCAL'
1433 ! include 'COMMON.INTERACT'
1434 ! include 'COMMON.IOUNITS'
1435 ! include 'COMMON.NAMES'
1436 ! include 'COMMON.TIME1'
1437 real(kind=8),dimension(0:n_ene) :: energia_short,energia_long
1438 real(kind=8),dimension(3) :: L,vcm,incr
1439 real(kind=8),dimension(3,0:2*nres) :: dc_old0,d_t_old0,d_a_old0 !(3,0:maxres2) maxres2=2*maxres
1440 logical :: PRINT_AMTS_MSG = .false.
1441 integer :: count,rstcount !ilen,
1443 character(len=50) :: tytul
1444 integer :: maxcount_scale = 10
1445 !el common /gucio/ cm
1446 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1447 !el common /stochcalc/ stochforcvec
1448 integer :: itt,i,j,itsplit,itime
1450 !el common /cipiszcze/ itt
1452 real(kind=8) :: epdrift,tt0,epdriftmax
1455 if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres)) !(MAXRES6) maxres6=6*maxres
1459 if (large.and. mod(itime,ntwe).eq.0) then
1460 write (iout,*) "***************** RESPA itime",itime
1461 write (iout,*) "Cartesian and internal coordinates: step 0"
1463 call pdbout(0.0d0,"cipiszcze",iout)
1465 write (iout,*) "Accelerations from long-range forces"
1467 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1468 (d_a(j,i+nres),j=1,3)
1470 write (iout,*) "Velocities, step 0"
1472 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1473 (d_t(j,i+nres),j=1,3)
1478 ! Perform the initial RESPA step (increment velocities)
1479 ! write (iout,*) "*********************** RESPA ini"
1482 if (mod(itime,ntwe).eq.0 .and. large) then
1483 write (iout,*) "Velocities, end"
1485 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1486 (d_t(j,i+nres),j=1,3)
1490 ! Compute the short-range forces
1496 ! 7/2/2009 commented out
1498 ! call etotal_short(energia_short)
1501 ! 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
1504 d_a(j,i)=d_a_short(j,i)
1508 if (large.and. mod(itime,ntwe).eq.0) then
1509 write (iout,*) "energia_short",energia_short(0)
1510 write (iout,*) "Accelerations from short-range forces"
1512 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1513 (d_a(j,i+nres),j=1,3)
1518 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1520 t_enegrad=t_enegrad+tcpu()-tt0
1525 d_t_old(j,i)=d_t(j,i)
1526 d_a_old(j,i)=d_a(j,i)
1529 ! 6/30/08 A-MTS: attempt at increasing the split number
1532 dc_old0(j,i)=dc_old(j,i)
1533 d_t_old0(j,i)=d_t_old(j,i)
1534 d_a_old0(j,i)=d_a_old(j,i)
1537 if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
1538 if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
1545 ! write (iout,*) "itime",itime," ntime_split",ntime_split
1546 ! Split the time step
1547 d_time=d_time0/ntime_split
1548 ! Perform the short-range RESPA steps (velocity Verlet increments of
1549 ! positions and velocities using short-range forces)
1550 ! write (iout,*) "*********************** RESPA split"
1551 do itsplit=1,ntime_split
1554 else if (lang.eq.2 .or. lang.eq.3) then
1556 call stochastic_force(stochforcvec)
1559 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
1561 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1566 ! First step of the velocity Verlet algorithm
1571 else if (lang.eq.3) then
1573 call sd_verlet1_ciccotti
1575 else if (lang.eq.1) then
1580 ! Build the chain from the newly calculated coordinates
1581 call chainbuild_cart
1582 if (rattle) call rattle1
1584 if (large.and. mod(itime,ntwe).eq.0) then
1585 write (iout,*) "***** ITSPLIT",itsplit
1586 write (iout,*) "Cartesian and internal coordinates: step 1"
1587 call pdbout(0.0d0,"cipiszcze",iout)
1590 write (iout,*) "Velocities, step 1"
1592 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1593 (d_t(j,i+nres),j=1,3)
1602 ! Calculate energy and forces
1604 call etotal_short(energia_short)
1605 ! AL 4/17/17: Exit itime_split loop when energy goes infinite
1606 if (energia_short(0).gt.0.99e20 .or. isnan(energia_short(0)) ) then
1607 if (PRINT_AMTS_MSG) &
1608 write (iout,*) "Infinities/NaNs in energia_short",energia_short(0),"; increasing ntime_split to",ntime_split
1609 ntime_split=ntime_split*2
1610 if (ntime_split.gt.maxtime_split) then
1613 "Cannot rescue the run; aborting job. Retry with a smaller time step"
1615 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1618 "Cannot rescue the run; terminating. Retry with a smaller time step"
1624 if (large.and. mod(itime,ntwe).eq.0) &
1625 call enerprint(energia_short)
1628 t_eshort=t_eshort+MPI_Wtime()-tt0
1630 t_eshort=t_eshort+tcpu()-tt0
1634 ! Get the new accelerations
1636 ! 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
1639 d_a_short(j,i)=d_a(j,i)
1643 if (large.and. mod(itime,ntwe).eq.0) then
1644 write (iout,*)"energia_short",energia_short(0)
1645 write (iout,*) "Accelerations from short-range forces"
1647 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1648 (d_a(j,i+nres),j=1,3)
1653 ! Determine maximum acceleration and scale down the timestep if needed
1655 amax=amax/ntime_split**2
1656 call predict_edrift(epdrift)
1657 if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0) &
1658 write (iout,*) "amax",amax," damax",damax,&
1659 " epdrift",epdrift," epdriftmax",epdriftmax
1660 ! Exit loop and try with increased split number if the change of
1661 ! acceleration is too big
1662 if (amax.gt.damax .or. epdrift.gt.edriftmax) then
1663 if (ntime_split.lt.maxtime_split) then
1665 ntime_split=ntime_split*2
1666 ! AL 4/17/17: We should exit the itime_split loop when acceleration change is too big
1670 dc_old(j,i)=dc_old0(j,i)
1671 d_t_old(j,i)=d_t_old0(j,i)
1672 d_a_old(j,i)=d_a_old0(j,i)
1675 if (PRINT_AMTS_MSG) then
1676 write (iout,*) "acceleration/energy drift too large",amax,&
1677 epdrift," split increased to ",ntime_split," itime",itime,&
1683 "Uh-hu. Bumpy landscape. Maximum splitting number",&
1685 " already reached!!! Trying to carry on!"
1689 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1691 t_enegrad=t_enegrad+tcpu()-tt0
1693 ! Second step of the velocity Verlet algorithm
1698 else if (lang.eq.3) then
1700 call sd_verlet2_ciccotti
1702 else if (lang.eq.1) then
1707 if (rattle) call rattle2
1708 ! Backup the coordinates, velocities, and accelerations
1712 d_t_old(j,i)=d_t(j,i)
1713 d_a_old(j,i)=d_a(j,i)
1720 ! Restore the time step
1722 ! Compute long-range forces
1729 call etotal_long(energia_long)
1730 if (energia_long(0).gt.0.99e20 .or. isnan(energia_long(0))) then
1733 "Infinitied/NaNs in energia_long, Aborting MPI job."
1735 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1737 write (iout,*) "Infinitied/NaNs in energia_long, terminating."
1741 if (large.and. mod(itime,ntwe).eq.0) &
1742 call enerprint(energia_long)
1745 t_elong=t_elong+MPI_Wtime()-tt0
1747 t_elong=t_elong+tcpu()-tt0
1753 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1755 t_enegrad=t_enegrad+tcpu()-tt0
1757 ! Compute accelerations from long-range forces
1759 if (large.and. mod(itime,ntwe).eq.0) then
1760 write (iout,*) "energia_long",energia_long(0)
1761 write (iout,*) "Cartesian and internal coordinates: step 2"
1763 call pdbout(0.0d0,"cipiszcze",iout)
1765 write (iout,*) "Accelerations from long-range forces"
1767 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1768 (d_a(j,i+nres),j=1,3)
1770 write (iout,*) "Velocities, step 2"
1772 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1773 (d_t(j,i+nres),j=1,3)
1777 ! Compute the final RESPA step (increment velocities)
1778 ! write (iout,*) "*********************** RESPA fin"
1780 ! Compute the complete potential energy
1782 potEcomp(i)=energia_short(i)+energia_long(i)
1784 potE=potEcomp(0)-potEcomp(20)
1785 ! potE=energia_short(0)+energia_long(0)
1788 ! Calculate the kinetic and the total energy and the kinetic temperature
1791 ! Couple the system to Berendsen bath if needed
1792 if (tbf .and. lang.eq.0) then
1795 kinetic_T=2.0d0/(dimen3*Rb)*EK
1796 ! Backup the coordinates, velocities, and accelerations
1798 if (mod(itime,ntwe).eq.0 .and. large) then
1799 write (iout,*) "Velocities, end"
1801 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1802 (d_t(j,i+nres),j=1,3)
1807 end subroutine RESPA_step
1808 !-----------------------------------------------------------------------------
1809 subroutine RESPA_vel
1810 ! First and last RESPA step (incrementing velocities using long-range
1813 ! implicit real*8 (a-h,o-z)
1814 ! include 'DIMENSIONS'
1815 ! include 'COMMON.CONTROL'
1816 ! include 'COMMON.VAR'
1817 ! include 'COMMON.MD'
1818 ! include 'COMMON.CHAIN'
1819 ! include 'COMMON.DERIV'
1820 ! include 'COMMON.GEO'
1821 ! include 'COMMON.LOCAL'
1822 ! include 'COMMON.INTERACT'
1823 ! include 'COMMON.IOUNITS'
1824 ! include 'COMMON.NAMES'
1825 integer :: i,j,inres,mnum
1828 d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
1832 d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
1837 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
1838 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
1839 .and.(mnum.ne.5)) then
1842 d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
1847 end subroutine RESPA_vel
1848 !-----------------------------------------------------------------------------
1850 ! Applying velocity Verlet algorithm - step 1 to coordinates
1852 ! implicit real*8 (a-h,o-z)
1853 ! include 'DIMENSIONS'
1854 ! include 'COMMON.CONTROL'
1855 ! include 'COMMON.VAR'
1856 ! include 'COMMON.MD'
1857 ! include 'COMMON.CHAIN'
1858 ! include 'COMMON.DERIV'
1859 ! include 'COMMON.GEO'
1860 ! include 'COMMON.LOCAL'
1861 ! include 'COMMON.INTERACT'
1862 ! include 'COMMON.IOUNITS'
1863 ! include 'COMMON.NAMES'
1864 real(kind=8) :: adt,adt2
1865 integer :: i,j,inres,mnum
1868 write (iout,*) "VELVERLET1 START: DC"
1870 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),&
1871 (dc(j,i+nres),j=1,3)
1875 adt=d_a_old(j,0)*d_time
1877 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1878 d_t_new(j,0)=d_t_old(j,0)+adt2
1879 d_t(j,0)=d_t_old(j,0)+adt
1883 adt=d_a_old(j,i)*d_time
1885 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1886 d_t_new(j,i)=d_t_old(j,i)+adt2
1887 d_t(j,i)=d_t_old(j,i)+adt
1892 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
1893 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
1894 .and.(mnum.ne.5)) then
1897 adt=d_a_old(j,inres)*d_time
1899 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1900 d_t_new(j,inres)=d_t_old(j,inres)+adt2
1901 d_t(j,inres)=d_t_old(j,inres)+adt
1906 write (iout,*) "VELVERLET1 END: DC"
1908 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),&
1909 (dc(j,i+nres),j=1,3)
1913 end subroutine verlet1
1914 !-----------------------------------------------------------------------------
1916 ! Step 2 of the velocity Verlet algorithm: update velocities
1918 ! implicit real*8 (a-h,o-z)
1919 ! include 'DIMENSIONS'
1920 ! include 'COMMON.CONTROL'
1921 ! include 'COMMON.VAR'
1922 ! include 'COMMON.MD'
1923 ! include 'COMMON.CHAIN'
1924 ! include 'COMMON.DERIV'
1925 ! include 'COMMON.GEO'
1926 ! include 'COMMON.LOCAL'
1927 ! include 'COMMON.INTERACT'
1928 ! include 'COMMON.IOUNITS'
1929 ! include 'COMMON.NAMES'
1930 integer :: i,j,inres,mnum
1933 d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
1937 d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
1942 ! iti=iabs(itype(i,mnum))
1943 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
1944 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
1945 .and.(mnum.ne.5)) then
1948 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
1953 end subroutine verlet2
1954 !-----------------------------------------------------------------------------
1955 subroutine sddir_precalc
1956 ! Applying velocity Verlet algorithm - step 1 to coordinates
1957 ! implicit real*8 (a-h,o-z)
1958 ! include 'DIMENSIONS'
1964 ! include 'COMMON.CONTROL'
1965 ! include 'COMMON.VAR'
1966 ! include 'COMMON.MD'
1968 ! include 'COMMON.LANGEVIN'
1970 ! include 'COMMON.LANGEVIN.lang0'
1972 ! include 'COMMON.CHAIN'
1973 ! include 'COMMON.DERIV'
1974 ! include 'COMMON.GEO'
1975 ! include 'COMMON.LOCAL'
1976 ! include 'COMMON.INTERACT'
1977 ! include 'COMMON.IOUNITS'
1978 ! include 'COMMON.NAMES'
1979 ! include 'COMMON.TIME1'
1980 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1981 !el common /stochcalc/ stochforcvec
1982 real(kind=8) :: time00
1984 ! Compute friction and stochastic forces
1989 time_fric=time_fric+MPI_Wtime()-time00
1991 call stochastic_force(stochforcvec)
1992 time_stoch=time_stoch+MPI_Wtime()-time00
1995 ! Compute the acceleration due to friction forces (d_af_work) and stochastic
1996 ! forces (d_as_work)
1998 call ginv_mult(fric_work, d_af_work)
1999 call ginv_mult(stochforcvec, d_as_work)
2001 end subroutine sddir_precalc
2002 !-----------------------------------------------------------------------------
2003 subroutine sddir_verlet1
2004 ! Applying velocity Verlet algorithm - step 1 to velocities
2007 ! implicit real*8 (a-h,o-z)
2008 ! include 'DIMENSIONS'
2009 ! include 'COMMON.CONTROL'
2010 ! include 'COMMON.VAR'
2011 ! include 'COMMON.MD'
2013 ! include 'COMMON.LANGEVIN'
2015 ! include 'COMMON.LANGEVIN.lang0'
2017 ! include 'COMMON.CHAIN'
2018 ! include 'COMMON.DERIV'
2019 ! include 'COMMON.GEO'
2020 ! include 'COMMON.LOCAL'
2021 ! include 'COMMON.INTERACT'
2022 ! include 'COMMON.IOUNITS'
2023 ! include 'COMMON.NAMES'
2024 ! Revised 3/31/05 AL: correlation between random contributions to
2025 ! position and velocity increments included.
2026 real(kind=8) :: sqrt13 = 0.57735026918962576451d0 ! 1/sqrt(3)
2027 real(kind=8) :: adt,adt2
2028 integer :: i,j,ind,inres,mnum
2030 ! Add the contribution from BOTH friction and stochastic force to the
2031 ! coordinates, but ONLY the contribution from the friction forces to velocities
2034 adt=(d_a_old(j,0)+d_af_work(j))*d_time
2035 adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
2036 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
2037 d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
2038 d_t(j,0)=d_t_old(j,0)+adt
2043 adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
2044 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
2045 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
2046 d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
2047 d_t(j,i)=d_t_old(j,i)+adt
2053 ! iti=iabs(itype(i,mnum))
2054 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2055 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2056 .and.(mnum.ne.5)) then
2059 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
2060 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
2061 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
2062 d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
2063 d_t(j,inres)=d_t_old(j,inres)+adt
2069 end subroutine sddir_verlet1
2070 !-----------------------------------------------------------------------------
2071 subroutine sddir_verlet2
2072 ! Calculating the adjusted velocities for accelerations
2075 ! implicit real*8 (a-h,o-z)
2076 ! include 'DIMENSIONS'
2077 ! include 'COMMON.CONTROL'
2078 ! include 'COMMON.VAR'
2079 ! include 'COMMON.MD'
2081 ! include 'COMMON.LANGEVIN'
2083 ! include 'COMMON.LANGEVIN.lang0'
2085 ! include 'COMMON.CHAIN'
2086 ! include 'COMMON.DERIV'
2087 ! include 'COMMON.GEO'
2088 ! include 'COMMON.LOCAL'
2089 ! include 'COMMON.INTERACT'
2090 ! include 'COMMON.IOUNITS'
2091 ! include 'COMMON.NAMES'
2092 real(kind=8),dimension(6*nres) :: stochforcvec,d_as_work1 !(MAXRES6) maxres6=6*maxres
2093 real(kind=8) :: cos60 = 0.5d0, sin60 = 0.86602540378443864676d0
2094 integer :: i,j,ind,inres,mnum
2095 ! Revised 3/31/05 AL: correlation between random contributions to
2096 ! position and velocity increments included.
2097 ! The correlation coefficients are calculated at low-friction limit.
2098 ! Also, friction forces are now not calculated with new velocities.
2100 ! call friction_force
2101 call stochastic_force(stochforcvec)
2103 ! Compute the acceleration due to friction forces (d_af_work) and stochastic
2104 ! forces (d_as_work)
2106 call ginv_mult(stochforcvec, d_as_work1)
2112 d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j)) &
2113 +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
2118 d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j)) &
2119 +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
2125 ! iti=iabs(itype(i,mnum))
2126 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2127 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2128 .and.(mnum.ne.5)) then
2131 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres) &
2132 +d_af_work(ind+j))+sin60*d_as_work(ind+j) &
2133 +cos60*d_as_work1(ind+j))*d_time
2139 end subroutine sddir_verlet2
2140 !-----------------------------------------------------------------------------
2141 subroutine max_accel
2143 ! Find the maximum difference in the accelerations of the the sites
2144 ! at the beginning and the end of the time step.
2147 ! implicit real*8 (a-h,o-z)
2148 ! include 'DIMENSIONS'
2149 ! include 'COMMON.CONTROL'
2150 ! include 'COMMON.VAR'
2151 ! include 'COMMON.MD'
2152 ! include 'COMMON.CHAIN'
2153 ! include 'COMMON.DERIV'
2154 ! include 'COMMON.GEO'
2155 ! include 'COMMON.LOCAL'
2156 ! include 'COMMON.INTERACT'
2157 ! include 'COMMON.IOUNITS'
2158 real(kind=8),dimension(3) :: aux,accel,accel_old
2159 real(kind=8) :: dacc
2163 ! aux(j)=d_a(j,0)-d_a_old(j,0)
2164 accel_old(j)=d_a_old(j,0)
2171 ! 7/3/08 changed to asymmetric difference
2173 ! accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
2174 accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
2175 accel(j)=accel(j)+0.5d0*d_a(j,i)
2176 ! if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
2177 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
2178 dacc=dabs(accel(j)-accel_old(j))
2179 ! write (iout,*) i,dacc
2180 if (dacc.gt.amax) amax=dacc
2188 accel_old(j)=d_a_old(j,0)
2193 accel_old(j)=accel_old(j)+d_a_old(j,1)
2194 accel(j)=accel(j)+d_a(j,1)
2199 ! iti=iabs(itype(i,mnum))
2200 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2201 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2202 .and.(mnum.ne.5)) then
2204 ! accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
2205 accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
2206 accel(j)=accel(j)+d_a(j,i+nres)
2210 ! if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
2211 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
2212 dacc=dabs(accel(j)-accel_old(j))
2213 ! write (iout,*) "side-chain",i,dacc
2214 if (dacc.gt.amax) amax=dacc
2218 accel_old(j)=accel_old(j)+d_a_old(j,i)
2219 accel(j)=accel(j)+d_a(j,i)
2220 ! aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
2224 end subroutine max_accel
2225 !-----------------------------------------------------------------------------
2226 subroutine predict_edrift(epdrift)
2228 ! Predict the drift of the potential energy
2231 use control_data, only: lmuca
2232 ! implicit real*8 (a-h,o-z)
2233 ! include 'DIMENSIONS'
2234 ! include 'COMMON.CONTROL'
2235 ! include 'COMMON.VAR'
2236 ! include 'COMMON.MD'
2237 ! include 'COMMON.CHAIN'
2238 ! include 'COMMON.DERIV'
2239 ! include 'COMMON.GEO'
2240 ! include 'COMMON.LOCAL'
2241 ! include 'COMMON.INTERACT'
2242 ! include 'COMMON.IOUNITS'
2243 ! include 'COMMON.MUCA'
2244 real(kind=8) :: epdrift,epdriftij
2246 ! Drift of the potential energy
2252 epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
2253 if (lmuca) epdriftij=epdriftij*factor
2254 ! write (iout,*) "back",i,j,epdriftij
2255 if (epdriftij.gt.epdrift) epdrift=epdriftij
2259 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1.and.&
2260 molnum(i).ne.5) then
2263 dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
2264 if (lmuca) epdriftij=epdriftij*factor
2265 ! write (iout,*) "side",i,j,epdriftij
2266 if (epdriftij.gt.epdrift) epdrift=epdriftij
2270 epdrift=0.5d0*epdrift*d_time*d_time
2271 ! write (iout,*) "epdrift",epdrift
2273 end subroutine predict_edrift
2274 !-----------------------------------------------------------------------------
2275 subroutine verlet_bath
2277 ! Coupling to the thermostat by using the Berendsen algorithm
2280 ! implicit real*8 (a-h,o-z)
2281 ! include 'DIMENSIONS'
2282 ! include 'COMMON.CONTROL'
2283 ! include 'COMMON.VAR'
2284 ! include 'COMMON.MD'
2285 ! include 'COMMON.CHAIN'
2286 ! include 'COMMON.DERIV'
2287 ! include 'COMMON.GEO'
2288 ! include 'COMMON.LOCAL'
2289 ! include 'COMMON.INTERACT'
2290 ! include 'COMMON.IOUNITS'
2291 ! include 'COMMON.NAMES'
2292 real(kind=8) :: T_half,fact
2293 integer :: i,j,inres,mnum
2295 T_half=2.0d0/(dimen3*Rb)*EK
2296 fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
2297 ! write(iout,*) "T_half", T_half
2298 ! write(iout,*) "EK", EK
2299 ! write(iout,*) "fact", fact
2301 d_t(j,0)=fact*d_t(j,0)
2305 d_t(j,i)=fact*d_t(j,i)
2310 ! iti=iabs(itype(i,mnum))
2311 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2312 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2313 .and.(mnum.ne.5)) then
2316 d_t(j,inres)=fact*d_t(j,inres)
2321 end subroutine verlet_bath
2322 !-----------------------------------------------------------------------------
2324 ! Set up the initial conditions of a MD simulation
2327 use control, only:tcpu
2328 !el use io_basic, only:ilen
2331 use minimm, only:minim_dc,minimize,sc_move
2332 use io_config, only:readrst
2333 use io, only:statout
2334 ! implicit real*8 (a-h,o-z)
2335 ! include 'DIMENSIONS'
2338 character(len=16) :: form
2339 integer :: IERROR,ERRCODE
2341 integer :: iranmin,itrial,itmp
2342 ! include 'COMMON.SETUP'
2343 ! include 'COMMON.CONTROL'
2344 ! include 'COMMON.VAR'
2345 ! include 'COMMON.MD'
2347 ! include 'COMMON.LANGEVIN'
2349 ! include 'COMMON.LANGEVIN.lang0'
2351 ! include 'COMMON.CHAIN'
2352 ! include 'COMMON.DERIV'
2353 ! include 'COMMON.GEO'
2354 ! include 'COMMON.LOCAL'
2355 ! include 'COMMON.INTERACT'
2356 ! include 'COMMON.IOUNITS'
2357 ! include 'COMMON.NAMES'
2358 ! include 'COMMON.REMD'
2359 real(kind=8),dimension(0:n_ene) :: energia_long,energia_short
2360 real(kind=8),dimension(3) :: vcm,incr,L
2361 real(kind=8) :: xv,sigv,lowb,highb
2362 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
2363 character(len=256) :: qstr
2366 character(len=50) :: tytul
2367 logical :: file_exist
2368 !el common /gucio/ cm
2369 integer :: i,j,ipos,iq,iw,nft_sc,iretcode,nfun,ierr,mnum,itime
2370 real(kind=8) :: etot,tt0
2374 ! write(iout,*) "d_time", d_time
2375 ! Compute the standard deviations of stochastic forces for Langevin dynamics
2376 ! if the friction coefficients do not depend on surface area
2377 if (lang.gt.0 .and. .not.surfarea) then
2380 stdforcp(i)=stdfp(mnum)*dsqrt(gamp(mnum))
2384 stdforcsc(i)=stdfsc(iabs(itype(i,mnum)),mnum) &
2385 *dsqrt(gamsc(iabs(itype(i,mnum)),mnum))
2388 ! Open the pdb file for snapshotshots
2391 if (ilen(tmpdir).gt.0) &
2392 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2393 liczba(:ilen(liczba))//".pdb")
2395 file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2399 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) &
2400 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2401 liczba(:ilen(liczba))//".x")
2402 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2405 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) &
2406 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2407 liczba(:ilen(liczba))//".cx")
2408 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2414 if (ilen(tmpdir).gt.0) &
2415 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
2416 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
2418 if (ilen(tmpdir).gt.0) &
2419 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
2420 cartname=prefix(:ilen(prefix))//"_MD.cx"
2424 write (qstr,'(256(1h ))')
2427 iq = qinfrag(i,iset)*10
2428 iw = wfrag(i,iset)/100
2430 if(me.eq.king.or..not.out1file) &
2431 write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
2432 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
2437 iq = qinpair(i,iset)*10
2438 iw = wpair(i,iset)/100
2440 if(me.eq.king.or..not.out1file) &
2441 write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
2442 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
2446 ! pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
2448 ! cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
2450 ! cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
2452 ! statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
2456 if (restart1file) then
2458 inquire(file=mremd_rst_name,exist=file_exist)
2459 write (*,*) me," Before broadcast: file_exist",file_exist
2461 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,&
2464 write (*,*) me," After broadcast: file_exist",file_exist
2465 ! inquire(file=mremd_rst_name,exist=file_exist)
2466 if(me.eq.king.or..not.out1file) &
2467 write(iout,*) "Initial state read by master and distributed"
2469 if (ilen(tmpdir).gt.0) &
2470 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_' &
2471 //liczba(:ilen(liczba))//'.rst')
2472 inquire(file=rest2name,exist=file_exist)
2475 if(.not.restart1file) then
2476 if(me.eq.king.or..not.out1file) &
2477 write(iout,*) "Initial state will be read from file ",&
2478 rest2name(:ilen(rest2name))
2481 call rescale_weights(t_bath)
2483 if(me.eq.king.or..not.out1file)then
2484 if (restart1file) then
2485 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),&
2488 write(iout,*) "File ",rest2name(:ilen(rest2name)),&
2491 write(iout,*) "Initial velocities randomly generated"
2498 ! Generate initial velocities
2499 if(me.eq.king.or..not.out1file) &
2500 write(iout,*) "Initial velocities randomly generated"
2505 ! rest2name = prefix(:ilen(prefix))//'.rst'
2506 if(me.eq.king.or..not.out1file)then
2507 write (iout,*) "Initial velocities"
2509 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2510 (d_t(j,i+nres),j=1,3)
2512 ! Zeroing the total angular momentum of the system
2513 write(iout,*) "Calling the zero-angular momentum subroutine"
2516 ! Getting the potential energy and forces and velocities and accelerations
2518 ! write (iout,*) "velocity of the center of the mass:"
2519 ! write (iout,*) (vcm(j),j=1,3)
2521 d_t(j,0)=d_t(j,0)-vcm(j)
2523 ! Removing the velocity of the center of mass
2525 if(me.eq.king.or..not.out1file)then
2526 write (iout,*) "vcm right after adjustment:"
2527 write (iout,*) (vcm(j),j=1,3)
2531 if(iranconf.ne.0 .or.indpdb.gt.0.and..not.unres_pdb .or.preminim) then
2533 print *, 'Calling OVERLAP_SC'
2534 call overlap_sc(fail)
2535 print *,'after OVERLAP'
2538 print *,'call SC_MOVE'
2539 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
2540 print *,'SC_move',nft_sc,etot
2541 if(me.eq.king.or..not.out1file) &
2542 write(iout,*) 'SC_move',nft_sc,etot
2546 print *, 'Calling MINIM_DC'
2547 call minim_dc(etot,iretcode,nfun)
2549 call geom_to_var(nvar,varia)
2550 print *,'Calling MINIMIZE.'
2551 call minimize(etot,varia,iretcode,nfun)
2552 call var_to_geom(nvar,varia)
2554 if(me.eq.king.or..not.out1file) &
2555 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2557 if(iranconf.ne.0) then
2558 !c 8/22/17 AL Loop to produce a low-energy random conformation
2561 if(me.eq.king.or..not.out1file) &
2562 write (iout,*) 'Calling OVERLAP_SC'
2563 call overlap_sc(fail)
2564 endif !endif overlap
2567 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
2568 print *,'SC_move',nft_sc,etot
2569 if(me.eq.king.or..not.out1file) &
2570 write(iout,*) 'SC_move',nft_sc,etot
2574 print *, 'Calling MINIM_DC'
2575 call minim_dc(etot,iretcode,nfun)
2576 call int_from_cart1(.false.)
2578 call geom_to_var(nvar,varia)
2579 print *,'Calling MINIMIZE.'
2580 call minimize(etot,varia,iretcode,nfun)
2581 call var_to_geom(nvar,varia)
2583 if(me.eq.king.or..not.out1file) &
2584 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2586 if (isnan(etot) .or. etot.gt.4.0d6) then
2587 write (iout,*) "Energy too large",etot, &
2588 " trying another random conformation"
2591 call gen_rand_conf(itmp,*30)
2593 30 write (iout,*) 'Failed to generate random conformation', &
2595 write (*,*) 'Processor:',me, &
2596 ' Failed to generate random conformation',&
2605 write (iout,'(a,i3,a)') 'Processor:',me, &
2606 ' error in generating random conformation.'
2607 write (*,'(a,i3,a)') 'Processor:',me, &
2608 ' error in generating random conformation.'
2611 ! call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
2612 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
2622 write (iout,'(a,i3,a)') 'Processor:',me, &
2623 ' failed to generate a low-energy random conformation.'
2624 write (*,'(a,i3,a,f10.3)') 'Processor:',me, &
2625 ' failed to generate a low-energy random conformation.',etot
2629 ! call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
2630 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
2637 call chainbuild_cart
2642 kinetic_T=2.0d0/(dimen3*Rb)*EK
2643 if(me.eq.king.or..not.out1file)then
2653 write(iout,*) "before ETOTAL"
2654 call etotal(potEcomp)
2655 if (large) call enerprint(potEcomp)
2658 t_etotal=t_etotal+MPI_Wtime()-tt0
2660 t_etotal=t_etotal+tcpu()-tt0
2667 if (amax*d_time .gt. dvmax) then
2668 d_time=d_time*dvmax/amax
2669 if(me.eq.king.or..not.out1file) write (iout,*) &
2670 "Time step reduced to",d_time,&
2671 " because of too large initial acceleration."
2673 if(me.eq.king.or..not.out1file)then
2674 write(iout,*) "Potential energy and its components"
2675 call enerprint(potEcomp)
2676 ! write(iout,*) (potEcomp(i),i=0,n_ene)
2678 potE=potEcomp(0)-potEcomp(20)
2682 if (ntwe.ne.0) call statout(itime)
2683 if(me.eq.king.or..not.out1file) &
2684 write (iout,'(/a/3(a25,1pe14.5/))') "Initial:", &
2685 " Kinetic energy",EK," Potential energy",potE, &
2686 " Total energy",totE," Maximum acceleration ", &
2689 write (iout,*) "Initial coordinates"
2691 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),&
2694 write (iout,*) "Initial dC"
2696 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),&
2697 (dc(j,i+nres),j=1,3)
2699 write (iout,*) "Initial velocities"
2700 write (iout,"(13x,' backbone ',23x,' side chain')")
2702 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2703 (d_t(j,i+nres),j=1,3)
2705 write (iout,*) "Initial accelerations"
2707 ! write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2708 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),&
2709 (d_a(j,i+nres),j=1,3)
2715 d_t_old(j,i)=d_t(j,i)
2716 d_a_old(j,i)=d_a(j,i)
2718 ! write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
2727 call etotal_short(energia_short)
2728 if (large) call enerprint(potEcomp)
2731 t_eshort=t_eshort+MPI_Wtime()-tt0
2733 t_eshort=t_eshort+tcpu()-tt0
2738 if(.not.out1file .and. large) then
2739 write (iout,*) "energia_long",energia_long(0),&
2740 " energia_short",energia_short(0),&
2741 " total",energia_long(0)+energia_short(0)
2742 write (iout,*) "Initial fast-force accelerations"
2744 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
2745 (d_a(j,i+nres),j=1,3)
2748 ! 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
2751 d_a_short(j,i)=d_a(j,i)
2760 call etotal_long(energia_long)
2761 if (large) call enerprint(potEcomp)
2764 t_elong=t_elong+MPI_Wtime()-tt0
2766 t_elong=t_elong+tcpu()-tt0
2771 if(.not.out1file .and. large) then
2772 write (iout,*) "energia_long",energia_long(0)
2773 write (iout,*) "Initial slow-force accelerations"
2775 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
2776 (d_a(j,i+nres),j=1,3)
2780 t_enegrad=t_enegrad+MPI_Wtime()-tt0
2782 t_enegrad=t_enegrad+tcpu()-tt0
2786 end subroutine init_MD
2787 !-----------------------------------------------------------------------------
2788 subroutine random_vel
2790 ! implicit real*8 (a-h,o-z)
2792 use random, only:anorm_distr
2794 ! include 'DIMENSIONS'
2795 ! include 'COMMON.CONTROL'
2796 ! include 'COMMON.VAR'
2797 ! include 'COMMON.MD'
2799 ! include 'COMMON.LANGEVIN'
2801 ! include 'COMMON.LANGEVIN.lang0'
2803 ! include 'COMMON.CHAIN'
2804 ! include 'COMMON.DERIV'
2805 ! include 'COMMON.GEO'
2806 ! include 'COMMON.LOCAL'
2807 ! include 'COMMON.INTERACT'
2808 ! include 'COMMON.IOUNITS'
2809 ! include 'COMMON.NAMES'
2810 ! include 'COMMON.TIME1'
2811 real(kind=8) :: xv,sigv,lowb,highb ,Ek1
2814 real(kind=8) ,allocatable, dimension(:) :: DDU1,DDU2,DL2,DL1,xsolv,DML,rs
2815 real(kind=8) :: sumx
2817 real(kind=8) ,allocatable, dimension(:) :: rsold
2818 real (kind=8),allocatable,dimension(:,:) :: matold
2822 integer :: i,j,ii,k,ind,mark,imark,mnum
2823 ! Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
2824 ! First generate velocities in the eigenspace of the G matrix
2825 ! write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
2832 sigv=dsqrt((Rb*t_bath)/geigen(i))
2835 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
2837 write (iout,*) "i",i," ii",ii," geigen",geigen(i),&
2838 " d_t_work_new",d_t_work_new(ii)
2849 Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
2852 write (iout,*) "Ek from eigenvectors",Ek1
2853 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
2857 ! Transform velocities to UNRES coordinate space
2858 allocate (DL1(2*nres))
2859 allocate (DDU1(2*nres))
2860 allocate (DL2(2*nres))
2861 allocate (DDU2(2*nres))
2862 allocate (xsolv(2*nres))
2863 allocate (DML(2*nres))
2864 allocate (rs(2*nres))
2866 allocate (rsold(2*nres))
2867 allocate (matold(2*nres,2*nres))
2869 matold(1,1)=DMorig(1)
2870 matold(1,2)=DU1orig(1)
2871 matold(1,3)=DU2orig(1)
2872 write (*,*) DMorig(1),DU1orig(1),DU2orig(1)
2877 matold(i,j)=DMorig(i)
2878 matold(i,j-1)=DU1orig(i-1)
2880 matold(i,j-2)=DU2orig(i-2)
2888 matold(i,j+1)=DU1orig(i)
2894 matold(i,j+2)=DU2orig(i)
2898 matold(dimen,dimen)=DMorig(dimen)
2899 matold(dimen,dimen-1)=DU1orig(dimen-1)
2900 matold(dimen,dimen-2)=DU2orig(dimen-2)
2901 write(iout,*) "old gmatrix"
2902 call matout(dimen,dimen,2*nres,2*nres,matold)
2906 ! Find the ith eigenvector of the pentadiagonal inertiq matrix
2910 DML(j)=DMorig(j)-geigen(i)
2913 DML(j-1)=DMorig(j)-geigen(i)
2918 DDU1(imark-1)=DU2orig(imark-1)
2919 do j=imark+1,dimen-1
2920 DDU1(j-1)=DU1orig(j)
2928 DDU2(j)=DU2orig(j+1)
2937 write (iout,*) "DL2,DL1,DML,DDU1,DDU2"
2938 write(iout,'(10f10.5)') (DL2(k),k=3,dimen-1)
2939 write(iout,'(10f10.5)') (DL1(k),k=2,dimen-1)
2940 write(iout,'(10f10.5)') (DML(k),k=1,dimen-1)
2941 write(iout,'(10f10.5)') (DDU1(k),k=1,dimen-2)
2942 write(iout,'(10f10.5)') (DDU2(k),k=1,dimen-3)
2945 if (imark.gt.2) rs(imark-2)=-DU2orig(imark-2)
2946 if (imark.gt.1) rs(imark-1)=-DU1orig(imark-1)
2947 if (imark.lt.dimen) rs(imark)=-DU1orig(imark)
2948 if (imark.lt.dimen-1) rs(imark+1)=-DU2orig(imark)
2952 ! write (iout,*) "Vector rs"
2954 ! write (iout,*) j,rs(j)
2957 call FDIAG(dimen-1,DL2,DL1,DML,DDU1,DDU2,rs,xsolv,mark)
2964 sumx=-geigen(i)*xsolv(j)
2966 sumx=sumx+matold(j,k)*xsolv(k)
2969 sumx=sumx+matold(j,k)*xsolv(k-1)
2971 write(iout,'(i5,3f10.5)') j,sumx,rsold(j),sumx-rsold(j)
2974 sumx=-geigen(i)*xsolv(j-1)
2976 sumx=sumx+matold(j,k)*xsolv(k)
2979 sumx=sumx+matold(j,k)*xsolv(k-1)
2981 write(iout,'(i5,3f10.5)') j-1,sumx,rsold(j-1),sumx-rsold(j-1)
2985 "Solution of equations system",i," for eigenvalue",geigen(i)
2987 write(iout,'(i5,f10.5)') j,xsolv(j)
2990 do j=dimen-1,imark,-1
2995 write (iout,*) "Un-normalized eigenvector",i," for eigenvalue",geigen(i)
2997 write(iout,'(i5,f10.5)') j,xsolv(j)
3000 ! Normalize ith eigenvector
3003 sumx=sumx+xsolv(j)**2
3007 xsolv(j)=xsolv(j)/sumx
3010 write (iout,*) "Eigenvector",i," for eigenvalue",geigen(i)
3012 write(iout,'(i5,3f10.5)') j,xsolv(j)
3015 ! All done at this point for eigenvector i, exit loop
3023 write (iout,*) "Unable to find eigenvector",i
3026 ! write (iout,*) "i=",i
3028 ! write (iout,*) "k=",k
3031 ! write(iout,*) "ind",ind," ind1",3*(i-1)+k
3032 d_t_work(ind)=d_t_work(ind) &
3033 +xsolv(j)*d_t_work_new(3*(i-1)+k)
3036 enddo ! i (loop over the eigenvectors)
3039 write (iout,*) "d_t_work"
3041 write (iout,"(i5,f10.5)") i,d_t_work(i)
3046 ! if (itype(i,1).eq.10) then
3048 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
3049 iti=iabs(itype(i,mnum))
3050 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3051 if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
3052 .or.(mnum.eq.5)) then
3059 ! write (iout,*) "k",k," ii+k",ii+k," ii+j+k",ii+j+k,"EK1",Ek1
3060 Ek1=Ek1+0.5d0*mp(mnum)*((d_t_work(ii+j+k)+d_t_work_new(ii+k))/2)**2+&
3061 0.5d0*0.25d0*IP(mnum)*(d_t_work(ii+j+k)-d_t_work_new(ii+k))**2
3064 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3065 .and.(mnum.ne.5)) ii=ii+3
3066 write (iout,*) "i",i," itype",itype(i,mnum)," mass",msc(itype(i,mnum),mnum)
3067 write (iout,*) "ii",ii
3070 write (iout,*) "k",k," ii",ii,"EK1",EK1
3071 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3073 Ek1=Ek1+0.5d0*Isc(iabs(itype(i,mnum)),mnum)*(d_t_work(ii)-d_t_work(ii-3))**2
3074 Ek1=Ek1+0.5d0*msc(iabs(itype(i,mnum)),mnum)*d_t_work(ii)**2
3076 write (iout,*) "i",i," ii",ii
3078 write (iout,*) "Ek from d_t_work",Ek1
3079 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
3081 if(allocated(DDU1)) deallocate(DDU1)
3082 if(allocated(DDU2)) deallocate(DDU2)
3083 if(allocated(DL2)) deallocate(DL2)
3084 if(allocated(DL1)) deallocate(DL1)
3085 if(allocated(xsolv)) deallocate(xsolv)
3086 if(allocated(DML)) deallocate(DML)
3087 if(allocated(rs)) deallocate(rs)
3089 if(allocated(matold)) deallocate(matold)
3090 if(allocated(rsold)) deallocate(rsold)
3095 d_t(k,j)=d_t_work(ind)
3099 if (itype(j,1).ne.10 .and. itype(j,mnum).ne.ntyp1_molec(mnum)&
3100 .and.(mnum.ne.5)) then
3102 d_t(k,j+nres)=d_t_work(ind)
3108 write (iout,*) "Random velocities in the Calpha,SC space"
3110 write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3)
3113 write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3)
3120 ! if (itype(i,1).eq.10) then
3122 if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
3123 .or.(mnum.eq.5)) then
3125 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
3129 d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
3130 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
3135 write (iout,*)"Random velocities in the virtual-bond-vector space"
3137 write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3)
3140 write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3)
3143 write (iout,*) "Ek from d_t_work",Ek1
3144 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
3152 d_t_work(ind)=d_t_work(ind) &
3153 +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
3155 ! write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
3159 ! Transfer to the d_t vector
3161 d_t(j,0)=d_t_work(j)
3167 d_t(j,i)=d_t_work(ind)
3172 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3173 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3174 .and.(mnum.ne.5)) then
3177 d_t(j,i+nres)=d_t_work(ind)
3183 ! write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",&
3184 ! 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
3186 ! write(iout,*) "end init MD"
3188 end subroutine random_vel
3189 !-----------------------------------------------------------------------------
3191 subroutine sd_verlet_p_setup
3192 ! Sets up the parameters of stochastic Verlet algorithm
3193 ! implicit real*8 (a-h,o-z)
3194 ! include 'DIMENSIONS'
3195 use control, only: tcpu
3200 ! include 'COMMON.CONTROL'
3201 ! include 'COMMON.VAR'
3202 ! include 'COMMON.MD'
3204 ! include 'COMMON.LANGEVIN'
3206 ! include 'COMMON.LANGEVIN.lang0'
3208 ! include 'COMMON.CHAIN'
3209 ! include 'COMMON.DERIV'
3210 ! include 'COMMON.GEO'
3211 ! include 'COMMON.LOCAL'
3212 ! include 'COMMON.INTERACT'
3213 ! include 'COMMON.IOUNITS'
3214 ! include 'COMMON.NAMES'
3215 ! include 'COMMON.TIME1'
3216 real(kind=8),dimension(6*nres) :: emgdt !(MAXRES6) maxres6=6*maxres
3217 real(kind=8) :: pterm,vterm,rho,rhoc,vsig
3218 real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
3219 prand_vec,vrand_vec1,vrand_vec2 !(MAXRES6) maxres6=6*maxres
3220 logical :: lprn = .false.
3221 real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
3222 real(kind=8) :: ktm,gdt,egdt,gdt2,gdt3,gdt4,gdt5,gdt6,gdt7,gdt8,&
3224 integer :: i,maxres2
3231 ! AL 8/17/04 Code adapted from tinker
3233 ! Get the frictional and random terms for stochastic dynamics in the
3234 ! eigenspace of mass-scaled UNRES friction matrix
3238 gdt = fricgam(i) * d_time
3240 ! Stochastic dynamics reduces to simple MD for zero friction
3242 if (gdt .le. zero) then
3243 pfric_vec(i) = 1.0d0
3244 vfric_vec(i) = d_time
3245 afric_vec(i) = 0.5d0 * d_time * d_time
3246 prand_vec(i) = 0.0d0
3247 vrand_vec1(i) = 0.0d0
3248 vrand_vec2(i) = 0.0d0
3250 ! Analytical expressions when friction coefficient is large
3253 if (gdt .ge. gdt_radius) then
3256 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
3257 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
3258 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
3259 vterm = 1.0d0 - egdt**2
3260 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
3262 ! Use series expansions when friction coefficient is small
3273 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0 &
3274 - gdt5/120.0d0 + gdt6/720.0d0 &
3275 - gdt7/5040.0d0 + gdt8/40320.0d0 &
3276 - gdt9/362880.0d0) / fricgam(i)**2
3277 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
3278 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
3279 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0 &
3280 + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0 &
3281 + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0 &
3282 + 127.0d0*gdt9/90720.0d0
3283 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0 &
3284 - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0 &
3285 - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0 &
3286 - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
3287 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0 &
3288 - 17.0d0*gdt2/1280.0d0 &
3289 + 17.0d0*gdt3/6144.0d0 &
3290 + 40967.0d0*gdt4/34406400.0d0 &
3291 - 57203.0d0*gdt5/275251200.0d0 &
3292 - 1429487.0d0*gdt6/13212057600.0d0)
3295 ! Compute the scaling factors of random terms for the nonzero friction case
3297 ktm = 0.5d0*d_time/fricgam(i)
3298 psig = dsqrt(ktm*pterm) / fricgam(i)
3299 vsig = dsqrt(ktm*vterm)
3300 rhoc = dsqrt(1.0d0 - rho*rho)
3302 vrand_vec1(i) = vsig * rho
3303 vrand_vec2(i) = vsig * rhoc
3308 "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",&
3311 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),&
3312 afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
3316 ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
3319 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
3320 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
3321 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
3322 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
3323 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
3324 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
3327 t_sdsetup=t_sdsetup+MPI_Wtime()
3329 t_sdsetup=t_sdsetup+tcpu()-tt0
3332 end subroutine sd_verlet_p_setup
3333 !-----------------------------------------------------------------------------
3334 subroutine eigtransf1(n,ndim,ab,d,c)
3338 real(kind=8) :: ab(ndim,ndim,n),c(ndim,n),d(ndim)
3344 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
3349 end subroutine eigtransf1
3350 !-----------------------------------------------------------------------------
3351 subroutine eigtransf(n,ndim,a,b,d,c)
3355 real(kind=8) :: a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
3361 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
3366 end subroutine eigtransf
3367 !-----------------------------------------------------------------------------
3368 subroutine sd_verlet1
3370 ! Applying stochastic velocity Verlet algorithm - step 1 to velocities
3372 ! implicit real*8 (a-h,o-z)
3373 ! include 'DIMENSIONS'
3374 ! include 'COMMON.CONTROL'
3375 ! include 'COMMON.VAR'
3376 ! include 'COMMON.MD'
3378 ! include 'COMMON.LANGEVIN'
3380 ! include 'COMMON.LANGEVIN.lang0'
3382 ! include 'COMMON.CHAIN'
3383 ! include 'COMMON.DERIV'
3384 ! include 'COMMON.GEO'
3385 ! include 'COMMON.LOCAL'
3386 ! include 'COMMON.INTERACT'
3387 ! include 'COMMON.IOUNITS'
3388 ! include 'COMMON.NAMES'
3389 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
3390 !el common /stochcalc/ stochforcvec
3391 logical :: lprn = .false.
3392 real(kind=8) :: ddt1,ddt2
3393 integer :: i,j,ind,inres
3395 ! write (iout,*) "dc_old"
3397 ! write (iout,'(i5,3f10.5,5x,3f10.5)')
3398 ! & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3401 dc_work(j)=dc_old(j,0)
3402 d_t_work(j)=d_t_old(j,0)
3403 d_a_work(j)=d_a_old(j,0)
3408 dc_work(ind+j)=dc_old(j,i)
3409 d_t_work(ind+j)=d_t_old(j,i)
3410 d_a_work(ind+j)=d_a_old(j,i)
3416 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3417 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3418 .and.(mnum.ne.5)) then
3420 dc_work(ind+j)=dc_old(j,i+nres)
3421 d_t_work(ind+j)=d_t_old(j,i+nres)
3422 d_a_work(ind+j)=d_a_old(j,i+nres)
3430 "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",&
3434 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
3435 vfric_mat(i,j),afric_mat(i,j),&
3436 prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3444 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) &
3445 +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3446 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3447 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3449 d_t_work_new(i)=ddt1+0.5d0*ddt2
3450 d_t_work(i)=ddt1+ddt2
3455 d_t(j,0)=d_t_work(j)
3460 dc(j,i)=dc_work(ind+j)
3461 d_t(j,i)=d_t_work(ind+j)
3467 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3468 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3469 .and.(mnum.ne.5)) then
3472 dc(j,inres)=dc_work(ind+j)
3473 d_t(j,inres)=d_t_work(ind+j)
3479 end subroutine sd_verlet1
3480 !-----------------------------------------------------------------------------
3481 subroutine sd_verlet2
3483 ! Calculating the adjusted velocities for accelerations
3485 ! implicit real*8 (a-h,o-z)
3486 ! include 'DIMENSIONS'
3487 ! include 'COMMON.CONTROL'
3488 ! include 'COMMON.VAR'
3489 ! include 'COMMON.MD'
3491 ! include 'COMMON.LANGEVIN'
3493 ! include 'COMMON.LANGEVIN.lang0'
3495 ! include 'COMMON.CHAIN'
3496 ! include 'COMMON.DERIV'
3497 ! include 'COMMON.GEO'
3498 ! include 'COMMON.LOCAL'
3499 ! include 'COMMON.INTERACT'
3500 ! include 'COMMON.IOUNITS'
3501 ! include 'COMMON.NAMES'
3502 !el real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV !(MAXRES6) maxres6=6*maxres
3503 real(kind=8),dimension(6*nres) :: stochforcvecV !(MAXRES6) maxres6=6*maxres
3504 !el common /stochcalc/ stochforcvec
3506 real(kind=8) :: ddt1,ddt2
3507 integer :: i,j,ind,inres
3508 ! Compute the stochastic forces which contribute to velocity change
3510 call stochastic_force(stochforcvecV)
3517 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3518 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+ &
3519 vrand_mat2(i,j)*stochforcvecV(j)
3521 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
3525 d_t(j,0)=d_t_work(j)
3530 d_t(j,i)=d_t_work(ind+j)
3536 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3537 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3538 .and.(mnum.ne.5)) then
3541 d_t(j,inres)=d_t_work(ind+j)
3547 end subroutine sd_verlet2
3548 !-----------------------------------------------------------------------------
3549 subroutine sd_verlet_ciccotti_setup
3551 ! Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
3553 ! implicit real*8 (a-h,o-z)
3554 ! include 'DIMENSIONS'
3555 use control, only: tcpu
3560 ! include 'COMMON.CONTROL'
3561 ! include 'COMMON.VAR'
3562 ! include 'COMMON.MD'
3564 ! include 'COMMON.LANGEVIN'
3566 ! include 'COMMON.LANGEVIN.lang0'
3568 ! include 'COMMON.CHAIN'
3569 ! include 'COMMON.DERIV'
3570 ! include 'COMMON.GEO'
3571 ! include 'COMMON.LOCAL'
3572 ! include 'COMMON.INTERACT'
3573 ! include 'COMMON.IOUNITS'
3574 ! include 'COMMON.NAMES'
3575 ! include 'COMMON.TIME1'
3576 real(kind=8),dimension(6*nres) :: emgdt !(MAXRES6) maxres6=6*maxres
3577 real(kind=8) :: pterm,vterm,rho,rhoc,vsig
3578 real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
3579 prand_vec,vrand_vec1,vrand_vec2 !(MAXRES6) maxres6=6*maxres
3580 logical :: lprn = .false.
3581 real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
3582 real(kind=8) :: ktm,gdt,egdt,tt0
3583 integer :: i,maxres2
3590 ! AL 8/17/04 Code adapted from tinker
3592 ! Get the frictional and random terms for stochastic dynamics in the
3593 ! eigenspace of mass-scaled UNRES friction matrix
3597 write (iout,*) "i",i," fricgam",fricgam(i)
3598 gdt = fricgam(i) * d_time
3600 ! Stochastic dynamics reduces to simple MD for zero friction
3602 if (gdt .le. zero) then
3603 pfric_vec(i) = 1.0d0
3604 vfric_vec(i) = d_time
3605 afric_vec(i) = 0.5d0*d_time*d_time
3606 prand_vec(i) = afric_vec(i)
3607 vrand_vec2(i) = vfric_vec(i)
3609 ! Analytical expressions when friction coefficient is large
3614 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
3615 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
3616 prand_vec(i) = afric_vec(i)
3617 vrand_vec2(i) = vfric_vec(i)
3619 ! Compute the scaling factors of random terms for the nonzero friction case
3621 ! ktm = 0.5d0*d_time/fricgam(i)
3622 ! psig = dsqrt(ktm*pterm) / fricgam(i)
3623 ! vsig = dsqrt(ktm*vterm)
3624 ! prand_vec(i) = psig*afric_vec(i)
3625 ! vrand_vec2(i) = vsig*vfric_vec(i)
3630 "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",&
3633 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),&
3634 afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
3638 ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
3640 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
3641 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
3642 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
3643 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
3644 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
3646 t_sdsetup=t_sdsetup+MPI_Wtime()
3648 t_sdsetup=t_sdsetup+tcpu()-tt0
3651 end subroutine sd_verlet_ciccotti_setup
3652 !-----------------------------------------------------------------------------
3653 subroutine sd_verlet1_ciccotti
3655 ! Applying stochastic velocity Verlet algorithm - step 1 to velocities
3656 ! implicit real*8 (a-h,o-z)
3658 ! include 'DIMENSIONS'
3662 ! include 'COMMON.CONTROL'
3663 ! include 'COMMON.VAR'
3664 ! include 'COMMON.MD'
3666 ! include 'COMMON.LANGEVIN'
3668 ! include 'COMMON.LANGEVIN.lang0'
3670 ! include 'COMMON.CHAIN'
3671 ! include 'COMMON.DERIV'
3672 ! include 'COMMON.GEO'
3673 ! include 'COMMON.LOCAL'
3674 ! include 'COMMON.INTERACT'
3675 ! include 'COMMON.IOUNITS'
3676 ! include 'COMMON.NAMES'
3677 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
3678 !el common /stochcalc/ stochforcvec
3679 logical :: lprn = .false.
3680 real(kind=8) :: ddt1,ddt2
3681 integer :: i,j,ind,inres
3682 ! write (iout,*) "dc_old"
3684 ! write (iout,'(i5,3f10.5,5x,3f10.5)')
3685 ! & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3688 dc_work(j)=dc_old(j,0)
3689 d_t_work(j)=d_t_old(j,0)
3690 d_a_work(j)=d_a_old(j,0)
3695 dc_work(ind+j)=dc_old(j,i)
3696 d_t_work(ind+j)=d_t_old(j,i)
3697 d_a_work(ind+j)=d_a_old(j,i)
3702 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3704 dc_work(ind+j)=dc_old(j,i+nres)
3705 d_t_work(ind+j)=d_t_old(j,i+nres)
3706 d_a_work(ind+j)=d_a_old(j,i+nres)
3715 "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",&
3719 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
3720 vfric_mat(i,j),afric_mat(i,j),&
3721 prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3729 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) &
3730 +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3731 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3732 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3734 d_t_work_new(i)=ddt1+0.5d0*ddt2
3735 d_t_work(i)=ddt1+ddt2
3740 d_t(j,0)=d_t_work(j)
3745 dc(j,i)=dc_work(ind+j)
3746 d_t(j,i)=d_t_work(ind+j)
3752 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3753 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3754 .and.(mnum.ne.5)) then
3757 dc(j,inres)=dc_work(ind+j)
3758 d_t(j,inres)=d_t_work(ind+j)
3764 end subroutine sd_verlet1_ciccotti
3765 !-----------------------------------------------------------------------------
3766 subroutine sd_verlet2_ciccotti
3768 ! Calculating the adjusted velocities for accelerations
3770 ! implicit real*8 (a-h,o-z)
3771 ! include 'DIMENSIONS'
3772 ! include 'COMMON.CONTROL'
3773 ! include 'COMMON.VAR'
3774 ! include 'COMMON.MD'
3776 ! include 'COMMON.LANGEVIN'
3778 ! include 'COMMON.LANGEVIN.lang0'
3780 ! include 'COMMON.CHAIN'
3781 ! include 'COMMON.DERIV'
3782 ! include 'COMMON.GEO'
3783 ! include 'COMMON.LOCAL'
3784 ! include 'COMMON.INTERACT'
3785 ! include 'COMMON.IOUNITS'
3786 ! include 'COMMON.NAMES'
3787 !el real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV !(MAXRES6) maxres6=6*maxres
3788 real(kind=8),dimension(6*nres) :: stochforcvecV !(MAXRES6) maxres6=6*maxres
3789 !el common /stochcalc/ stochforcvec
3790 real(kind=8) :: ddt1,ddt2
3791 integer :: i,j,ind,inres
3793 ! Compute the stochastic forces which contribute to velocity change
3795 call stochastic_force(stochforcvecV)
3802 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3803 ! ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
3804 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
3806 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
3810 d_t(j,0)=d_t_work(j)
3815 d_t(j,i)=d_t_work(ind+j)
3821 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3823 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3826 d_t(j,inres)=d_t_work(ind+j)
3832 end subroutine sd_verlet2_ciccotti
3834 !-----------------------------------------------------------------------------
3836 !-----------------------------------------------------------------------------
3837 subroutine inertia_tensor
3839 ! Calculating the intertia tensor for the entire protein in order to
3840 ! remove the perpendicular components of velocity matrix which cause
3841 ! the molecule to rotate.
3844 ! implicit real*8 (a-h,o-z)
3845 ! include 'DIMENSIONS'
3846 ! include 'COMMON.CONTROL'
3847 ! include 'COMMON.VAR'
3848 ! include 'COMMON.MD'
3849 ! include 'COMMON.CHAIN'
3850 ! include 'COMMON.DERIV'
3851 ! include 'COMMON.GEO'
3852 ! include 'COMMON.LOCAL'
3853 ! include 'COMMON.INTERACT'
3854 ! include 'COMMON.IOUNITS'
3855 ! include 'COMMON.NAMES'
3857 real(kind=8),dimension(3,3) :: Im,Imcp,eigvec,Id
3858 real(kind=8),dimension(3) :: pr,eigval,L,vp,vrot
3859 real(kind=8) :: M_SC,mag,mag2,M_PEP
3860 real(kind=8),dimension(3,0:nres) :: vpp !(3,0:MAXRES)
3861 real(kind=8),dimension(3) :: vs_p,pp,incr,v
3862 real(kind=8),dimension(3,3) :: pr1,pr2
3864 !el common /gucio/ cm
3865 integer :: iti,inres,i,j,k,mnum
3876 ! calculating the center of the mass of the protein
3880 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
3881 if (itype(i,mnum).eq.ntyp1_molec(mnum)) cycle
3882 M_PEP=M_PEP+mp(mnum)
3884 cm(j)=cm(j)+(c(j,i)+0.5d0*dc(j,i))*mp(mnum)
3893 if (mnum.eq.5) cycle
3894 iti=iabs(itype(i,mnum))
3895 M_SC=M_SC+msc(iabs(iti),mnum)
3898 cm(j)=cm(j)+msc(iabs(iti),mnum)*c(j,inres)
3902 cm(j)=cm(j)/(M_SC+M_PEP)
3907 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
3909 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
3911 Im(1,1)=Im(1,1)+mp(mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
3912 Im(1,2)=Im(1,2)-mp(mnum)*pr(1)*pr(2)
3913 Im(1,3)=Im(1,3)-mp(mnum)*pr(1)*pr(3)
3914 Im(2,3)=Im(2,3)-mp(mnum)*pr(2)*pr(3)
3915 Im(2,2)=Im(2,2)+mp(mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
3916 Im(3,3)=Im(3,3)+mp(mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
3921 iti=iabs(itype(i,mnum))
3922 if (mnum.eq.5) cycle
3925 pr(j)=c(j,inres)-cm(j)
3927 Im(1,1)=Im(1,1)+msc(iabs(iti),mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
3928 Im(1,2)=Im(1,2)-msc(iabs(iti),mnum)*pr(1)*pr(2)
3929 Im(1,3)=Im(1,3)-msc(iabs(iti),mnum)*pr(1)*pr(3)
3930 Im(2,3)=Im(2,3)-msc(iabs(iti),mnum)*pr(2)*pr(3)
3931 Im(2,2)=Im(2,2)+msc(iabs(iti),mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
3932 Im(3,3)=Im(3,3)+msc(iabs(iti),mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
3937 Im(1,1)=Im(1,1)+Ip(mnum)*(1-dc_norm(1,i)*dc_norm(1,i))* &
3938 vbld(i+1)*vbld(i+1)*0.25d0
3939 Im(1,2)=Im(1,2)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(2,i))* &
3940 vbld(i+1)*vbld(i+1)*0.25d0
3941 Im(1,3)=Im(1,3)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(3,i))* &
3942 vbld(i+1)*vbld(i+1)*0.25d0
3943 Im(2,3)=Im(2,3)+Ip(mnum)*(-dc_norm(2,i)*dc_norm(3,i))* &
3944 vbld(i+1)*vbld(i+1)*0.25d0
3945 Im(2,2)=Im(2,2)+Ip(mnum)*(1-dc_norm(2,i)*dc_norm(2,i))* &
3946 vbld(i+1)*vbld(i+1)*0.25d0
3947 Im(3,3)=Im(3,3)+Ip(mnum)*(1-dc_norm(3,i)*dc_norm(3,i))* &
3948 vbld(i+1)*vbld(i+1)*0.25d0
3954 ! if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
3955 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3956 .and.(mnum.ne.5)) then
3957 iti=iabs(itype(i,mnum))
3959 Im(1,1)=Im(1,1)+Isc(iti,mnum)*(1-dc_norm(1,inres)* &
3960 dc_norm(1,inres))*vbld(inres)*vbld(inres)
3961 Im(1,2)=Im(1,2)-Isc(iti,mnum)*(dc_norm(1,inres)* &
3962 dc_norm(2,inres))*vbld(inres)*vbld(inres)
3963 Im(1,3)=Im(1,3)-Isc(iti,mnum)*(dc_norm(1,inres)* &
3964 dc_norm(3,inres))*vbld(inres)*vbld(inres)
3965 Im(2,3)=Im(2,3)-Isc(iti,mnum)*(dc_norm(2,inres)* &
3966 dc_norm(3,inres))*vbld(inres)*vbld(inres)
3967 Im(2,2)=Im(2,2)+Isc(iti,mnum)*(1-dc_norm(2,inres)* &
3968 dc_norm(2,inres))*vbld(inres)*vbld(inres)
3969 Im(3,3)=Im(3,3)+Isc(iti,mnum)*(1-dc_norm(3,inres)* &
3970 dc_norm(3,inres))*vbld(inres)*vbld(inres)
3975 ! write(iout,*) "The angular momentum before adjustment:"
3976 ! write(iout,*) (L(j),j=1,3)
3982 ! Copying the Im matrix for the djacob subroutine
3990 ! Finding the eigenvectors and eignvalues of the inertia tensor
3991 call djacob(3,3,10000,1.0d-10,Imcp,eigvec,eigval)
3992 ! write (iout,*) "Eigenvalues & Eigenvectors"
3993 ! write (iout,'(5x,3f10.5)') (eigval(i),i=1,3)
3996 ! write (iout,'(i5,3f10.5)') i,(eigvec(i,j),j=1,3)
3998 ! Constructing the diagonalized matrix
4000 if (dabs(eigval(i)).gt.1.0d-15) then
4001 Id(i,i)=1.0d0/eigval(i)
4008 Imcp(i,j)=eigvec(j,i)
4014 pr1(i,j)=pr1(i,j)+Id(i,k)*Imcp(k,j)
4021 pr2(i,j)=pr2(i,j)+eigvec(i,k)*pr1(k,j)
4025 ! Calculating the total rotational velocity of the molecule
4028 vrot(i)=vrot(i)+pr2(i,j)*L(j)
4031 ! Resetting the velocities
4033 call vecpr(vrot(1),dc(1,i),vp)
4035 d_t(j,i)=d_t(j,i)-vp(j)
4040 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4041 .and.(mnum.ne.5)) then
4042 ! if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
4043 ! if(itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
4045 call vecpr(vrot(1),dc(1,inres),vp)
4047 d_t(j,inres)=d_t(j,inres)-vp(j)
4052 ! write(iout,*) "The angular momentum after adjustment:"
4053 ! write(iout,*) (L(j),j=1,3)
4056 end subroutine inertia_tensor
4057 !-----------------------------------------------------------------------------
4058 subroutine angmom(cm,L)
4061 ! implicit real*8 (a-h,o-z)
4062 ! include 'DIMENSIONS'
4063 ! include 'COMMON.CONTROL'
4064 ! include 'COMMON.VAR'
4065 ! include 'COMMON.MD'
4066 ! include 'COMMON.CHAIN'
4067 ! include 'COMMON.DERIV'
4068 ! include 'COMMON.GEO'
4069 ! include 'COMMON.LOCAL'
4070 ! include 'COMMON.INTERACT'
4071 ! include 'COMMON.IOUNITS'
4072 ! include 'COMMON.NAMES'
4073 real(kind=8) :: mscab
4074 real(kind=8),dimension(3) :: L,cm,pr,vp,vrot,incr,v,pp
4075 integer :: iti,inres,i,j,mnum
4076 ! Calculate the angular momentum
4085 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
4087 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
4090 v(j)=incr(j)+0.5d0*d_t(j,i)
4093 incr(j)=incr(j)+d_t(j,i)
4095 call vecpr(pr(1),v(1),vp)
4097 L(j)=L(j)+mp(mnum)*vp(j)
4101 pp(j)=0.5d0*d_t(j,i)
4103 call vecpr(pr(1),pp(1),vp)
4105 L(j)=L(j)+Ip(mnum)*vp(j)
4113 iti=iabs(itype(i,mnum))
4121 pr(j)=c(j,inres)-cm(j)
4123 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4124 .and.(mnum.ne.5)) then
4126 v(j)=incr(j)+d_t(j,inres)
4133 call vecpr(pr(1),v(1),vp)
4134 ! write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),&
4135 ! " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
4137 L(j)=L(j)+mscab*vp(j)
4139 ! write (iout,*) "L",(l(j),j=1,3)
4140 if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4141 .and.(mnum.ne.5)) then
4143 v(j)=incr(j)+d_t(j,inres)
4145 call vecpr(dc(1,inres),d_t(1,inres),vp)
4147 L(j)=L(j)+Isc(iti,mnum)*vp(j)
4151 incr(j)=incr(j)+d_t(j,i)
4155 end subroutine angmom
4156 !-----------------------------------------------------------------------------
4157 subroutine vcm_vel(vcm)
4160 ! implicit real*8 (a-h,o-z)
4161 ! include 'DIMENSIONS'
4162 ! include 'COMMON.VAR'
4163 ! include 'COMMON.MD'
4164 ! include 'COMMON.CHAIN'
4165 ! include 'COMMON.DERIV'
4166 ! include 'COMMON.GEO'
4167 ! include 'COMMON.LOCAL'
4168 ! include 'COMMON.INTERACT'
4169 ! include 'COMMON.IOUNITS'
4170 real(kind=8),dimension(3) :: vcm,vv
4171 real(kind=8) :: summas,amas
4181 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
4183 summas=summas+mp(mnum)
4185 vcm(j)=vcm(j)+mp(mnum)*(vv(j)+0.5d0*d_t(j,i))
4189 amas=msc(iabs(itype(i,mnum)),mnum)
4194 if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4195 .and.(mnum.ne.5)) then
4197 vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres))
4201 vcm(j)=vcm(j)+amas*vv(j)
4205 vv(j)=vv(j)+d_t(j,i)
4208 ! write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas
4210 vcm(j)=vcm(j)/summas
4213 end subroutine vcm_vel
4214 !-----------------------------------------------------------------------------
4216 !-----------------------------------------------------------------------------
4218 ! RATTLE algorithm for velocity Verlet - step 1, UNRES
4222 ! implicit real*8 (a-h,o-z)
4223 ! include 'DIMENSIONS'
4225 ! include 'COMMON.CONTROL'
4226 ! include 'COMMON.VAR'
4227 ! include 'COMMON.MD'
4229 ! include 'COMMON.LANGEVIN'
4231 ! include 'COMMON.LANGEVIN.lang0'
4233 ! include 'COMMON.CHAIN'
4234 ! include 'COMMON.DERIV'
4235 ! include 'COMMON.GEO'
4236 ! include 'COMMON.LOCAL'
4237 ! include 'COMMON.INTERACT'
4238 ! include 'COMMON.IOUNITS'
4239 ! include 'COMMON.NAMES'
4240 ! include 'COMMON.TIME1'
4241 !el real(kind=8) :: gginv(2*nres,2*nres),&
4242 !el gdc(3,2*nres,2*nres)
4243 real(kind=8) :: dC_uncor(3,2*nres) !,&
4244 !el real(kind=8) :: Cmat(2*nres,2*nres)
4245 real(kind=8) :: x(2*nres),xcorr(3,2*nres) !maxres2=2*maxres
4246 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
4247 !el common /przechowalnia/ nbond
4248 integer :: max_rattle = 5
4249 logical :: lprn = .false., lprn1 = .false., not_done
4250 real(kind=8) :: tol_rattle = 1.0d-5
4252 integer :: ii,i,j,jj,l,ind,ind1,nres2
4255 !el /common/ przechowalnia
4257 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4258 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4259 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4261 if (lprn) write (iout,*) "RATTLE1"
4265 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4266 .and.(mnum.ne.5)) nbond=nbond+1
4268 ! Make a folded form of the Ginv-matrix
4281 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4286 if (itype(k,1).ne.10 .and. itype(k,mnum).ne.ntyp1_molec(mnum)&
4287 .and.(mnum.ne.5)) then
4291 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4299 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4310 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4314 if (itype(k,1).ne.10) then
4318 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4326 write (iout,*) "Matrix GGinv"
4327 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
4333 if (iter.gt.max_rattle) then
4334 write (iout,*) "Error - too many iterations in RATTLE."
4337 ! Calculate the matrix C = GG**(-1) dC_old o dC
4342 dC_uncor(j,ind1)=dC(j,i)
4346 if (itype(i,1).ne.10) then
4349 dC_uncor(j,ind1)=dC(j,i+nres)
4358 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
4362 if (itype(k,1).ne.10) then
4365 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
4370 ! Calculate deviations from standard virtual-bond lengths
4374 x(ind)=vbld(i+1)**2-vbl**2
4377 if (itype(i,1).ne.10) then
4379 x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4383 write (iout,*) "Coordinates and violations"
4385 write(iout,'(i5,3f10.5,5x,e15.5)') &
4386 i,(dC_uncor(j,i),j=1,3),x(i)
4388 write (iout,*) "Velocities and violations"
4392 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4393 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4397 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4398 .and.(mnum.ne.5)) then
4401 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4402 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4403 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4406 ! write (iout,*) "gdc"
4408 ! write (iout,*) "i",i
4410 ! write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4416 if (dabs(x(i)).gt.xmax) then
4420 if (xmax.lt.tol_rattle) then
4424 ! Calculate the matrix of the system of equations
4429 Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
4434 write (iout,*) "Matrix Cmat"
4435 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4437 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4438 ! Add constraint term to positions
4445 xx = xx+x(ii)*gdc(j,ind,ii)
4449 d_t_new(j,i)=d_t_new(j,i)-xx/d_time
4453 if (itype(i,1).ne.10) then
4458 xx = xx+x(ii)*gdc(j,ind,ii)
4461 dC(j,i+nres)=dC(j,i+nres)-xx
4462 d_t_new(j,i+nres)=d_t_new(j,i+nres)-xx/d_time
4466 ! Rebuild the chain using the new coordinates
4467 call chainbuild_cart
4469 write (iout,*) "New coordinates, Lagrange multipliers,",&
4470 " and differences between actual and standard bond lengths"
4474 xx=vbld(i+1)**2-vbl**2
4475 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4476 i,(dC(j,i),j=1,3),x(ind),xx
4480 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4483 xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4484 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4485 i,(dC(j,i+nres),j=1,3),x(ind),xx
4488 write (iout,*) "Velocities and violations"
4492 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4493 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4496 if (itype(i,1).ne.10) then
4498 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4499 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4500 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4507 10 write (iout,*) "Error - singularity in solving the system",&
4508 " of equations for Lagrange multipliers."
4512 "RATTLE inactive; use -DRATTLE switch at compile time."
4515 end subroutine rattle1
4516 !-----------------------------------------------------------------------------
4518 ! RATTLE algorithm for velocity Verlet - step 2, UNRES
4522 ! implicit real*8 (a-h,o-z)
4523 ! include 'DIMENSIONS'
4525 ! include 'COMMON.CONTROL'
4526 ! include 'COMMON.VAR'
4527 ! include 'COMMON.MD'
4529 ! include 'COMMON.LANGEVIN'
4531 ! include 'COMMON.LANGEVIN.lang0'
4533 ! include 'COMMON.CHAIN'
4534 ! include 'COMMON.DERIV'
4535 ! include 'COMMON.GEO'
4536 ! include 'COMMON.LOCAL'
4537 ! include 'COMMON.INTERACT'
4538 ! include 'COMMON.IOUNITS'
4539 ! include 'COMMON.NAMES'
4540 ! include 'COMMON.TIME1'
4541 !el real(kind=8) :: gginv(2*nres,2*nres),&
4542 !el gdc(3,2*nres,2*nres)
4543 real(kind=8) :: dC_uncor(3,2*nres) !,&
4544 !el Cmat(2*nres,2*nres)
4545 real(kind=8) :: x(2*nres) !maxres2=2*maxres
4546 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
4547 !el common /przechowalnia/ nbond
4548 integer :: max_rattle = 5
4549 logical :: lprn = .false., lprn1 = .false., not_done
4550 real(kind=8) :: tol_rattle = 1.0d-5
4554 !el /common/ przechowalnia
4555 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4556 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4557 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4559 if (lprn) write (iout,*) "RATTLE2"
4560 if (lprn) write (iout,*) "Velocity correction"
4561 ! Calculate the matrix G dC
4567 gdc(j,i,ind)=GGinv(i,ind)*dC(j,k)
4572 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4573 .and.(mnum.ne.5)) then
4576 gdc(j,i,ind)=GGinv(i,ind)*dC(j,k+nres)
4582 ! write (iout,*) "gdc"
4584 ! write (iout,*) "i",i
4586 ! write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4590 ! Calculate the matrix of the system of equations
4597 Cmat(ind,j)=Cmat(ind,j)+dC(k,i)*gdc(k,ind,j)
4603 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4604 .and.(mnum.ne.5)) then
4609 Cmat(ind,j)=Cmat(ind,j)+dC(k,i+nres)*gdc(k,ind,j)
4614 ! Calculate the scalar product dC o d_t_new
4618 x(ind)=scalar(d_t(1,i),dC(1,i))
4622 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4623 .and.(mnum.ne.5)) then
4625 x(ind)=scalar(d_t(1,i+nres),dC(1,i+nres))
4629 write (iout,*) "Velocities and violations"
4633 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4634 i,ind,(d_t(j,i),j=1,3),x(ind)
4638 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4639 .and.(mnum.ne.5)) then
4641 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4642 i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind)
4648 if (dabs(x(i)).gt.xmax) then
4652 if (xmax.lt.tol_rattle) then
4657 write (iout,*) "Matrix Cmat"
4658 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4660 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4661 ! Add constraint term to velocities
4668 xx = xx+x(ii)*gdc(j,ind,ii)
4670 d_t(j,i)=d_t(j,i)-xx
4675 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4676 .and.(mnum.ne.5)) then
4681 xx = xx+x(ii)*gdc(j,ind,ii)
4683 d_t(j,i+nres)=d_t(j,i+nres)-xx
4689 "New velocities, Lagrange multipliers violations"
4693 if (lprn) write (iout,'(2i5,3f10.5,5x,2e15.5)') &
4694 i,ind,(d_t(j,i),j=1,3),x(ind),scalar(d_t(1,i),dC(1,i))
4698 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4701 write (iout,'(2i5,3f10.5,5x,2e15.5)') &
4702 i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind),&
4703 scalar(d_t(1,i+nres),dC(1,i+nres))
4709 10 write (iout,*) "Error - singularity in solving the system",&
4710 " of equations for Lagrange multipliers."
4714 "RATTLE inactive; use -DRATTLE option at compile time."
4717 end subroutine rattle2
4718 !-----------------------------------------------------------------------------
4719 subroutine rattle_brown
4720 ! RATTLE/LINCS algorithm for Brownian dynamics, UNRES
4724 ! implicit real*8 (a-h,o-z)
4725 ! include 'DIMENSIONS'
4727 ! include 'COMMON.CONTROL'
4728 ! include 'COMMON.VAR'
4729 ! include 'COMMON.MD'
4731 ! include 'COMMON.LANGEVIN'
4733 ! include 'COMMON.LANGEVIN.lang0'
4735 ! include 'COMMON.CHAIN'
4736 ! include 'COMMON.DERIV'
4737 ! include 'COMMON.GEO'
4738 ! include 'COMMON.LOCAL'
4739 ! include 'COMMON.INTERACT'
4740 ! include 'COMMON.IOUNITS'
4741 ! include 'COMMON.NAMES'
4742 ! include 'COMMON.TIME1'
4743 !el real(kind=8) :: gginv(2*nres,2*nres),&
4744 !el gdc(3,2*nres,2*nres)
4745 real(kind=8) :: dC_uncor(3,2*nres) !,&
4746 !el real(kind=8) :: Cmat(2*nres,2*nres)
4747 real(kind=8) :: x(2*nres) !maxres2=2*maxres
4748 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
4749 !el common /przechowalnia/ nbond
4750 integer :: max_rattle = 5
4751 logical :: lprn = .false., lprn1 = .false., not_done
4752 real(kind=8) :: tol_rattle = 1.0d-5
4756 !el /common/ przechowalnia
4757 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4758 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4759 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4762 if (lprn) write (iout,*) "RATTLE_BROWN"
4765 if (itype(i,1).ne.10) nbond=nbond+1
4767 ! Make a folded form of the Ginv-matrix
4780 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4784 if (itype(k,1).ne.10) then
4788 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4795 if (itype(i,1).ne.10) then
4805 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4809 if (itype(k,1).ne.10) then
4813 if (j.eq.1 .and. l.eq.1)GGinv(ii,jj)=fricmat(ind,ind1)
4821 write (iout,*) "Matrix GGinv"
4822 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
4828 if (iter.gt.max_rattle) then
4829 write (iout,*) "Error - too many iterations in RATTLE."
4832 ! Calculate the matrix C = GG**(-1) dC_old o dC
4837 dC_uncor(j,ind1)=dC(j,i)
4841 if (itype(i,1).ne.10) then
4844 dC_uncor(j,ind1)=dC(j,i+nres)
4853 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
4857 if (itype(k,1).ne.10) then
4860 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
4865 ! Calculate deviations from standard virtual-bond lengths
4869 x(ind)=vbld(i+1)**2-vbl**2
4872 if (itype(i,1).ne.10) then
4874 x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4878 write (iout,*) "Coordinates and violations"
4880 write(iout,'(i5,3f10.5,5x,e15.5)') &
4881 i,(dC_uncor(j,i),j=1,3),x(i)
4883 write (iout,*) "Velocities and violations"
4887 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4888 i,ind,(d_t(j,i),j=1,3),scalar(d_t(1,i),dC_old(1,i))
4891 if (itype(i,1).ne.10) then
4893 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4894 i+nres,ind,(d_t(j,i+nres),j=1,3),&
4895 scalar(d_t(1,i+nres),dC_old(1,i+nres))
4898 write (iout,*) "gdc"
4900 write (iout,*) "i",i
4902 write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4908 if (dabs(x(i)).gt.xmax) then
4912 if (xmax.lt.tol_rattle) then
4916 ! Calculate the matrix of the system of equations
4921 Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
4926 write (iout,*) "Matrix Cmat"
4927 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4929 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4930 ! Add constraint term to positions
4937 xx = xx+x(ii)*gdc(j,ind,ii)
4940 d_t(j,i)=d_t(j,i)+xx/d_time
4945 if (itype(i,1).ne.10) then
4950 xx = xx+x(ii)*gdc(j,ind,ii)
4953 d_t(j,i+nres)=d_t(j,i+nres)+xx/d_time
4954 dC(j,i+nres)=dC(j,i+nres)+xx
4958 ! Rebuild the chain using the new coordinates
4959 call chainbuild_cart
4961 write (iout,*) "New coordinates, Lagrange multipliers,",&
4962 " and differences between actual and standard bond lengths"
4966 xx=vbld(i+1)**2-vbl**2
4967 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4968 i,(dC(j,i),j=1,3),x(ind),xx
4971 if (itype(i,1).ne.10) then
4973 xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4974 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4975 i,(dC(j,i+nres),j=1,3),x(ind),xx
4978 write (iout,*) "Velocities and violations"
4982 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4983 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4986 if (itype(i,1).ne.10) then
4988 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4989 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4990 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4997 10 write (iout,*) "Error - singularity in solving the system",&
4998 " of equations for Lagrange multipliers."
5002 "RATTLE inactive; use -DRATTLE option at compile time"
5005 end subroutine rattle_brown
5006 !-----------------------------------------------------------------------------
5008 !-----------------------------------------------------------------------------
5009 subroutine friction_force
5014 ! implicit real*8 (a-h,o-z)
5015 ! include 'DIMENSIONS'
5016 ! include 'COMMON.VAR'
5017 ! include 'COMMON.CHAIN'
5018 ! include 'COMMON.DERIV'
5019 ! include 'COMMON.GEO'
5020 ! include 'COMMON.LOCAL'
5021 ! include 'COMMON.INTERACT'
5022 ! include 'COMMON.MD'
5024 ! include 'COMMON.LANGEVIN'
5026 ! include 'COMMON.LANGEVIN.lang0'
5028 ! include 'COMMON.IOUNITS'
5029 !el real(kind=8),dimension(6*nres) :: gamvec !(MAXRES6) maxres6=6*maxres
5030 !el common /syfek/ gamvec
5031 real(kind=8) :: vv(3),vvtot(3,nres),v_work(6*nres) !,&
5032 !el ginvfric(2*nres,2*nres) !maxres2=2*maxres
5033 !el common /przechowalnia/ ginvfric
5035 logical :: lprn = .false., checkmode = .false.
5036 integer :: i,j,ind,k,nres2,nres6,mnum
5040 if(.not.allocated(gamvec)) allocate(gamvec(nres6)) !(MAXRES6)
5041 if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
5049 d_t_work(j)=d_t(j,0)
5054 d_t_work(ind+j)=d_t(j,i)
5060 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5061 .and.(mnum.ne.5)) then
5063 d_t_work(ind+j)=d_t(j,i+nres)
5069 call fricmat_mult(d_t_work,fric_work)
5071 if (.not.checkmode) return
5074 write (iout,*) "d_t_work and fric_work"
5076 write (iout,'(i3,2e15.5)') i,d_t_work(i),fric_work(i)
5080 friction(j,0)=fric_work(j)
5085 friction(j,i)=fric_work(ind+j)
5091 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5092 .and.(mnum.ne.5)) then
5093 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
5095 friction(j,i+nres)=fric_work(ind+j)
5101 write(iout,*) "Friction backbone"
5103 write(iout,'(i5,3e15.5,5x,3e15.5)') &
5104 i,(friction(j,i),j=1,3),(d_t(j,i),j=1,3)
5106 write(iout,*) "Friction side chain"
5108 write(iout,'(i5,3e15.5,5x,3e15.5)') &
5109 i,(friction(j,i+nres),j=1,3),(d_t(j,i+nres),j=1,3)
5118 vvtot(j,i)=vv(j)+0.5d0*d_t(j,i)
5119 vvtot(j,i+nres)=vv(j)+d_t(j,i+nres)
5120 vv(j)=vv(j)+d_t(j,i)
5123 write (iout,*) "vvtot backbone and sidechain"
5125 write (iout,'(i5,3e15.5,5x,3e15.5)') i,(vvtot(j,i),j=1,3),&
5126 (vvtot(j,i+nres),j=1,3)
5131 v_work(ind+j)=vvtot(j,i)
5137 v_work(ind+j)=vvtot(j,i+nres)
5141 write (iout,*) "v_work gamvec and site-based friction forces"
5143 write (iout,'(i5,3e15.5)') i,v_work(i),gamvec(i),&
5147 ! fric_work1(i)=0.0d0
5149 ! fric_work1(i)=fric_work1(i)-A(j,i)*gamvec(j)*v_work(j)
5152 ! write (iout,*) "fric_work and fric_work1"
5154 ! write (iout,'(i5,2e15.5)') i,fric_work(i),fric_work1(i)
5160 ginvfric(i,j)=ginvfric(i,j)+ginv(i,k)*fricmat(k,j)
5164 write (iout,*) "ginvfric"
5166 write (iout,'(i5,100f8.3)') i,(ginvfric(i,j),j=1,dimen)
5168 write (iout,*) "symmetry check"
5171 write (iout,*) i,j,ginvfric(i,j)-ginvfric(j,i)
5176 end subroutine friction_force
5177 !-----------------------------------------------------------------------------
5178 subroutine setup_fricmat
5182 use control_data, only:time_Bcast
5183 use control, only:tcpu
5185 ! implicit real*8 (a-h,o-z)
5189 real(kind=8) :: time00
5191 ! include 'DIMENSIONS'
5192 ! include 'COMMON.VAR'
5193 ! include 'COMMON.CHAIN'
5194 ! include 'COMMON.DERIV'
5195 ! include 'COMMON.GEO'
5196 ! include 'COMMON.LOCAL'
5197 ! include 'COMMON.INTERACT'
5198 ! include 'COMMON.MD'
5199 ! include 'COMMON.SETUP'
5200 ! include 'COMMON.TIME1'
5201 ! integer licznik /0/
5204 ! include 'COMMON.LANGEVIN'
5206 ! include 'COMMON.LANGEVIN.lang0'
5208 ! include 'COMMON.IOUNITS'
5210 integer :: i,j,ind,ind1,m
5211 logical :: lprn = .false.
5212 real(kind=8) :: dtdi !el ,gamvec(2*nres)
5213 !el real(kind=8),dimension(2*nres,2*nres) :: ginvfric,fcopy
5214 ! real(kind=8),allocatable,dimension(:,:) :: fcopy
5215 !el real(kind=8),dimension(2*nres*(2*nres+1)/2) :: Ghalf !(mmaxres2) (mmaxres2=(maxres2*(maxres2+1)/2))
5216 !el common /syfek/ gamvec
5217 real(kind=8) :: work(8*2*nres)
5218 integer :: iwork(2*nres)
5219 !el common /przechowalnia/ ginvfric,Ghalf,fcopy
5220 integer :: ii,iti,k,l,nzero,nres2,nres6,ierr,mnum
5224 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
5225 if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
5226 if (fg_rank.ne.king) goto 10
5231 if(.not.allocated(gamvec)) allocate(gamvec(nres2)) !(MAXRES2)
5232 if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
5233 if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
5234 !el allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
5235 if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !maxres2=2*maxres
5237 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
5238 ! Zeroing out fricmat
5244 ! Load the friction coefficients corresponding to peptide groups
5249 gamvec(ind1)=gamp(mnum)
5251 ! Load the friction coefficients corresponding to side chains
5255 gamsc(ntyp1_molec(j),j)=1.0d0
5262 gamvec(ii)=gamsc(iabs(iti),mnum)
5264 if (surfarea) call sdarea(gamvec)
5266 ! write (iout,*) "Matrix A and vector gamma"
5268 ! write (iout,'(i2,$)') i
5270 ! write (iout,'(f4.1,$)') A(i,j)
5272 ! write (iout,'(f8.3)') gamvec(i)
5276 write (iout,*) "Vector gamvec"
5278 write (iout,'(i5,f10.5)') i, gamvec(i)
5282 ! The friction matrix
5287 dtdi=dtdi+A(j,k)*A(j,i)*gamvec(j)
5294 write (iout,'(//a)') "Matrix fricmat"
5295 call matout2(dimen,dimen,nres2,nres2,fricmat)
5297 if (lang.eq.2 .or. lang.eq.3) then
5298 ! Mass-scale the friction matrix if non-direct integration will be performed
5304 Ginvfric(i,j)=Ginvfric(i,j)+ &
5305 Gsqrm(i,k)*Gsqrm(l,j)*fricmat(k,l)
5310 ! Diagonalize the friction matrix
5315 Ghalf(ind)=Ginvfric(i,j)
5318 call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
5321 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
5322 " mass-scaled friction matrix"
5323 call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
5325 ! Precompute matrices for tinker stochastic integrator
5332 mt1(i,j)=mt1(i,j)+fricvec(k,i)*gsqrm(k,j)
5333 mt2(i,j)=mt2(i,j)+fricvec(k,i)*gsqrp(k,j)
5339 else if (lang.eq.4) then
5340 ! Diagonalize the friction matrix
5345 Ghalf(ind)=fricmat(i,j)
5348 call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
5351 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
5353 call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
5355 ! Determine the number of zero eigenvalues of the friction matrix
5356 nzero=max0(dimen-dimen1,0)
5357 ! do while (fricgam(nzero+1).le.1.0d-5 .and. nzero.lt.dimen)
5360 write (iout,*) "Number of zero eigenvalues:",nzero
5365 fricmat(i,j)=fricmat(i,j) &
5366 +fricvec(i,k)*fricvec(j,k)/fricgam(k)
5371 write (iout,'(//a)') "Generalized inverse of fricmat"
5372 call matout(dimen,dimen,nres6,nres6,fricmat)
5377 if (nfgtasks.gt.1) then
5378 if (fg_rank.eq.0) then
5379 ! The matching BROADCAST for fg processors is called in ERGASTULUM
5385 call MPI_Bcast(10,1,MPI_INTEGER,king,FG_COMM,IERROR)
5387 time_Bcast=time_Bcast+MPI_Wtime()-time00
5389 time_Bcast=time_Bcast+tcpu()-time00
5391 ! print *,"Processor",myrank,
5392 ! & " BROADCAST iorder in SETUP_FRICMAT"
5395 write (iout,*) "setup_fricmat licznik"!,licznik !sp
5401 ! Scatter the friction matrix
5402 call MPI_Scatterv(fricmat(1,1),nginv_counts(0),&
5403 nginv_start(0),MPI_DOUBLE_PRECISION,fcopy(1,1),&
5404 myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
5407 time_scatter=time_scatter+MPI_Wtime()-time00
5408 time_scatter_fmat=time_scatter_fmat+MPI_Wtime()-time00
5410 time_scatter=time_scatter+tcpu()-time00
5411 time_scatter_fmat=time_scatter_fmat+tcpu()-time00
5415 do j=1,2*my_ng_count
5416 fricmat(j,i)=fcopy(i,j)
5419 ! write (iout,*) "My chunk of fricmat"
5420 ! call MATOUT2(my_ng_count,dimen,maxres2,maxres2,fcopy)
5424 end subroutine setup_fricmat
5425 !-----------------------------------------------------------------------------
5426 subroutine stochastic_force(stochforcvec)
5429 use random, only:anorm_distr
5430 ! implicit real*8 (a-h,o-z)
5431 ! include 'DIMENSIONS'
5432 use control, only: tcpu
5437 ! include 'COMMON.VAR'
5438 ! include 'COMMON.CHAIN'
5439 ! include 'COMMON.DERIV'
5440 ! include 'COMMON.GEO'
5441 ! include 'COMMON.LOCAL'
5442 ! include 'COMMON.INTERACT'
5443 ! include 'COMMON.MD'
5444 ! include 'COMMON.TIME1'
5446 ! include 'COMMON.LANGEVIN'
5448 ! include 'COMMON.LANGEVIN.lang0'
5450 ! include 'COMMON.IOUNITS'
5452 real(kind=8) :: x,sig,lowb,highb
5453 real(kind=8) :: ff(3),force(3,0:2*nres),zeta2,lowb2
5454 real(kind=8) :: highb2,sig2,forcvec(6*nres),stochforcvec(6*nres)
5455 real(kind=8) :: time00
5456 logical :: lprn = .false.
5457 integer :: i,j,ind,mnum
5461 stochforc(j,i)=0.0d0
5471 ! Compute the stochastic forces acting on bodies. Store in force.
5477 force(j,i)=anorm_distr(x,sig,lowb,highb)
5485 force(j,i+nres)=anorm_distr(x,sig2,lowb2,highb2)
5489 time_fsample=time_fsample+MPI_Wtime()-time00
5491 time_fsample=time_fsample+tcpu()-time00
5493 ! Compute the stochastic forces acting on virtual-bond vectors.
5499 stochforc(j,i)=ff(j)+0.5d0*force(j,i)
5502 ff(j)=ff(j)+force(j,i)
5504 ! if (itype(i+1,1).ne.ntyp1) then
5506 if (itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
5508 stochforc(j,i)=stochforc(j,i)+force(j,i+nres+1)
5509 ff(j)=ff(j)+force(j,i+nres+1)
5514 stochforc(j,0)=ff(j)+force(j,nnt+nres)
5518 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5519 .and.(mnum.ne.5)) then
5520 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
5522 stochforc(j,i+nres)=force(j,i+nres)
5528 stochforcvec(j)=stochforc(j,0)
5533 stochforcvec(ind+j)=stochforc(j,i)
5539 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5540 .and.(mnum.ne.5)) then
5541 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
5543 stochforcvec(ind+j)=stochforc(j,i+nres)
5549 write (iout,*) "stochforcvec"
5551 write(iout,'(i5,e15.5)') i,stochforcvec(i)
5553 write(iout,*) "Stochastic forces backbone"
5555 write(iout,'(i5,3e15.5)') i,(stochforc(j,i),j=1,3)
5557 write(iout,*) "Stochastic forces side chain"
5559 write(iout,'(i5,3e15.5)') &
5560 i,(stochforc(j,i+nres),j=1,3)
5568 write (iout,*) i,ind
5570 forcvec(ind+j)=force(j,i)
5575 write (iout,*) i,ind
5577 forcvec(j+ind)=force(j,i+nres)
5582 write (iout,*) "forcvec"
5586 write (iout,'(2i3,2f10.5)') i,j,force(j,i),&
5593 write (iout,'(2i3,2f10.5)') i,j,force(j,i+nres),&
5602 end subroutine stochastic_force
5603 !-----------------------------------------------------------------------------
5604 subroutine sdarea(gamvec)
5606 ! Scale the friction coefficients according to solvent accessible surface areas
5607 ! Code adapted from TINKER
5611 ! implicit real*8 (a-h,o-z)
5612 ! include 'DIMENSIONS'
5613 ! include 'COMMON.CONTROL'
5614 ! include 'COMMON.VAR'
5615 ! include 'COMMON.MD'
5617 ! include 'COMMON.LANGEVIN'
5619 ! include 'COMMON.LANGEVIN.lang0'
5621 ! include 'COMMON.CHAIN'
5622 ! include 'COMMON.DERIV'
5623 ! include 'COMMON.GEO'
5624 ! include 'COMMON.LOCAL'
5625 ! include 'COMMON.INTERACT'
5626 ! include 'COMMON.IOUNITS'
5627 ! include 'COMMON.NAMES'
5628 real(kind=8),dimension(2*nres) :: radius,gamvec !(maxres2)
5629 real(kind=8),parameter :: twosix = 1.122462048309372981d0
5630 logical :: lprn = .false.
5631 real(kind=8) :: probe,area,ratio
5632 integer :: i,j,ind,iti,mnum
5634 ! determine new friction coefficients every few SD steps
5636 ! set the atomic radii to estimates of sigma values
5638 ! print *,"Entered sdarea"
5644 ! Load peptide group radii
5647 radius(i)=pstok(mnum)
5649 ! Load side chain radii
5653 radius(i+nres)=restok(iti,mnum)
5656 ! write (iout,*) "i",i," radius",radius(i)
5659 radius(i) = radius(i) / twosix
5660 if (radius(i) .ne. 0.0d0) radius(i) = radius(i) + probe
5663 ! scale atomic friction coefficients by accessible area
5665 if (lprn) write (iout,*) &
5666 "Original gammas, surface areas, scaling factors, new gammas, ",&
5667 "std's of stochastic forces"
5670 if (radius(i).gt.0.0d0) then
5671 call surfatom (i,area,radius)
5672 ratio = dmax1(area/(4.0d0*pi*radius(i)**2),1.0d-1)
5673 if (lprn) write (iout,'(i5,3f10.5,$)') &
5674 i,gamvec(ind+1),area,ratio
5677 gamvec(ind) = ratio * gamvec(ind)
5680 stdforcp(i)=stdfp(mnum)*dsqrt(gamvec(ind))
5681 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcp(i)
5685 if (radius(i+nres).gt.0.0d0) then
5686 call surfatom (i+nres,area,radius)
5687 ratio = dmax1(area/(4.0d0*pi*radius(i+nres)**2),1.0d-1)
5688 if (lprn) write (iout,'(i5,3f10.5,$)') &
5689 i,gamvec(ind+1),area,ratio
5692 gamvec(ind) = ratio * gamvec(ind)
5695 stdforcsc(i)=stdfsc(itype(i,mnum),mnum)*dsqrt(gamvec(ind))
5696 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcsc(i)
5701 end subroutine sdarea
5702 !-----------------------------------------------------------------------------
5704 !-----------------------------------------------------------------------------
5707 ! ###################################################
5708 ! ## COPYRIGHT (C) 1996 by Jay William Ponder ##
5709 ! ## All Rights Reserved ##
5710 ! ###################################################
5712 ! ################################################################
5714 ! ## subroutine surfatom -- exposed surface area of an atom ##
5716 ! ################################################################
5719 ! "surfatom" performs an analytical computation of the surface
5720 ! area of a specified atom; a simplified version of "surface"
5722 ! literature references:
5724 ! T. J. Richmond, "Solvent Accessible Surface Area and
5725 ! Excluded Volume in Proteins", Journal of Molecular Biology,
5728 ! L. Wesson and D. Eisenberg, "Atomic Solvation Parameters
5729 ! Applied to Molecular Dynamics of Proteins in Solution",
5730 ! Protein Science, 1, 227-235 (1992)
5732 ! variables and parameters:
5734 ! ir number of atom for which area is desired
5735 ! area accessible surface area of the atom
5736 ! radius radii of each of the individual atoms
5739 subroutine surfatom(ir,area,radius)
5741 ! implicit real*8 (a-h,o-z)
5742 ! include 'DIMENSIONS'
5744 ! include 'COMMON.GEO'
5745 ! include 'COMMON.IOUNITS'
5747 integer :: nsup,nstart_sup
5748 ! double precision c,dc,dc_old,d_c_work,xloc,xrot,dc_norm
5749 ! common /chain/ c(3,maxres2+2),dc(3,0:maxres2),dc_old(3,0:maxres2),
5750 ! & xloc(3,maxres),xrot(3,maxres),dc_norm(3,0:maxres2),
5751 ! & dc_work(MAXRES6),nres,nres0
5752 integer,parameter :: maxarc=300
5756 integer :: mi,ni,narc
5757 integer :: key(maxarc)
5758 integer :: intag(maxarc)
5759 integer :: intag1(maxarc)
5760 real(kind=8) :: area,arcsum
5761 real(kind=8) :: arclen,exang
5762 real(kind=8) :: delta,delta2
5763 real(kind=8) :: eps,rmove
5764 real(kind=8) :: xr,yr,zr
5765 real(kind=8) :: rr,rrsq
5766 real(kind=8) :: rplus,rminus
5767 real(kind=8) :: axx,axy,axz
5768 real(kind=8) :: ayx,ayy
5769 real(kind=8) :: azx,azy,azz
5770 real(kind=8) :: uxj,uyj,uzj
5771 real(kind=8) :: tx,ty,tz
5772 real(kind=8) :: txb,tyb,td
5773 real(kind=8) :: tr2,tr,txr,tyr
5774 real(kind=8) :: tk1,tk2
5775 real(kind=8) :: thec,the,t,tb
5776 real(kind=8) :: txk,tyk,tzk
5777 real(kind=8) :: t1,ti,tf,tt
5778 real(kind=8) :: txj,tyj,tzj
5779 real(kind=8) :: ccsq,cc,xysq
5780 real(kind=8) :: bsqk,bk,cosine
5781 real(kind=8) :: dsqj,gi,pix2
5782 real(kind=8) :: therk,dk,gk
5783 real(kind=8) :: risqk,rik
5784 real(kind=8) :: radius(2*nres) !(maxatm) (maxatm=maxres2)
5785 real(kind=8) :: ri(maxarc),risq(maxarc)
5786 real(kind=8) :: ux(maxarc),uy(maxarc),uz(maxarc)
5787 real(kind=8) :: xc(maxarc),yc(maxarc),zc(maxarc)
5788 real(kind=8) :: xc1(maxarc),yc1(maxarc),zc1(maxarc)
5789 real(kind=8) :: dsq(maxarc),bsq(maxarc)
5790 real(kind=8) :: dsq1(maxarc),bsq1(maxarc)
5791 real(kind=8) :: arci(maxarc),arcf(maxarc)
5792 real(kind=8) :: ex(maxarc),lt(maxarc),gr(maxarc)
5793 real(kind=8) :: b(maxarc),b1(maxarc),bg(maxarc)
5794 real(kind=8) :: kent(maxarc),kout(maxarc)
5795 real(kind=8) :: ther(maxarc)
5796 logical :: moved,top
5797 logical :: omit(maxarc)
5800 ! maxatm = 2*nres !maxres2 maxres2=2*maxres
5801 ! maxlight = 8*maxatm
5804 ! maxtors = 4*maxatm
5806 ! zero out the surface area for the sphere of interest
5809 ! write (2,*) "ir",ir," radius",radius(ir)
5810 if (radius(ir) .eq. 0.0d0) return
5812 ! set the overlap significance and connectivity shift
5816 delta2 = delta * delta
5821 ! store coordinates and radius of the sphere of interest
5829 ! initialize values of some counters and summations
5838 ! test each sphere to see if it overlaps the sphere of interest
5841 if (i.eq.ir .or. radius(i).eq.0.0d0) goto 30
5842 rplus = rr + radius(i)
5844 if (abs(tx) .ge. rplus) goto 30
5846 if (abs(ty) .ge. rplus) goto 30
5848 if (abs(tz) .ge. rplus) goto 30
5850 ! check for sphere overlap by testing distance against radii
5852 xysq = tx*tx + ty*ty
5853 if (xysq .lt. delta2) then
5860 if (rplus-cc .le. delta) goto 30
5861 rminus = rr - radius(i)
5863 ! check to see if sphere of interest is completely buried
5865 if (cc-abs(rminus) .le. delta) then
5866 if (rminus .le. 0.0d0) goto 170
5870 ! check for too many overlaps with sphere of interest
5872 if (io .ge. maxarc) then
5874 20 format (/,' SURFATOM -- Increase the Value of MAXARC')
5878 ! get overlap between current sphere and sphere of interest
5887 gr(io) = (ccsq+rplus*rminus) / (2.0d0*rr*b1(io))
5893 ! case where no other spheres overlap the sphere of interest
5896 area = 4.0d0 * pi * rrsq
5900 ! case where only one sphere overlaps the sphere of interest
5903 area = pix2 * (1.0d0 + gr(1))
5904 area = mod(area,4.0d0*pi) * rrsq
5908 ! case where many spheres intersect the sphere of interest;
5909 ! sort the intersecting spheres by their degree of overlap
5911 call sort2 (io,gr,key)
5914 intag(i) = intag1(k)
5923 ! get radius of each overlap circle on surface of the sphere
5928 risq(i) = rrsq - gi*gi
5929 ri(i) = sqrt(risq(i))
5930 ther(i) = 0.5d0*pi - asin(min(1.0d0,max(-1.0d0,gr(i))))
5933 ! find boundary of inaccessible area on sphere of interest
5936 if (.not. omit(k)) then
5943 ! check to see if J circle is intersecting K circle;
5944 ! get distance between circle centers and sum of radii
5947 if (omit(j)) goto 60
5948 cc = (txk*xc(j)+tyk*yc(j)+tzk*zc(j))/(bk*b(j))
5949 cc = acos(min(1.0d0,max(-1.0d0,cc)))
5950 td = therk + ther(j)
5952 ! check to see if circles enclose separate regions
5954 if (cc .ge. td) goto 60
5956 ! check for circle J completely inside circle K
5958 if (cc+ther(j) .lt. therk) goto 40
5960 ! check for circles that are essentially parallel
5962 if (cc .gt. delta) goto 50
5967 ! check to see if sphere of interest is completely buried
5970 if (pix2-cc .le. td) goto 170
5976 ! find T value of circle intersections
5979 if (omit(k)) goto 110
5994 ! rotation matrix elements
6006 if (.not. omit(j)) then
6011 ! rotate spheres so K vector colinear with z-axis
6013 uxj = txj*axx + tyj*axy - tzj*axz
6014 uyj = tyj*ayy - txj*ayx
6015 uzj = txj*azx + tyj*azy + tzj*azz
6016 cosine = min(1.0d0,max(-1.0d0,uzj/b(j)))
6017 if (acos(cosine) .lt. therk+ther(j)) then
6018 dsqj = uxj*uxj + uyj*uyj
6023 tr2 = risqk*dsqj - tb*tb
6029 ! get T values of intersection for K circle
6032 tb = min(1.0d0,max(-1.0d0,tb))
6034 if (tyb-txr .lt. 0.0d0) tk1 = pix2 - tk1
6036 tb = min(1.0d0,max(-1.0d0,tb))
6038 if (tyb+txr .lt. 0.0d0) tk2 = pix2 - tk2
6039 thec = (rrsq*uzj-gk*bg(j)) / (rik*ri(j)*b(j))
6040 if (abs(thec) .lt. 1.0d0) then
6042 else if (thec .ge. 1.0d0) then
6044 else if (thec .le. -1.0d0) then
6048 ! see if "tk1" is entry or exit point; check t=0 point;
6049 ! "ti" is exit point, "tf" is entry point
6051 cosine = min(1.0d0,max(-1.0d0, &
6052 (uzj*gk-uxj*rik)/(b(j)*rr)))
6053 if ((acos(cosine)-ther(j))*(tk2-tk1) .le. 0.0d0) then
6061 if (narc .ge. maxarc) then
6063 70 format (/,' SURFATOM -- Increase the Value',&
6067 if (tf .le. ti) then
6088 ! special case; K circle without intersections
6090 if (narc .le. 0) goto 90
6092 ! general case; sum up arclength and set connectivity code
6094 call sort2 (narc,arci,key)
6099 if (narc .gt. 1) then
6102 if (t .lt. arci(j)) then
6103 arcsum = arcsum + arci(j) - t
6104 exang = exang + ex(ni)
6106 if (jb .ge. maxarc) then
6108 80 format (/,' SURFATOM -- Increase the Value',&
6113 kent(jb) = maxarc*i + k
6115 kout(jb) = maxarc*k + i
6124 arcsum = arcsum + pix2 - t
6126 exang = exang + ex(ni)
6129 kent(jb) = maxarc*i + k
6131 kout(jb) = maxarc*k + i
6138 arclen = arclen + gr(k)*arcsum
6141 if (arclen .eq. 0.0d0) goto 170
6142 if (jb .eq. 0) goto 150
6144 ! find number of independent boundaries and check connectivity
6148 if (kout(k) .ne. 0) then
6155 if (m .eq. kent(ii)) then
6158 if (j .eq. jb) goto 150
6170 ! attempt to fix connectivity error by moving atom slightly
6174 140 format (/,' SURFATOM -- Connectivity Error at Atom',i6)
6183 ! compute the exposed surface area for the sphere of interest
6186 area = ib*pix2 + exang + arclen
6187 area = mod(area,4.0d0*pi) * rrsq
6189 ! attempt to fix negative area by moving atom slightly
6191 if (area .lt. 0.0d0) then
6194 160 format (/,' SURFATOM -- Negative Area at Atom',i6)
6205 end subroutine surfatom
6206 !----------------------------------------------------------------
6207 !----------------------------------------------------------------
6208 subroutine alloc_MD_arrays
6209 !EL Allocation of arrays used by MD module
6211 integer :: nres2,nres6
6214 !----------------------
6218 allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
6219 allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
6220 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
6221 allocate(fricvec(nres2,nres2))
6222 allocate(pfric_mat(nres2,nres2),vfric_mat(nres2,nres2))
6223 allocate(afric_mat(nres2,nres2),prand_mat(nres2,nres2))
6224 allocate(vrand_mat1(nres2,nres2),vrand_mat2(nres2,nres2)) !(MAXRES2,MAXRES2)
6225 allocate(pfric0_mat(nres2,nres2,0:maxflag_stoch))
6226 allocate(afric0_mat(nres2,nres2,0:maxflag_stoch))
6227 allocate(vfric0_mat(nres2,nres2,0:maxflag_stoch))
6228 allocate(prand0_mat(nres2,nres2,0:maxflag_stoch))
6229 allocate(vrand0_mat1(nres2,nres2,0:maxflag_stoch))
6230 allocate(vrand0_mat2(nres2,nres2,0:maxflag_stoch)) !(MAXRES2,MAXRES2,0:maxflag_stoch)
6231 allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
6233 allocate(mt1(nres2,nres2),mt2(nres2,nres2),mt3(nres2,nres2)) !(maxres2,maxres2)
6234 !----------------------
6236 ! commom.langevin.lang0
6238 allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
6239 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
6240 allocate(fricvec(nres2,nres2)) !(MAXRES2,MAXRES2)
6241 allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
6242 allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
6245 if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2))
6246 !----------------------
6247 ! commom.hairpin in CSA module
6248 !----------------------
6249 ! common.mce in MCM_MD module
6250 !----------------------
6252 ! common /mdgrad/ in module.energy
6253 ! common /back_constr/ in module.energy
6254 ! common /qmeas/ in module.energy
6257 allocate(potEcomp(0:n_ene+4)) !(0:n_ene+4)
6259 allocate(d_t(3,0:nres2),d_a(3,0:nres2),d_t_old(3,0:nres2)) !(3,0:MAXRES2)
6260 allocate(d_a_work(nres6)) !(6*MAXRES)
6262 allocate(DM(nres2),DU1(nres2),DU2(nres2))
6263 allocate(DMorig(nres2),DU1orig(nres2),DU2orig(nres2))
6265 write (iout,*) "Before A Allocation",nfgtasks-1
6267 allocate(Gmat(nres2,nres2),A(nres2,nres2))
6268 if(.not.allocated(Ginv)) allocate(Ginv(nres2,nres2)) !in control: ergastulum
6269 allocate(Gsqrp(nres2,nres2),Gsqrm(nres2,nres2),Gvec(nres2,nres2)) !(maxres2,maxres2)
6271 allocate(Geigen(nres2)) !(maxres2)
6272 if(.not.allocated(vtot)) allocate(vtot(nres2)) !(maxres2)
6273 ! common /inertia/ in io_conf: parmread
6274 ! real(kind=8),dimension(:),allocatable :: ISC,msc !(ntyp+1)
6275 ! common /langevin/in io read_MDpar
6276 ! real(kind=8),dimension(:),allocatable :: gamsc !(ntyp1)
6277 ! real(kind=8),dimension(:),allocatable :: stdfsc !(ntyp)
6278 ! in io_conf: parmread
6279 ! real(kind=8),dimension(:),allocatable :: restok !(ntyp+1)
6280 ! common /mdpmpi/ in control: ergastulum
6281 if(.not.allocated(ng_start)) allocate(ng_start(0:nfgtasks-1))
6282 if(.not.allocated(ng_counts)) allocate(ng_counts(0:nfgtasks-1))
6283 if(.not.allocated(nginv_counts)) allocate(nginv_counts(0:nfgtasks-1)) !(0:MaxProcs-1)
6284 if(.not.allocated(nginv_start)) allocate(nginv_start(0:nfgtasks)) !(0:MaxProcs)
6285 !----------------------
6286 ! common.muca in read_muca
6287 ! common /double_muca/
6288 ! real(kind=8) :: elow,ehigh,factor,hbin,factor_min
6289 ! real(kind=8),dimension(:),allocatable :: emuca,nemuca,&
6290 ! nemuca2,hist !(4*maxres)
6291 ! common /integer_muca/
6292 ! integer :: nmuca,imtime,muca_smooth
6294 ! real(kind=8),dimension(:),allocatable :: elowi,ehighi !(maxprocs)
6295 !----------------------
6297 ! common /mdgrad/ in module.energy
6298 ! common /back_constr/ in module.energy
6299 ! common /qmeas/ in module.energy
6303 allocate(d_t_work(nres6),d_t_work_new(nres6),d_af_work(nres6))
6304 allocate(d_as_work(nres6),kinetic_force(nres6)) !(MAXRES6)
6305 allocate(d_t_new(3,0:nres2),d_a_old(3,0:nres2),d_a_short(3,0:nres2)) !,d_a !(3,0:MAXRES2)
6306 allocate(stdforcp(nres),stdforcsc(nres)) !(MAXRES)
6307 !----------------------
6309 allocate(D_ban(nres6)) !(MAXRES6) maxres6=6*maxres
6310 ! common /stochcalc/ stochforcvec
6311 allocate(stochforcvec(nres6)) !(MAXRES6) maxres6=6*maxres
6312 !----------------------
6314 end subroutine alloc_MD_arrays
6315 !-----------------------------------------------------------------------------
6316 !-----------------------------------------------------------------------------