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 ! write(iout,*) "adt",adt,"ads",adt2
2062 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
2063 d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
2064 d_t(j,inres)=d_t_old(j,inres)+adt
2070 end subroutine sddir_verlet1
2071 !-----------------------------------------------------------------------------
2072 subroutine sddir_verlet2
2073 ! Calculating the adjusted velocities for accelerations
2076 ! implicit real*8 (a-h,o-z)
2077 ! include 'DIMENSIONS'
2078 ! include 'COMMON.CONTROL'
2079 ! include 'COMMON.VAR'
2080 ! include 'COMMON.MD'
2082 ! include 'COMMON.LANGEVIN'
2084 ! include 'COMMON.LANGEVIN.lang0'
2086 ! include 'COMMON.CHAIN'
2087 ! include 'COMMON.DERIV'
2088 ! include 'COMMON.GEO'
2089 ! include 'COMMON.LOCAL'
2090 ! include 'COMMON.INTERACT'
2091 ! include 'COMMON.IOUNITS'
2092 ! include 'COMMON.NAMES'
2093 real(kind=8),dimension(6*nres) :: stochforcvec,d_as_work1 !(MAXRES6) maxres6=6*maxres
2094 real(kind=8) :: cos60 = 0.5d0, sin60 = 0.86602540378443864676d0
2095 integer :: i,j,ind,inres,mnum
2096 ! Revised 3/31/05 AL: correlation between random contributions to
2097 ! position and velocity increments included.
2098 ! The correlation coefficients are calculated at low-friction limit.
2099 ! Also, friction forces are now not calculated with new velocities.
2101 ! call friction_force
2102 call stochastic_force(stochforcvec)
2104 ! Compute the acceleration due to friction forces (d_af_work) and stochastic
2105 ! forces (d_as_work)
2107 call ginv_mult(stochforcvec, d_as_work1)
2113 d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j)) &
2114 +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
2119 d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j)) &
2120 +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
2126 ! iti=iabs(itype(i,mnum))
2127 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2128 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2129 .and.(mnum.ne.5)) then
2132 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres) &
2133 +d_af_work(ind+j))+sin60*d_as_work(ind+j) &
2134 +cos60*d_as_work1(ind+j))*d_time
2140 end subroutine sddir_verlet2
2141 !-----------------------------------------------------------------------------
2142 subroutine max_accel
2144 ! Find the maximum difference in the accelerations of the the sites
2145 ! at the beginning and the end of the time step.
2148 ! implicit real*8 (a-h,o-z)
2149 ! include 'DIMENSIONS'
2150 ! include 'COMMON.CONTROL'
2151 ! include 'COMMON.VAR'
2152 ! include 'COMMON.MD'
2153 ! include 'COMMON.CHAIN'
2154 ! include 'COMMON.DERIV'
2155 ! include 'COMMON.GEO'
2156 ! include 'COMMON.LOCAL'
2157 ! include 'COMMON.INTERACT'
2158 ! include 'COMMON.IOUNITS'
2159 real(kind=8),dimension(3) :: aux,accel,accel_old
2160 real(kind=8) :: dacc
2164 ! aux(j)=d_a(j,0)-d_a_old(j,0)
2165 accel_old(j)=d_a_old(j,0)
2172 ! 7/3/08 changed to asymmetric difference
2174 ! accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
2175 accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
2176 accel(j)=accel(j)+0.5d0*d_a(j,i)
2177 ! if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
2178 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
2179 dacc=dabs(accel(j)-accel_old(j))
2180 ! write (iout,*) i,dacc
2181 if (dacc.gt.amax) amax=dacc
2189 accel_old(j)=d_a_old(j,0)
2194 accel_old(j)=accel_old(j)+d_a_old(j,1)
2195 accel(j)=accel(j)+d_a(j,1)
2200 ! iti=iabs(itype(i,mnum))
2201 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2202 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2203 .and.(mnum.ne.5)) then
2205 ! accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
2206 accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
2207 accel(j)=accel(j)+d_a(j,i+nres)
2211 ! if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
2212 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
2213 dacc=dabs(accel(j)-accel_old(j))
2214 ! write (iout,*) "side-chain",i,dacc
2215 if (dacc.gt.amax) amax=dacc
2219 accel_old(j)=accel_old(j)+d_a_old(j,i)
2220 accel(j)=accel(j)+d_a(j,i)
2221 ! aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
2225 end subroutine max_accel
2226 !-----------------------------------------------------------------------------
2227 subroutine predict_edrift(epdrift)
2229 ! Predict the drift of the potential energy
2232 use control_data, only: lmuca
2233 ! implicit real*8 (a-h,o-z)
2234 ! include 'DIMENSIONS'
2235 ! include 'COMMON.CONTROL'
2236 ! include 'COMMON.VAR'
2237 ! include 'COMMON.MD'
2238 ! include 'COMMON.CHAIN'
2239 ! include 'COMMON.DERIV'
2240 ! include 'COMMON.GEO'
2241 ! include 'COMMON.LOCAL'
2242 ! include 'COMMON.INTERACT'
2243 ! include 'COMMON.IOUNITS'
2244 ! include 'COMMON.MUCA'
2245 real(kind=8) :: epdrift,epdriftij
2247 ! Drift of the potential energy
2253 epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
2254 if (lmuca) epdriftij=epdriftij*factor
2255 ! write (iout,*) "back",i,j,epdriftij
2256 if (epdriftij.gt.epdrift) epdrift=epdriftij
2260 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1.and.&
2261 molnum(i).ne.5) then
2264 dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
2265 if (lmuca) epdriftij=epdriftij*factor
2266 ! write (iout,*) "side",i,j,epdriftij
2267 if (epdriftij.gt.epdrift) epdrift=epdriftij
2271 epdrift=0.5d0*epdrift*d_time*d_time
2272 ! write (iout,*) "epdrift",epdrift
2274 end subroutine predict_edrift
2275 !-----------------------------------------------------------------------------
2276 subroutine verlet_bath
2278 ! Coupling to the thermostat by using the Berendsen algorithm
2281 ! implicit real*8 (a-h,o-z)
2282 ! include 'DIMENSIONS'
2283 ! include 'COMMON.CONTROL'
2284 ! include 'COMMON.VAR'
2285 ! include 'COMMON.MD'
2286 ! include 'COMMON.CHAIN'
2287 ! include 'COMMON.DERIV'
2288 ! include 'COMMON.GEO'
2289 ! include 'COMMON.LOCAL'
2290 ! include 'COMMON.INTERACT'
2291 ! include 'COMMON.IOUNITS'
2292 ! include 'COMMON.NAMES'
2293 real(kind=8) :: T_half,fact
2294 integer :: i,j,inres,mnum
2296 T_half=2.0d0/(dimen3*Rb)*EK
2297 fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
2298 ! write(iout,*) "T_half", T_half
2299 ! write(iout,*) "EK", EK
2300 ! write(iout,*) "fact", fact
2302 d_t(j,0)=fact*d_t(j,0)
2306 d_t(j,i)=fact*d_t(j,i)
2311 ! iti=iabs(itype(i,mnum))
2312 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2313 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2314 .and.(mnum.ne.5)) then
2317 d_t(j,inres)=fact*d_t(j,inres)
2322 end subroutine verlet_bath
2323 !-----------------------------------------------------------------------------
2325 ! Set up the initial conditions of a MD simulation
2328 use control, only:tcpu
2329 !el use io_basic, only:ilen
2332 use minimm, only:minim_dc,minimize,sc_move
2333 use io_config, only:readrst
2334 use io, only:statout
2335 ! implicit real*8 (a-h,o-z)
2336 ! include 'DIMENSIONS'
2339 character(len=16) :: form
2340 integer :: IERROR,ERRCODE
2342 integer :: iranmin,itrial,itmp
2343 ! include 'COMMON.SETUP'
2344 ! include 'COMMON.CONTROL'
2345 ! include 'COMMON.VAR'
2346 ! include 'COMMON.MD'
2348 ! include 'COMMON.LANGEVIN'
2350 ! include 'COMMON.LANGEVIN.lang0'
2352 ! include 'COMMON.CHAIN'
2353 ! include 'COMMON.DERIV'
2354 ! include 'COMMON.GEO'
2355 ! include 'COMMON.LOCAL'
2356 ! include 'COMMON.INTERACT'
2357 ! include 'COMMON.IOUNITS'
2358 ! include 'COMMON.NAMES'
2359 ! include 'COMMON.REMD'
2360 real(kind=8),dimension(0:n_ene) :: energia_long,energia_short
2361 real(kind=8),dimension(3) :: vcm,incr,L
2362 real(kind=8) :: xv,sigv,lowb,highb
2363 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
2364 character(len=256) :: qstr
2367 character(len=50) :: tytul
2368 logical :: file_exist
2369 !el common /gucio/ cm
2370 integer :: i,j,ipos,iq,iw,nft_sc,iretcode,nfun,ierr,mnum,itime
2371 real(kind=8) :: etot,tt0
2375 ! write(iout,*) "d_time", d_time
2376 ! Compute the standard deviations of stochastic forces for Langevin dynamics
2377 ! if the friction coefficients do not depend on surface area
2378 if (lang.gt.0 .and. .not.surfarea) then
2381 stdforcp(i)=stdfp(mnum)*dsqrt(gamp(mnum))
2385 stdforcsc(i)=stdfsc(iabs(itype(i,mnum)),mnum) &
2386 *dsqrt(gamsc(iabs(itype(i,mnum)),mnum))
2390 ! Open the pdb file for snapshotshots
2393 if (ilen(tmpdir).gt.0) &
2394 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2395 liczba(:ilen(liczba))//".pdb")
2397 file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2401 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) &
2402 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2403 liczba(:ilen(liczba))//".x")
2404 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2407 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) &
2408 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2409 liczba(:ilen(liczba))//".cx")
2410 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2416 if (ilen(tmpdir).gt.0) &
2417 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
2418 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
2420 if (ilen(tmpdir).gt.0) &
2421 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
2422 cartname=prefix(:ilen(prefix))//"_MD.cx"
2426 write (qstr,'(256(1h ))')
2429 iq = qinfrag(i,iset)*10
2430 iw = wfrag(i,iset)/100
2432 if(me.eq.king.or..not.out1file) &
2433 write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
2434 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
2439 iq = qinpair(i,iset)*10
2440 iw = wpair(i,iset)/100
2442 if(me.eq.king.or..not.out1file) &
2443 write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
2444 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
2448 ! pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
2450 ! cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
2452 ! cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
2454 ! statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
2458 if (restart1file) then
2460 inquire(file=mremd_rst_name,exist=file_exist)
2461 write (*,*) me," Before broadcast: file_exist",file_exist
2463 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,&
2466 write (*,*) me," After broadcast: file_exist",file_exist
2467 ! inquire(file=mremd_rst_name,exist=file_exist)
2468 if(me.eq.king.or..not.out1file) &
2469 write(iout,*) "Initial state read by master and distributed"
2471 if (ilen(tmpdir).gt.0) &
2472 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_' &
2473 //liczba(:ilen(liczba))//'.rst')
2474 inquire(file=rest2name,exist=file_exist)
2477 if(.not.restart1file) then
2478 if(me.eq.king.or..not.out1file) &
2479 write(iout,*) "Initial state will be read from file ",&
2480 rest2name(:ilen(rest2name))
2483 call rescale_weights(t_bath)
2485 if(me.eq.king.or..not.out1file)then
2486 if (restart1file) then
2487 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),&
2490 write(iout,*) "File ",rest2name(:ilen(rest2name)),&
2493 write(iout,*) "Initial velocities randomly generated"
2500 ! Generate initial velocities
2501 if(me.eq.king.or..not.out1file) &
2502 write(iout,*) "Initial velocities randomly generated"
2507 ! rest2name = prefix(:ilen(prefix))//'.rst'
2508 if(me.eq.king.or..not.out1file)then
2509 write (iout,*) "Initial velocities"
2511 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2512 (d_t(j,i+nres),j=1,3)
2514 ! Zeroing the total angular momentum of the system
2515 write(iout,*) "Calling the zero-angular momentum subroutine"
2518 ! Getting the potential energy and forces and velocities and accelerations
2520 ! write (iout,*) "velocity of the center of the mass:"
2521 ! write (iout,*) (vcm(j),j=1,3)
2523 d_t(j,0)=d_t(j,0)-vcm(j)
2525 ! Removing the velocity of the center of mass
2527 if(me.eq.king.or..not.out1file)then
2528 write (iout,*) "vcm right after adjustment:"
2529 write (iout,*) (vcm(j),j=1,3)
2531 if ((.not.rest).or.(forceminim)) then
2532 if (forceminim) call chainbuild_cart
2533 if(iranconf.ne.0 .or.indpdb.gt.0.and..not.unres_pdb .or.preminim) then
2535 print *, 'Calling OVERLAP_SC'
2536 call overlap_sc(fail)
2537 print *,'after OVERLAP'
2540 print *,'call SC_MOVE'
2541 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
2542 print *,'SC_move',nft_sc,etot
2543 if(me.eq.king.or..not.out1file) &
2544 write(iout,*) 'SC_move',nft_sc,etot
2548 print *, 'Calling MINIM_DC'
2549 call minim_dc(etot,iretcode,nfun)
2551 call geom_to_var(nvar,varia)
2552 print *,'Calling MINIMIZE.'
2553 call minimize(etot,varia,iretcode,nfun)
2554 call var_to_geom(nvar,varia)
2556 if(me.eq.king.or..not.out1file) &
2557 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2559 if(iranconf.ne.0) then
2560 !c 8/22/17 AL Loop to produce a low-energy random conformation
2563 if(me.eq.king.or..not.out1file) &
2564 write (iout,*) 'Calling OVERLAP_SC'
2565 call overlap_sc(fail)
2566 endif !endif overlap
2569 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
2570 print *,'SC_move',nft_sc,etot
2571 if(me.eq.king.or..not.out1file) &
2572 write(iout,*) 'SC_move',nft_sc,etot
2576 print *, 'Calling MINIM_DC'
2577 call minim_dc(etot,iretcode,nfun)
2578 call int_from_cart1(.false.)
2580 call geom_to_var(nvar,varia)
2581 print *,'Calling MINIMIZE.'
2582 call minimize(etot,varia,iretcode,nfun)
2583 call var_to_geom(nvar,varia)
2585 if(me.eq.king.or..not.out1file) &
2586 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2588 if (isnan(etot) .or. etot.gt.4.0d6) then
2589 write (iout,*) "Energy too large",etot, &
2590 " trying another random conformation"
2593 call gen_rand_conf(itmp,*30)
2595 30 write (iout,*) 'Failed to generate random conformation', &
2597 write (*,*) 'Processor:',me, &
2598 ' Failed to generate random conformation',&
2607 write (iout,'(a,i3,a)') 'Processor:',me, &
2608 ' error in generating random conformation.'
2609 write (*,'(a,i3,a)') 'Processor:',me, &
2610 ' error in generating random conformation.'
2613 ! call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
2614 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
2624 write (iout,'(a,i3,a)') 'Processor:',me, &
2625 ' failed to generate a low-energy random conformation.'
2626 write (*,'(a,i3,a,f10.3)') 'Processor:',me, &
2627 ' failed to generate a low-energy random conformation.',etot
2631 ! call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
2632 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
2639 call chainbuild_cart
2644 kinetic_T=2.0d0/(dimen3*Rb)*EK
2645 if(me.eq.king.or..not.out1file)then
2655 write(iout,*) "before ETOTAL"
2656 call etotal(potEcomp)
2657 if (large) call enerprint(potEcomp)
2660 t_etotal=t_etotal+MPI_Wtime()-tt0
2662 t_etotal=t_etotal+tcpu()-tt0
2669 if (amax*d_time .gt. dvmax) then
2670 d_time=d_time*dvmax/amax
2671 if(me.eq.king.or..not.out1file) write (iout,*) &
2672 "Time step reduced to",d_time,&
2673 " because of too large initial acceleration."
2675 if(me.eq.king.or..not.out1file)then
2676 write(iout,*) "Potential energy and its components"
2677 call enerprint(potEcomp)
2678 ! write(iout,*) (potEcomp(i),i=0,n_ene)
2680 potE=potEcomp(0)-potEcomp(20)
2684 if (ntwe.ne.0) call statout(itime)
2685 if(me.eq.king.or..not.out1file) &
2686 write (iout,'(/a/3(a25,1pe14.5/))') "Initial:", &
2687 " Kinetic energy",EK," Potential energy",potE, &
2688 " Total energy",totE," Maximum acceleration ", &
2691 write (iout,*) "Initial coordinates"
2693 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),&
2696 write (iout,*) "Initial dC"
2698 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),&
2699 (dc(j,i+nres),j=1,3)
2701 write (iout,*) "Initial velocities"
2702 write (iout,"(13x,' backbone ',23x,' side chain')")
2704 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2705 (d_t(j,i+nres),j=1,3)
2707 write (iout,*) "Initial accelerations"
2709 ! write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2710 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),&
2711 (d_a(j,i+nres),j=1,3)
2717 d_t_old(j,i)=d_t(j,i)
2718 d_a_old(j,i)=d_a(j,i)
2720 ! write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
2729 call etotal_short(energia_short)
2730 if (large) call enerprint(potEcomp)
2733 t_eshort=t_eshort+MPI_Wtime()-tt0
2735 t_eshort=t_eshort+tcpu()-tt0
2740 if(.not.out1file .and. large) then
2741 write (iout,*) "energia_long",energia_long(0),&
2742 " energia_short",energia_short(0),&
2743 " total",energia_long(0)+energia_short(0)
2744 write (iout,*) "Initial fast-force accelerations"
2746 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
2747 (d_a(j,i+nres),j=1,3)
2750 ! 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
2753 d_a_short(j,i)=d_a(j,i)
2762 call etotal_long(energia_long)
2763 if (large) call enerprint(potEcomp)
2766 t_elong=t_elong+MPI_Wtime()-tt0
2768 t_elong=t_elong+tcpu()-tt0
2773 if(.not.out1file .and. large) then
2774 write (iout,*) "energia_long",energia_long(0)
2775 write (iout,*) "Initial slow-force accelerations"
2777 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
2778 (d_a(j,i+nres),j=1,3)
2782 t_enegrad=t_enegrad+MPI_Wtime()-tt0
2784 t_enegrad=t_enegrad+tcpu()-tt0
2788 end subroutine init_MD
2789 !-----------------------------------------------------------------------------
2790 subroutine random_vel
2792 ! implicit real*8 (a-h,o-z)
2794 use random, only:anorm_distr
2796 ! include 'DIMENSIONS'
2797 ! include 'COMMON.CONTROL'
2798 ! include 'COMMON.VAR'
2799 ! include 'COMMON.MD'
2801 ! include 'COMMON.LANGEVIN'
2803 ! include 'COMMON.LANGEVIN.lang0'
2805 ! include 'COMMON.CHAIN'
2806 ! include 'COMMON.DERIV'
2807 ! include 'COMMON.GEO'
2808 ! include 'COMMON.LOCAL'
2809 ! include 'COMMON.INTERACT'
2810 ! include 'COMMON.IOUNITS'
2811 ! include 'COMMON.NAMES'
2812 ! include 'COMMON.TIME1'
2813 real(kind=8) :: xv,sigv,lowb,highb ,Ek1
2816 real(kind=8) ,allocatable, dimension(:) :: DDU1,DDU2,DL2,DL1,xsolv,DML,rs
2817 real(kind=8) :: sumx
2819 real(kind=8) ,allocatable, dimension(:) :: rsold
2820 real (kind=8),allocatable,dimension(:,:) :: matold
2824 integer :: i,j,ii,k,ind,mark,imark,mnum
2825 ! Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
2826 ! First generate velocities in the eigenspace of the G matrix
2827 ! write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
2834 sigv=dsqrt((Rb*t_bath)/geigen(i))
2837 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
2839 write (iout,*) "i",i," ii",ii," geigen",geigen(i),&
2840 " d_t_work_new",d_t_work_new(ii)
2851 Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
2854 write (iout,*) "Ek from eigenvectors",Ek1
2855 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
2859 ! Transform velocities to UNRES coordinate space
2860 allocate (DL1(2*nres))
2861 allocate (DDU1(2*nres))
2862 allocate (DL2(2*nres))
2863 allocate (DDU2(2*nres))
2864 allocate (xsolv(2*nres))
2865 allocate (DML(2*nres))
2866 allocate (rs(2*nres))
2868 allocate (rsold(2*nres))
2869 allocate (matold(2*nres,2*nres))
2871 matold(1,1)=DMorig(1)
2872 matold(1,2)=DU1orig(1)
2873 matold(1,3)=DU2orig(1)
2874 write (*,*) DMorig(1),DU1orig(1),DU2orig(1)
2879 matold(i,j)=DMorig(i)
2880 matold(i,j-1)=DU1orig(i-1)
2882 matold(i,j-2)=DU2orig(i-2)
2890 matold(i,j+1)=DU1orig(i)
2896 matold(i,j+2)=DU2orig(i)
2900 matold(dimen,dimen)=DMorig(dimen)
2901 matold(dimen,dimen-1)=DU1orig(dimen-1)
2902 matold(dimen,dimen-2)=DU2orig(dimen-2)
2903 write(iout,*) "old gmatrix"
2904 call matout(dimen,dimen,2*nres,2*nres,matold)
2908 ! Find the ith eigenvector of the pentadiagonal inertiq matrix
2912 DML(j)=DMorig(j)-geigen(i)
2915 DML(j-1)=DMorig(j)-geigen(i)
2920 DDU1(imark-1)=DU2orig(imark-1)
2921 do j=imark+1,dimen-1
2922 DDU1(j-1)=DU1orig(j)
2930 DDU2(j)=DU2orig(j+1)
2939 write (iout,*) "DL2,DL1,DML,DDU1,DDU2"
2940 write(iout,'(10f10.5)') (DL2(k),k=3,dimen-1)
2941 write(iout,'(10f10.5)') (DL1(k),k=2,dimen-1)
2942 write(iout,'(10f10.5)') (DML(k),k=1,dimen-1)
2943 write(iout,'(10f10.5)') (DDU1(k),k=1,dimen-2)
2944 write(iout,'(10f10.5)') (DDU2(k),k=1,dimen-3)
2947 if (imark.gt.2) rs(imark-2)=-DU2orig(imark-2)
2948 if (imark.gt.1) rs(imark-1)=-DU1orig(imark-1)
2949 if (imark.lt.dimen) rs(imark)=-DU1orig(imark)
2950 if (imark.lt.dimen-1) rs(imark+1)=-DU2orig(imark)
2954 ! write (iout,*) "Vector rs"
2956 ! write (iout,*) j,rs(j)
2959 call FDIAG(dimen-1,DL2,DL1,DML,DDU1,DDU2,rs,xsolv,mark)
2966 sumx=-geigen(i)*xsolv(j)
2968 sumx=sumx+matold(j,k)*xsolv(k)
2971 sumx=sumx+matold(j,k)*xsolv(k-1)
2973 write(iout,'(i5,3f10.5)') j,sumx,rsold(j),sumx-rsold(j)
2976 sumx=-geigen(i)*xsolv(j-1)
2978 sumx=sumx+matold(j,k)*xsolv(k)
2981 sumx=sumx+matold(j,k)*xsolv(k-1)
2983 write(iout,'(i5,3f10.5)') j-1,sumx,rsold(j-1),sumx-rsold(j-1)
2987 "Solution of equations system",i," for eigenvalue",geigen(i)
2989 write(iout,'(i5,f10.5)') j,xsolv(j)
2992 do j=dimen-1,imark,-1
2997 write (iout,*) "Un-normalized eigenvector",i," for eigenvalue",geigen(i)
2999 write(iout,'(i5,f10.5)') j,xsolv(j)
3002 ! Normalize ith eigenvector
3005 sumx=sumx+xsolv(j)**2
3009 xsolv(j)=xsolv(j)/sumx
3012 write (iout,*) "Eigenvector",i," for eigenvalue",geigen(i)
3014 write(iout,'(i5,3f10.5)') j,xsolv(j)
3017 ! All done at this point for eigenvector i, exit loop
3025 write (iout,*) "Unable to find eigenvector",i
3028 ! write (iout,*) "i=",i
3030 ! write (iout,*) "k=",k
3033 ! write(iout,*) "ind",ind," ind1",3*(i-1)+k
3034 d_t_work(ind)=d_t_work(ind) &
3035 +xsolv(j)*d_t_work_new(3*(i-1)+k)
3038 enddo ! i (loop over the eigenvectors)
3041 write (iout,*) "d_t_work"
3043 write (iout,"(i5,f10.5)") i,d_t_work(i)
3048 ! if (itype(i,1).eq.10) then
3050 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
3051 iti=iabs(itype(i,mnum))
3052 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3053 if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
3054 .or.(mnum.eq.5)) then
3061 ! write (iout,*) "k",k," ii+k",ii+k," ii+j+k",ii+j+k,"EK1",Ek1
3062 Ek1=Ek1+0.5d0*mp(mnum)*((d_t_work(ii+j+k)+d_t_work_new(ii+k))/2)**2+&
3063 0.5d0*0.25d0*IP(mnum)*(d_t_work(ii+j+k)-d_t_work_new(ii+k))**2
3066 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3067 .and.(mnum.ne.5)) ii=ii+3
3068 write (iout,*) "i",i," itype",itype(i,mnum)," mass",msc(itype(i,mnum),mnum)
3069 write (iout,*) "ii",ii
3072 write (iout,*) "k",k," ii",ii,"EK1",EK1
3073 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3075 Ek1=Ek1+0.5d0*Isc(iabs(itype(i,mnum)),mnum)*(d_t_work(ii)-d_t_work(ii-3))**2
3076 Ek1=Ek1+0.5d0*msc(iabs(itype(i,mnum)),mnum)*d_t_work(ii)**2
3078 write (iout,*) "i",i," ii",ii
3080 write (iout,*) "Ek from d_t_work",Ek1
3081 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
3083 if(allocated(DDU1)) deallocate(DDU1)
3084 if(allocated(DDU2)) deallocate(DDU2)
3085 if(allocated(DL2)) deallocate(DL2)
3086 if(allocated(DL1)) deallocate(DL1)
3087 if(allocated(xsolv)) deallocate(xsolv)
3088 if(allocated(DML)) deallocate(DML)
3089 if(allocated(rs)) deallocate(rs)
3091 if(allocated(matold)) deallocate(matold)
3092 if(allocated(rsold)) deallocate(rsold)
3097 d_t(k,j)=d_t_work(ind)
3101 if (itype(j,1).ne.10 .and. itype(j,mnum).ne.ntyp1_molec(mnum)&
3102 .and.(mnum.ne.5)) then
3104 d_t(k,j+nres)=d_t_work(ind)
3110 write (iout,*) "Random velocities in the Calpha,SC space"
3112 write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3)
3115 write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3)
3122 ! if (itype(i,1).eq.10) then
3124 if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
3125 .or.(mnum.eq.5)) then
3127 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
3131 d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
3132 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
3137 write (iout,*)"Random velocities in the virtual-bond-vector space"
3139 write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3)
3142 write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3)
3145 write (iout,*) "Ek from d_t_work",Ek1
3146 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
3154 d_t_work(ind)=d_t_work(ind) &
3155 +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
3157 ! write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
3161 ! Transfer to the d_t vector
3163 d_t(j,0)=d_t_work(j)
3169 d_t(j,i)=d_t_work(ind)
3174 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3175 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3176 .and.(mnum.ne.5)) then
3179 d_t(j,i+nres)=d_t_work(ind)
3185 ! write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",&
3186 ! 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
3188 ! write(iout,*) "end init MD"
3190 end subroutine random_vel
3191 !-----------------------------------------------------------------------------
3193 subroutine sd_verlet_p_setup
3194 ! Sets up the parameters of stochastic Verlet algorithm
3195 ! implicit real*8 (a-h,o-z)
3196 ! include 'DIMENSIONS'
3197 use control, only: tcpu
3202 ! include 'COMMON.CONTROL'
3203 ! include 'COMMON.VAR'
3204 ! include 'COMMON.MD'
3206 ! include 'COMMON.LANGEVIN'
3208 ! include 'COMMON.LANGEVIN.lang0'
3210 ! include 'COMMON.CHAIN'
3211 ! include 'COMMON.DERIV'
3212 ! include 'COMMON.GEO'
3213 ! include 'COMMON.LOCAL'
3214 ! include 'COMMON.INTERACT'
3215 ! include 'COMMON.IOUNITS'
3216 ! include 'COMMON.NAMES'
3217 ! include 'COMMON.TIME1'
3218 real(kind=8),dimension(6*nres) :: emgdt !(MAXRES6) maxres6=6*maxres
3219 real(kind=8) :: pterm,vterm,rho,rhoc,vsig
3220 real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
3221 prand_vec,vrand_vec1,vrand_vec2 !(MAXRES6) maxres6=6*maxres
3222 logical :: lprn = .false.
3223 real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
3224 real(kind=8) :: ktm,gdt,egdt,gdt2,gdt3,gdt4,gdt5,gdt6,gdt7,gdt8,&
3226 integer :: i,maxres2
3233 ! AL 8/17/04 Code adapted from tinker
3235 ! Get the frictional and random terms for stochastic dynamics in the
3236 ! eigenspace of mass-scaled UNRES friction matrix
3240 gdt = fricgam(i) * d_time
3242 ! Stochastic dynamics reduces to simple MD for zero friction
3244 if (gdt .le. zero) then
3245 pfric_vec(i) = 1.0d0
3246 vfric_vec(i) = d_time
3247 afric_vec(i) = 0.5d0 * d_time * d_time
3248 prand_vec(i) = 0.0d0
3249 vrand_vec1(i) = 0.0d0
3250 vrand_vec2(i) = 0.0d0
3252 ! Analytical expressions when friction coefficient is large
3255 if (gdt .ge. gdt_radius) then
3258 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
3259 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
3260 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
3261 vterm = 1.0d0 - egdt**2
3262 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
3264 ! Use series expansions when friction coefficient is small
3275 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0 &
3276 - gdt5/120.0d0 + gdt6/720.0d0 &
3277 - gdt7/5040.0d0 + gdt8/40320.0d0 &
3278 - gdt9/362880.0d0) / fricgam(i)**2
3279 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
3280 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
3281 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0 &
3282 + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0 &
3283 + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0 &
3284 + 127.0d0*gdt9/90720.0d0
3285 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0 &
3286 - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0 &
3287 - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0 &
3288 - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
3289 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0 &
3290 - 17.0d0*gdt2/1280.0d0 &
3291 + 17.0d0*gdt3/6144.0d0 &
3292 + 40967.0d0*gdt4/34406400.0d0 &
3293 - 57203.0d0*gdt5/275251200.0d0 &
3294 - 1429487.0d0*gdt6/13212057600.0d0)
3297 ! Compute the scaling factors of random terms for the nonzero friction case
3299 ktm = 0.5d0*d_time/fricgam(i)
3300 psig = dsqrt(ktm*pterm) / fricgam(i)
3301 vsig = dsqrt(ktm*vterm)
3302 rhoc = dsqrt(1.0d0 - rho*rho)
3304 vrand_vec1(i) = vsig * rho
3305 vrand_vec2(i) = vsig * rhoc
3310 "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",&
3313 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),&
3314 afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
3318 ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
3321 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
3322 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
3323 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
3324 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
3325 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
3326 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
3329 t_sdsetup=t_sdsetup+MPI_Wtime()
3331 t_sdsetup=t_sdsetup+tcpu()-tt0
3334 end subroutine sd_verlet_p_setup
3335 !-----------------------------------------------------------------------------
3336 subroutine eigtransf1(n,ndim,ab,d,c)
3340 real(kind=8) :: ab(ndim,ndim,n),c(ndim,n),d(ndim)
3346 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
3351 end subroutine eigtransf1
3352 !-----------------------------------------------------------------------------
3353 subroutine eigtransf(n,ndim,a,b,d,c)
3357 real(kind=8) :: a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
3363 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
3368 end subroutine eigtransf
3369 !-----------------------------------------------------------------------------
3370 subroutine sd_verlet1
3372 ! Applying stochastic velocity Verlet algorithm - step 1 to velocities
3374 ! implicit real*8 (a-h,o-z)
3375 ! include 'DIMENSIONS'
3376 ! include 'COMMON.CONTROL'
3377 ! include 'COMMON.VAR'
3378 ! include 'COMMON.MD'
3380 ! include 'COMMON.LANGEVIN'
3382 ! include 'COMMON.LANGEVIN.lang0'
3384 ! include 'COMMON.CHAIN'
3385 ! include 'COMMON.DERIV'
3386 ! include 'COMMON.GEO'
3387 ! include 'COMMON.LOCAL'
3388 ! include 'COMMON.INTERACT'
3389 ! include 'COMMON.IOUNITS'
3390 ! include 'COMMON.NAMES'
3391 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
3392 !el common /stochcalc/ stochforcvec
3393 logical :: lprn = .false.
3394 real(kind=8) :: ddt1,ddt2
3395 integer :: i,j,ind,inres
3397 ! write (iout,*) "dc_old"
3399 ! write (iout,'(i5,3f10.5,5x,3f10.5)')
3400 ! & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3403 dc_work(j)=dc_old(j,0)
3404 d_t_work(j)=d_t_old(j,0)
3405 d_a_work(j)=d_a_old(j,0)
3410 dc_work(ind+j)=dc_old(j,i)
3411 d_t_work(ind+j)=d_t_old(j,i)
3412 d_a_work(ind+j)=d_a_old(j,i)
3418 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3419 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3420 .and.(mnum.ne.5)) then
3422 dc_work(ind+j)=dc_old(j,i+nres)
3423 d_t_work(ind+j)=d_t_old(j,i+nres)
3424 d_a_work(ind+j)=d_a_old(j,i+nres)
3432 "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",&
3436 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
3437 vfric_mat(i,j),afric_mat(i,j),&
3438 prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3446 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) &
3447 +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3448 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3449 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3451 d_t_work_new(i)=ddt1+0.5d0*ddt2
3452 d_t_work(i)=ddt1+ddt2
3457 d_t(j,0)=d_t_work(j)
3462 dc(j,i)=dc_work(ind+j)
3463 d_t(j,i)=d_t_work(ind+j)
3469 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3470 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3471 .and.(mnum.ne.5)) then
3474 dc(j,inres)=dc_work(ind+j)
3475 d_t(j,inres)=d_t_work(ind+j)
3481 end subroutine sd_verlet1
3482 !-----------------------------------------------------------------------------
3483 subroutine sd_verlet2
3485 ! Calculating the adjusted velocities for accelerations
3487 ! implicit real*8 (a-h,o-z)
3488 ! include 'DIMENSIONS'
3489 ! include 'COMMON.CONTROL'
3490 ! include 'COMMON.VAR'
3491 ! include 'COMMON.MD'
3493 ! include 'COMMON.LANGEVIN'
3495 ! include 'COMMON.LANGEVIN.lang0'
3497 ! include 'COMMON.CHAIN'
3498 ! include 'COMMON.DERIV'
3499 ! include 'COMMON.GEO'
3500 ! include 'COMMON.LOCAL'
3501 ! include 'COMMON.INTERACT'
3502 ! include 'COMMON.IOUNITS'
3503 ! include 'COMMON.NAMES'
3504 !el real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV !(MAXRES6) maxres6=6*maxres
3505 real(kind=8),dimension(6*nres) :: stochforcvecV !(MAXRES6) maxres6=6*maxres
3506 !el common /stochcalc/ stochforcvec
3508 real(kind=8) :: ddt1,ddt2
3509 integer :: i,j,ind,inres
3510 ! Compute the stochastic forces which contribute to velocity change
3512 call stochastic_force(stochforcvecV)
3519 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3520 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+ &
3521 vrand_mat2(i,j)*stochforcvecV(j)
3523 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
3527 d_t(j,0)=d_t_work(j)
3532 d_t(j,i)=d_t_work(ind+j)
3538 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3539 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3540 .and.(mnum.ne.5)) then
3543 d_t(j,inres)=d_t_work(ind+j)
3549 end subroutine sd_verlet2
3550 !-----------------------------------------------------------------------------
3551 subroutine sd_verlet_ciccotti_setup
3553 ! Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
3555 ! implicit real*8 (a-h,o-z)
3556 ! include 'DIMENSIONS'
3557 use control, only: tcpu
3562 ! include 'COMMON.CONTROL'
3563 ! include 'COMMON.VAR'
3564 ! include 'COMMON.MD'
3566 ! include 'COMMON.LANGEVIN'
3568 ! include 'COMMON.LANGEVIN.lang0'
3570 ! include 'COMMON.CHAIN'
3571 ! include 'COMMON.DERIV'
3572 ! include 'COMMON.GEO'
3573 ! include 'COMMON.LOCAL'
3574 ! include 'COMMON.INTERACT'
3575 ! include 'COMMON.IOUNITS'
3576 ! include 'COMMON.NAMES'
3577 ! include 'COMMON.TIME1'
3578 real(kind=8),dimension(6*nres) :: emgdt !(MAXRES6) maxres6=6*maxres
3579 real(kind=8) :: pterm,vterm,rho,rhoc,vsig
3580 real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
3581 prand_vec,vrand_vec1,vrand_vec2 !(MAXRES6) maxres6=6*maxres
3582 logical :: lprn = .false.
3583 real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
3584 real(kind=8) :: ktm,gdt,egdt,tt0
3585 integer :: i,maxres2
3592 ! AL 8/17/04 Code adapted from tinker
3594 ! Get the frictional and random terms for stochastic dynamics in the
3595 ! eigenspace of mass-scaled UNRES friction matrix
3599 write (iout,*) "i",i," fricgam",fricgam(i)
3600 gdt = fricgam(i) * d_time
3602 ! Stochastic dynamics reduces to simple MD for zero friction
3604 if (gdt .le. zero) then
3605 pfric_vec(i) = 1.0d0
3606 vfric_vec(i) = d_time
3607 afric_vec(i) = 0.5d0*d_time*d_time
3608 prand_vec(i) = afric_vec(i)
3609 vrand_vec2(i) = vfric_vec(i)
3611 ! Analytical expressions when friction coefficient is large
3616 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
3617 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
3618 prand_vec(i) = afric_vec(i)
3619 vrand_vec2(i) = vfric_vec(i)
3621 ! Compute the scaling factors of random terms for the nonzero friction case
3623 ! ktm = 0.5d0*d_time/fricgam(i)
3624 ! psig = dsqrt(ktm*pterm) / fricgam(i)
3625 ! vsig = dsqrt(ktm*vterm)
3626 ! prand_vec(i) = psig*afric_vec(i)
3627 ! vrand_vec2(i) = vsig*vfric_vec(i)
3632 "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",&
3635 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),&
3636 afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
3640 ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
3642 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
3643 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
3644 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
3645 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
3646 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
3648 t_sdsetup=t_sdsetup+MPI_Wtime()
3650 t_sdsetup=t_sdsetup+tcpu()-tt0
3653 end subroutine sd_verlet_ciccotti_setup
3654 !-----------------------------------------------------------------------------
3655 subroutine sd_verlet1_ciccotti
3657 ! Applying stochastic velocity Verlet algorithm - step 1 to velocities
3658 ! implicit real*8 (a-h,o-z)
3660 ! include 'DIMENSIONS'
3664 ! include 'COMMON.CONTROL'
3665 ! include 'COMMON.VAR'
3666 ! include 'COMMON.MD'
3668 ! include 'COMMON.LANGEVIN'
3670 ! include 'COMMON.LANGEVIN.lang0'
3672 ! include 'COMMON.CHAIN'
3673 ! include 'COMMON.DERIV'
3674 ! include 'COMMON.GEO'
3675 ! include 'COMMON.LOCAL'
3676 ! include 'COMMON.INTERACT'
3677 ! include 'COMMON.IOUNITS'
3678 ! include 'COMMON.NAMES'
3679 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
3680 !el common /stochcalc/ stochforcvec
3681 logical :: lprn = .false.
3682 real(kind=8) :: ddt1,ddt2
3683 integer :: i,j,ind,inres
3684 ! write (iout,*) "dc_old"
3686 ! write (iout,'(i5,3f10.5,5x,3f10.5)')
3687 ! & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3690 dc_work(j)=dc_old(j,0)
3691 d_t_work(j)=d_t_old(j,0)
3692 d_a_work(j)=d_a_old(j,0)
3697 dc_work(ind+j)=dc_old(j,i)
3698 d_t_work(ind+j)=d_t_old(j,i)
3699 d_a_work(ind+j)=d_a_old(j,i)
3704 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3706 dc_work(ind+j)=dc_old(j,i+nres)
3707 d_t_work(ind+j)=d_t_old(j,i+nres)
3708 d_a_work(ind+j)=d_a_old(j,i+nres)
3717 "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",&
3721 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
3722 vfric_mat(i,j),afric_mat(i,j),&
3723 prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3731 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) &
3732 +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3733 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3734 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3736 d_t_work_new(i)=ddt1+0.5d0*ddt2
3737 d_t_work(i)=ddt1+ddt2
3742 d_t(j,0)=d_t_work(j)
3747 dc(j,i)=dc_work(ind+j)
3748 d_t(j,i)=d_t_work(ind+j)
3754 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3755 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3756 .and.(mnum.ne.5)) then
3759 dc(j,inres)=dc_work(ind+j)
3760 d_t(j,inres)=d_t_work(ind+j)
3766 end subroutine sd_verlet1_ciccotti
3767 !-----------------------------------------------------------------------------
3768 subroutine sd_verlet2_ciccotti
3770 ! Calculating the adjusted velocities for accelerations
3772 ! implicit real*8 (a-h,o-z)
3773 ! include 'DIMENSIONS'
3774 ! include 'COMMON.CONTROL'
3775 ! include 'COMMON.VAR'
3776 ! include 'COMMON.MD'
3778 ! include 'COMMON.LANGEVIN'
3780 ! include 'COMMON.LANGEVIN.lang0'
3782 ! include 'COMMON.CHAIN'
3783 ! include 'COMMON.DERIV'
3784 ! include 'COMMON.GEO'
3785 ! include 'COMMON.LOCAL'
3786 ! include 'COMMON.INTERACT'
3787 ! include 'COMMON.IOUNITS'
3788 ! include 'COMMON.NAMES'
3789 !el real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV !(MAXRES6) maxres6=6*maxres
3790 real(kind=8),dimension(6*nres) :: stochforcvecV !(MAXRES6) maxres6=6*maxres
3791 !el common /stochcalc/ stochforcvec
3792 real(kind=8) :: ddt1,ddt2
3793 integer :: i,j,ind,inres
3795 ! Compute the stochastic forces which contribute to velocity change
3797 call stochastic_force(stochforcvecV)
3804 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3805 ! ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
3806 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
3808 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
3812 d_t(j,0)=d_t_work(j)
3817 d_t(j,i)=d_t_work(ind+j)
3823 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3825 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3828 d_t(j,inres)=d_t_work(ind+j)
3834 end subroutine sd_verlet2_ciccotti
3836 !-----------------------------------------------------------------------------
3838 !-----------------------------------------------------------------------------
3839 subroutine inertia_tensor
3841 ! Calculating the intertia tensor for the entire protein in order to
3842 ! remove the perpendicular components of velocity matrix which cause
3843 ! the molecule to rotate.
3846 ! implicit real*8 (a-h,o-z)
3847 ! include 'DIMENSIONS'
3848 ! include 'COMMON.CONTROL'
3849 ! include 'COMMON.VAR'
3850 ! include 'COMMON.MD'
3851 ! include 'COMMON.CHAIN'
3852 ! include 'COMMON.DERIV'
3853 ! include 'COMMON.GEO'
3854 ! include 'COMMON.LOCAL'
3855 ! include 'COMMON.INTERACT'
3856 ! include 'COMMON.IOUNITS'
3857 ! include 'COMMON.NAMES'
3859 real(kind=8),dimension(3,3) :: Im,Imcp,eigvec,Id
3860 real(kind=8),dimension(3) :: pr,eigval,L,vp,vrot
3861 real(kind=8) :: M_SC,mag,mag2,M_PEP
3862 real(kind=8),dimension(3,0:nres) :: vpp !(3,0:MAXRES)
3863 real(kind=8),dimension(3) :: vs_p,pp,incr,v
3864 real(kind=8),dimension(3,3) :: pr1,pr2
3866 !el common /gucio/ cm
3867 integer :: iti,inres,i,j,k,mnum
3878 ! calculating the center of the mass of the protein
3882 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
3883 if (itype(i,mnum).eq.ntyp1_molec(mnum)) cycle
3884 M_PEP=M_PEP+mp(mnum)
3886 cm(j)=cm(j)+(c(j,i)+0.5d0*dc(j,i))*mp(mnum)
3895 if (mnum.eq.5) cycle
3896 iti=iabs(itype(i,mnum))
3897 M_SC=M_SC+msc(iabs(iti),mnum)
3900 cm(j)=cm(j)+msc(iabs(iti),mnum)*c(j,inres)
3904 cm(j)=cm(j)/(M_SC+M_PEP)
3909 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
3911 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
3913 Im(1,1)=Im(1,1)+mp(mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
3914 Im(1,2)=Im(1,2)-mp(mnum)*pr(1)*pr(2)
3915 Im(1,3)=Im(1,3)-mp(mnum)*pr(1)*pr(3)
3916 Im(2,3)=Im(2,3)-mp(mnum)*pr(2)*pr(3)
3917 Im(2,2)=Im(2,2)+mp(mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
3918 Im(3,3)=Im(3,3)+mp(mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
3923 iti=iabs(itype(i,mnum))
3924 if (mnum.eq.5) cycle
3927 pr(j)=c(j,inres)-cm(j)
3929 Im(1,1)=Im(1,1)+msc(iabs(iti),mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
3930 Im(1,2)=Im(1,2)-msc(iabs(iti),mnum)*pr(1)*pr(2)
3931 Im(1,3)=Im(1,3)-msc(iabs(iti),mnum)*pr(1)*pr(3)
3932 Im(2,3)=Im(2,3)-msc(iabs(iti),mnum)*pr(2)*pr(3)
3933 Im(2,2)=Im(2,2)+msc(iabs(iti),mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
3934 Im(3,3)=Im(3,3)+msc(iabs(iti),mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
3939 Im(1,1)=Im(1,1)+Ip(mnum)*(1-dc_norm(1,i)*dc_norm(1,i))* &
3940 vbld(i+1)*vbld(i+1)*0.25d0
3941 Im(1,2)=Im(1,2)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(2,i))* &
3942 vbld(i+1)*vbld(i+1)*0.25d0
3943 Im(1,3)=Im(1,3)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(3,i))* &
3944 vbld(i+1)*vbld(i+1)*0.25d0
3945 Im(2,3)=Im(2,3)+Ip(mnum)*(-dc_norm(2,i)*dc_norm(3,i))* &
3946 vbld(i+1)*vbld(i+1)*0.25d0
3947 Im(2,2)=Im(2,2)+Ip(mnum)*(1-dc_norm(2,i)*dc_norm(2,i))* &
3948 vbld(i+1)*vbld(i+1)*0.25d0
3949 Im(3,3)=Im(3,3)+Ip(mnum)*(1-dc_norm(3,i)*dc_norm(3,i))* &
3950 vbld(i+1)*vbld(i+1)*0.25d0
3956 ! if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
3957 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3958 .and.(mnum.ne.5)) then
3959 iti=iabs(itype(i,mnum))
3961 Im(1,1)=Im(1,1)+Isc(iti,mnum)*(1-dc_norm(1,inres)* &
3962 dc_norm(1,inres))*vbld(inres)*vbld(inres)
3963 Im(1,2)=Im(1,2)-Isc(iti,mnum)*(dc_norm(1,inres)* &
3964 dc_norm(2,inres))*vbld(inres)*vbld(inres)
3965 Im(1,3)=Im(1,3)-Isc(iti,mnum)*(dc_norm(1,inres)* &
3966 dc_norm(3,inres))*vbld(inres)*vbld(inres)
3967 Im(2,3)=Im(2,3)-Isc(iti,mnum)*(dc_norm(2,inres)* &
3968 dc_norm(3,inres))*vbld(inres)*vbld(inres)
3969 Im(2,2)=Im(2,2)+Isc(iti,mnum)*(1-dc_norm(2,inres)* &
3970 dc_norm(2,inres))*vbld(inres)*vbld(inres)
3971 Im(3,3)=Im(3,3)+Isc(iti,mnum)*(1-dc_norm(3,inres)* &
3972 dc_norm(3,inres))*vbld(inres)*vbld(inres)
3977 ! write(iout,*) "The angular momentum before adjustment:"
3978 ! write(iout,*) (L(j),j=1,3)
3984 ! Copying the Im matrix for the djacob subroutine
3992 ! Finding the eigenvectors and eignvalues of the inertia tensor
3993 call djacob(3,3,10000,1.0d-10,Imcp,eigvec,eigval)
3994 ! write (iout,*) "Eigenvalues & Eigenvectors"
3995 ! write (iout,'(5x,3f10.5)') (eigval(i),i=1,3)
3998 ! write (iout,'(i5,3f10.5)') i,(eigvec(i,j),j=1,3)
4000 ! Constructing the diagonalized matrix
4002 if (dabs(eigval(i)).gt.1.0d-15) then
4003 Id(i,i)=1.0d0/eigval(i)
4010 Imcp(i,j)=eigvec(j,i)
4016 pr1(i,j)=pr1(i,j)+Id(i,k)*Imcp(k,j)
4023 pr2(i,j)=pr2(i,j)+eigvec(i,k)*pr1(k,j)
4027 ! Calculating the total rotational velocity of the molecule
4030 vrot(i)=vrot(i)+pr2(i,j)*L(j)
4033 ! Resetting the velocities
4035 call vecpr(vrot(1),dc(1,i),vp)
4037 d_t(j,i)=d_t(j,i)-vp(j)
4042 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4043 .and.(mnum.ne.5)) then
4044 ! if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
4045 ! if(itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
4047 call vecpr(vrot(1),dc(1,inres),vp)
4049 d_t(j,inres)=d_t(j,inres)-vp(j)
4054 ! write(iout,*) "The angular momentum after adjustment:"
4055 ! write(iout,*) (L(j),j=1,3)
4058 end subroutine inertia_tensor
4059 !-----------------------------------------------------------------------------
4060 subroutine angmom(cm,L)
4063 ! implicit real*8 (a-h,o-z)
4064 ! include 'DIMENSIONS'
4065 ! include 'COMMON.CONTROL'
4066 ! include 'COMMON.VAR'
4067 ! include 'COMMON.MD'
4068 ! include 'COMMON.CHAIN'
4069 ! include 'COMMON.DERIV'
4070 ! include 'COMMON.GEO'
4071 ! include 'COMMON.LOCAL'
4072 ! include 'COMMON.INTERACT'
4073 ! include 'COMMON.IOUNITS'
4074 ! include 'COMMON.NAMES'
4075 real(kind=8) :: mscab
4076 real(kind=8),dimension(3) :: L,cm,pr,vp,vrot,incr,v,pp
4077 integer :: iti,inres,i,j,mnum
4078 ! Calculate the angular momentum
4087 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
4089 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
4092 v(j)=incr(j)+0.5d0*d_t(j,i)
4095 incr(j)=incr(j)+d_t(j,i)
4097 call vecpr(pr(1),v(1),vp)
4099 L(j)=L(j)+mp(mnum)*vp(j)
4103 pp(j)=0.5d0*d_t(j,i)
4105 call vecpr(pr(1),pp(1),vp)
4107 L(j)=L(j)+Ip(mnum)*vp(j)
4115 iti=iabs(itype(i,mnum))
4123 pr(j)=c(j,inres)-cm(j)
4125 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4126 .and.(mnum.ne.5)) then
4128 v(j)=incr(j)+d_t(j,inres)
4135 call vecpr(pr(1),v(1),vp)
4136 ! write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),&
4137 ! " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
4139 L(j)=L(j)+mscab*vp(j)
4141 ! write (iout,*) "L",(l(j),j=1,3)
4142 if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4143 .and.(mnum.ne.5)) then
4145 v(j)=incr(j)+d_t(j,inres)
4147 call vecpr(dc(1,inres),d_t(1,inres),vp)
4149 L(j)=L(j)+Isc(iti,mnum)*vp(j)
4153 incr(j)=incr(j)+d_t(j,i)
4157 end subroutine angmom
4158 !-----------------------------------------------------------------------------
4159 subroutine vcm_vel(vcm)
4162 ! implicit real*8 (a-h,o-z)
4163 ! include 'DIMENSIONS'
4164 ! include 'COMMON.VAR'
4165 ! include 'COMMON.MD'
4166 ! include 'COMMON.CHAIN'
4167 ! include 'COMMON.DERIV'
4168 ! include 'COMMON.GEO'
4169 ! include 'COMMON.LOCAL'
4170 ! include 'COMMON.INTERACT'
4171 ! include 'COMMON.IOUNITS'
4172 real(kind=8),dimension(3) :: vcm,vv
4173 real(kind=8) :: summas,amas
4183 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
4185 summas=summas+mp(mnum)
4187 vcm(j)=vcm(j)+mp(mnum)*(vv(j)+0.5d0*d_t(j,i))
4191 amas=msc(iabs(itype(i,mnum)),mnum)
4196 if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4197 .and.(mnum.ne.5)) then
4199 vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres))
4203 vcm(j)=vcm(j)+amas*vv(j)
4207 vv(j)=vv(j)+d_t(j,i)
4210 ! write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas
4212 vcm(j)=vcm(j)/summas
4215 end subroutine vcm_vel
4216 !-----------------------------------------------------------------------------
4218 !-----------------------------------------------------------------------------
4220 ! RATTLE algorithm for velocity Verlet - step 1, UNRES
4224 ! implicit real*8 (a-h,o-z)
4225 ! include 'DIMENSIONS'
4227 ! include 'COMMON.CONTROL'
4228 ! include 'COMMON.VAR'
4229 ! include 'COMMON.MD'
4231 ! include 'COMMON.LANGEVIN'
4233 ! include 'COMMON.LANGEVIN.lang0'
4235 ! include 'COMMON.CHAIN'
4236 ! include 'COMMON.DERIV'
4237 ! include 'COMMON.GEO'
4238 ! include 'COMMON.LOCAL'
4239 ! include 'COMMON.INTERACT'
4240 ! include 'COMMON.IOUNITS'
4241 ! include 'COMMON.NAMES'
4242 ! include 'COMMON.TIME1'
4243 !el real(kind=8) :: gginv(2*nres,2*nres),&
4244 !el gdc(3,2*nres,2*nres)
4245 real(kind=8) :: dC_uncor(3,2*nres) !,&
4246 !el real(kind=8) :: Cmat(2*nres,2*nres)
4247 real(kind=8) :: x(2*nres),xcorr(3,2*nres) !maxres2=2*maxres
4248 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
4249 !el common /przechowalnia/ nbond
4250 integer :: max_rattle = 5
4251 logical :: lprn = .false., lprn1 = .false., not_done
4252 real(kind=8) :: tol_rattle = 1.0d-5
4254 integer :: ii,i,j,jj,l,ind,ind1,nres2
4257 !el /common/ przechowalnia
4259 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4260 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4261 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4263 if (lprn) write (iout,*) "RATTLE1"
4267 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4268 .and.(mnum.ne.5)) nbond=nbond+1
4270 ! Make a folded form of the Ginv-matrix
4283 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4288 if (itype(k,1).ne.10 .and. itype(k,mnum).ne.ntyp1_molec(mnum)&
4289 .and.(mnum.ne.5)) then
4293 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4301 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4312 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4316 if (itype(k,1).ne.10) then
4320 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4328 write (iout,*) "Matrix GGinv"
4329 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
4335 if (iter.gt.max_rattle) then
4336 write (iout,*) "Error - too many iterations in RATTLE."
4339 ! Calculate the matrix C = GG**(-1) dC_old o dC
4344 dC_uncor(j,ind1)=dC(j,i)
4348 if (itype(i,1).ne.10) then
4351 dC_uncor(j,ind1)=dC(j,i+nres)
4360 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
4364 if (itype(k,1).ne.10) then
4367 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
4372 ! Calculate deviations from standard virtual-bond lengths
4376 x(ind)=vbld(i+1)**2-vbl**2
4379 if (itype(i,1).ne.10) then
4381 x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4385 write (iout,*) "Coordinates and violations"
4387 write(iout,'(i5,3f10.5,5x,e15.5)') &
4388 i,(dC_uncor(j,i),j=1,3),x(i)
4390 write (iout,*) "Velocities and violations"
4394 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4395 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4399 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4400 .and.(mnum.ne.5)) then
4403 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4404 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4405 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4408 ! write (iout,*) "gdc"
4410 ! write (iout,*) "i",i
4412 ! write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4418 if (dabs(x(i)).gt.xmax) then
4422 if (xmax.lt.tol_rattle) then
4426 ! Calculate the matrix of the system of equations
4431 Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
4436 write (iout,*) "Matrix Cmat"
4437 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4439 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4440 ! Add constraint term to positions
4447 xx = xx+x(ii)*gdc(j,ind,ii)
4451 d_t_new(j,i)=d_t_new(j,i)-xx/d_time
4455 if (itype(i,1).ne.10) then
4460 xx = xx+x(ii)*gdc(j,ind,ii)
4463 dC(j,i+nres)=dC(j,i+nres)-xx
4464 d_t_new(j,i+nres)=d_t_new(j,i+nres)-xx/d_time
4468 ! Rebuild the chain using the new coordinates
4469 call chainbuild_cart
4471 write (iout,*) "New coordinates, Lagrange multipliers,",&
4472 " and differences between actual and standard bond lengths"
4476 xx=vbld(i+1)**2-vbl**2
4477 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4478 i,(dC(j,i),j=1,3),x(ind),xx
4482 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4485 xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4486 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4487 i,(dC(j,i+nres),j=1,3),x(ind),xx
4490 write (iout,*) "Velocities and violations"
4494 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4495 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4498 if (itype(i,1).ne.10) then
4500 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4501 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4502 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4509 10 write (iout,*) "Error - singularity in solving the system",&
4510 " of equations for Lagrange multipliers."
4514 "RATTLE inactive; use -DRATTLE switch at compile time."
4517 end subroutine rattle1
4518 !-----------------------------------------------------------------------------
4520 ! RATTLE algorithm for velocity Verlet - step 2, UNRES
4524 ! implicit real*8 (a-h,o-z)
4525 ! include 'DIMENSIONS'
4527 ! include 'COMMON.CONTROL'
4528 ! include 'COMMON.VAR'
4529 ! include 'COMMON.MD'
4531 ! include 'COMMON.LANGEVIN'
4533 ! include 'COMMON.LANGEVIN.lang0'
4535 ! include 'COMMON.CHAIN'
4536 ! include 'COMMON.DERIV'
4537 ! include 'COMMON.GEO'
4538 ! include 'COMMON.LOCAL'
4539 ! include 'COMMON.INTERACT'
4540 ! include 'COMMON.IOUNITS'
4541 ! include 'COMMON.NAMES'
4542 ! include 'COMMON.TIME1'
4543 !el real(kind=8) :: gginv(2*nres,2*nres),&
4544 !el gdc(3,2*nres,2*nres)
4545 real(kind=8) :: dC_uncor(3,2*nres) !,&
4546 !el Cmat(2*nres,2*nres)
4547 real(kind=8) :: x(2*nres) !maxres2=2*maxres
4548 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
4549 !el common /przechowalnia/ nbond
4550 integer :: max_rattle = 5
4551 logical :: lprn = .false., lprn1 = .false., not_done
4552 real(kind=8) :: tol_rattle = 1.0d-5
4556 !el /common/ przechowalnia
4557 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4558 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4559 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4561 if (lprn) write (iout,*) "RATTLE2"
4562 if (lprn) write (iout,*) "Velocity correction"
4563 ! Calculate the matrix G dC
4569 gdc(j,i,ind)=GGinv(i,ind)*dC(j,k)
4574 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4575 .and.(mnum.ne.5)) then
4578 gdc(j,i,ind)=GGinv(i,ind)*dC(j,k+nres)
4584 ! write (iout,*) "gdc"
4586 ! write (iout,*) "i",i
4588 ! write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4592 ! Calculate the matrix of the system of equations
4599 Cmat(ind,j)=Cmat(ind,j)+dC(k,i)*gdc(k,ind,j)
4605 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4606 .and.(mnum.ne.5)) then
4611 Cmat(ind,j)=Cmat(ind,j)+dC(k,i+nres)*gdc(k,ind,j)
4616 ! Calculate the scalar product dC o d_t_new
4620 x(ind)=scalar(d_t(1,i),dC(1,i))
4624 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4625 .and.(mnum.ne.5)) then
4627 x(ind)=scalar(d_t(1,i+nres),dC(1,i+nres))
4631 write (iout,*) "Velocities and violations"
4635 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4636 i,ind,(d_t(j,i),j=1,3),x(ind)
4640 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4641 .and.(mnum.ne.5)) then
4643 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4644 i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind)
4650 if (dabs(x(i)).gt.xmax) then
4654 if (xmax.lt.tol_rattle) then
4659 write (iout,*) "Matrix Cmat"
4660 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4662 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4663 ! Add constraint term to velocities
4670 xx = xx+x(ii)*gdc(j,ind,ii)
4672 d_t(j,i)=d_t(j,i)-xx
4677 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4678 .and.(mnum.ne.5)) then
4683 xx = xx+x(ii)*gdc(j,ind,ii)
4685 d_t(j,i+nres)=d_t(j,i+nres)-xx
4691 "New velocities, Lagrange multipliers violations"
4695 if (lprn) write (iout,'(2i5,3f10.5,5x,2e15.5)') &
4696 i,ind,(d_t(j,i),j=1,3),x(ind),scalar(d_t(1,i),dC(1,i))
4700 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4703 write (iout,'(2i5,3f10.5,5x,2e15.5)') &
4704 i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind),&
4705 scalar(d_t(1,i+nres),dC(1,i+nres))
4711 10 write (iout,*) "Error - singularity in solving the system",&
4712 " of equations for Lagrange multipliers."
4716 "RATTLE inactive; use -DRATTLE option at compile time."
4719 end subroutine rattle2
4720 !-----------------------------------------------------------------------------
4721 subroutine rattle_brown
4722 ! RATTLE/LINCS algorithm for Brownian dynamics, UNRES
4726 ! implicit real*8 (a-h,o-z)
4727 ! include 'DIMENSIONS'
4729 ! include 'COMMON.CONTROL'
4730 ! include 'COMMON.VAR'
4731 ! include 'COMMON.MD'
4733 ! include 'COMMON.LANGEVIN'
4735 ! include 'COMMON.LANGEVIN.lang0'
4737 ! include 'COMMON.CHAIN'
4738 ! include 'COMMON.DERIV'
4739 ! include 'COMMON.GEO'
4740 ! include 'COMMON.LOCAL'
4741 ! include 'COMMON.INTERACT'
4742 ! include 'COMMON.IOUNITS'
4743 ! include 'COMMON.NAMES'
4744 ! include 'COMMON.TIME1'
4745 !el real(kind=8) :: gginv(2*nres,2*nres),&
4746 !el gdc(3,2*nres,2*nres)
4747 real(kind=8) :: dC_uncor(3,2*nres) !,&
4748 !el real(kind=8) :: Cmat(2*nres,2*nres)
4749 real(kind=8) :: x(2*nres) !maxres2=2*maxres
4750 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
4751 !el common /przechowalnia/ nbond
4752 integer :: max_rattle = 5
4753 logical :: lprn = .false., lprn1 = .false., not_done
4754 real(kind=8) :: tol_rattle = 1.0d-5
4758 !el /common/ przechowalnia
4759 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4760 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4761 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4764 if (lprn) write (iout,*) "RATTLE_BROWN"
4767 if (itype(i,1).ne.10) nbond=nbond+1
4769 ! Make a folded form of the Ginv-matrix
4782 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4786 if (itype(k,1).ne.10) then
4790 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4797 if (itype(i,1).ne.10) then
4807 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4811 if (itype(k,1).ne.10) then
4815 if (j.eq.1 .and. l.eq.1)GGinv(ii,jj)=fricmat(ind,ind1)
4823 write (iout,*) "Matrix GGinv"
4824 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
4830 if (iter.gt.max_rattle) then
4831 write (iout,*) "Error - too many iterations in RATTLE."
4834 ! Calculate the matrix C = GG**(-1) dC_old o dC
4839 dC_uncor(j,ind1)=dC(j,i)
4843 if (itype(i,1).ne.10) then
4846 dC_uncor(j,ind1)=dC(j,i+nres)
4855 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
4859 if (itype(k,1).ne.10) then
4862 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
4867 ! Calculate deviations from standard virtual-bond lengths
4871 x(ind)=vbld(i+1)**2-vbl**2
4874 if (itype(i,1).ne.10) then
4876 x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4880 write (iout,*) "Coordinates and violations"
4882 write(iout,'(i5,3f10.5,5x,e15.5)') &
4883 i,(dC_uncor(j,i),j=1,3),x(i)
4885 write (iout,*) "Velocities and violations"
4889 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4890 i,ind,(d_t(j,i),j=1,3),scalar(d_t(1,i),dC_old(1,i))
4893 if (itype(i,1).ne.10) then
4895 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4896 i+nres,ind,(d_t(j,i+nres),j=1,3),&
4897 scalar(d_t(1,i+nres),dC_old(1,i+nres))
4900 write (iout,*) "gdc"
4902 write (iout,*) "i",i
4904 write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4910 if (dabs(x(i)).gt.xmax) then
4914 if (xmax.lt.tol_rattle) then
4918 ! Calculate the matrix of the system of equations
4923 Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
4928 write (iout,*) "Matrix Cmat"
4929 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4931 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4932 ! Add constraint term to positions
4939 xx = xx+x(ii)*gdc(j,ind,ii)
4942 d_t(j,i)=d_t(j,i)+xx/d_time
4947 if (itype(i,1).ne.10) then
4952 xx = xx+x(ii)*gdc(j,ind,ii)
4955 d_t(j,i+nres)=d_t(j,i+nres)+xx/d_time
4956 dC(j,i+nres)=dC(j,i+nres)+xx
4960 ! Rebuild the chain using the new coordinates
4961 call chainbuild_cart
4963 write (iout,*) "New coordinates, Lagrange multipliers,",&
4964 " and differences between actual and standard bond lengths"
4968 xx=vbld(i+1)**2-vbl**2
4969 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4970 i,(dC(j,i),j=1,3),x(ind),xx
4973 if (itype(i,1).ne.10) then
4975 xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4976 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4977 i,(dC(j,i+nres),j=1,3),x(ind),xx
4980 write (iout,*) "Velocities and violations"
4984 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4985 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4988 if (itype(i,1).ne.10) then
4990 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4991 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4992 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4999 10 write (iout,*) "Error - singularity in solving the system",&
5000 " of equations for Lagrange multipliers."
5004 "RATTLE inactive; use -DRATTLE option at compile time"
5007 end subroutine rattle_brown
5008 !-----------------------------------------------------------------------------
5010 !-----------------------------------------------------------------------------
5011 subroutine friction_force
5016 ! implicit real*8 (a-h,o-z)
5017 ! include 'DIMENSIONS'
5018 ! include 'COMMON.VAR'
5019 ! include 'COMMON.CHAIN'
5020 ! include 'COMMON.DERIV'
5021 ! include 'COMMON.GEO'
5022 ! include 'COMMON.LOCAL'
5023 ! include 'COMMON.INTERACT'
5024 ! include 'COMMON.MD'
5026 ! include 'COMMON.LANGEVIN'
5028 ! include 'COMMON.LANGEVIN.lang0'
5030 ! include 'COMMON.IOUNITS'
5031 !el real(kind=8),dimension(6*nres) :: gamvec !(MAXRES6) maxres6=6*maxres
5032 !el common /syfek/ gamvec
5033 real(kind=8) :: vv(3),vvtot(3,nres),v_work(6*nres) !,&
5034 !el ginvfric(2*nres,2*nres) !maxres2=2*maxres
5035 !el common /przechowalnia/ ginvfric
5037 logical :: lprn, checkmode
5038 integer :: i,j,ind,k,nres2,nres6,mnum
5043 ! if (large) lprn=.true.
5044 ! if (large) checkmode=.true.
5045 if(.not.allocated(gamvec)) allocate(gamvec(nres6)) !(MAXRES6)
5046 if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
5054 d_t_work(j)=d_t(j,0)
5059 d_t_work(ind+j)=d_t(j,i)
5065 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5066 .and.(mnum.ne.5)) then
5068 d_t_work(ind+j)=d_t(j,i+nres)
5074 call fricmat_mult(d_t_work,fric_work)
5076 if (.not.checkmode) return
5079 write (iout,*) "d_t_work and fric_work"
5081 write (iout,'(i3,2e15.5)') i,d_t_work(i),fric_work(i)
5085 friction(j,0)=fric_work(j)
5090 friction(j,i)=fric_work(ind+j)
5096 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5097 .and.(mnum.ne.5)) then
5098 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
5100 friction(j,i+nres)=fric_work(ind+j)
5106 write(iout,*) "Friction backbone"
5108 write(iout,'(i5,3e15.5,5x,3e15.5)') &
5109 i,(friction(j,i),j=1,3),(d_t(j,i),j=1,3)
5111 write(iout,*) "Friction side chain"
5113 write(iout,'(i5,3e15.5,5x,3e15.5)') &
5114 i,(friction(j,i+nres),j=1,3),(d_t(j,i+nres),j=1,3)
5123 vvtot(j,i)=vv(j)+0.5d0*d_t(j,i)
5124 vvtot(j,i+nres)=vv(j)+d_t(j,i+nres)
5125 vv(j)=vv(j)+d_t(j,i)
5128 write (iout,*) "vvtot backbone and sidechain"
5130 write (iout,'(i5,3e15.5,5x,3e15.5)') i,(vvtot(j,i),j=1,3),&
5131 (vvtot(j,i+nres),j=1,3)
5136 v_work(ind+j)=vvtot(j,i)
5142 v_work(ind+j)=vvtot(j,i+nres)
5146 write (iout,*) "v_work gamvec and site-based friction forces"
5148 write (iout,'(i5,3e15.5)') i,v_work(i),gamvec(i),&
5152 ! fric_work1(i)=0.0d0
5154 ! fric_work1(i)=fric_work1(i)-A(j,i)*gamvec(j)*v_work(j)
5157 ! write (iout,*) "fric_work and fric_work1"
5159 ! write (iout,'(i5,2e15.5)') i,fric_work(i),fric_work1(i)
5165 ginvfric(i,j)=ginvfric(i,j)+ginv(i,k)*fricmat(k,j)
5169 write (iout,*) "ginvfric"
5171 write (iout,'(i5,100f8.3)') i,(ginvfric(i,j),j=1,dimen)
5173 write (iout,*) "symmetry check"
5176 write (iout,*) i,j,ginvfric(i,j)-ginvfric(j,i)
5181 end subroutine friction_force
5182 !-----------------------------------------------------------------------------
5183 subroutine setup_fricmat
5187 use control_data, only:time_Bcast
5188 use control, only:tcpu
5190 ! implicit real*8 (a-h,o-z)
5194 real(kind=8) :: time00
5196 ! include 'DIMENSIONS'
5197 ! include 'COMMON.VAR'
5198 ! include 'COMMON.CHAIN'
5199 ! include 'COMMON.DERIV'
5200 ! include 'COMMON.GEO'
5201 ! include 'COMMON.LOCAL'
5202 ! include 'COMMON.INTERACT'
5203 ! include 'COMMON.MD'
5204 ! include 'COMMON.SETUP'
5205 ! include 'COMMON.TIME1'
5206 ! integer licznik /0/
5209 ! include 'COMMON.LANGEVIN'
5211 ! include 'COMMON.LANGEVIN.lang0'
5213 ! include 'COMMON.IOUNITS'
5215 integer :: i,j,ind,ind1,m
5216 logical :: lprn = .false.
5217 real(kind=8) :: dtdi !el ,gamvec(2*nres)
5218 !el real(kind=8),dimension(2*nres,2*nres) :: ginvfric,fcopy
5219 ! real(kind=8),allocatable,dimension(:,:) :: fcopy
5220 !el real(kind=8),dimension(2*nres*(2*nres+1)/2) :: Ghalf !(mmaxres2) (mmaxres2=(maxres2*(maxres2+1)/2))
5221 !el common /syfek/ gamvec
5222 real(kind=8) :: work(8*2*nres)
5223 integer :: iwork(2*nres)
5224 !el common /przechowalnia/ ginvfric,Ghalf,fcopy
5225 integer :: ii,iti,k,l,nzero,nres2,nres6,ierr,mnum
5229 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
5230 if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
5231 if (fg_rank.ne.king) goto 10
5236 if(.not.allocated(gamvec)) allocate(gamvec(nres2)) !(MAXRES2)
5237 if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
5238 if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
5239 !el allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
5240 if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !maxres2=2*maxres
5242 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
5243 ! Zeroing out fricmat
5249 ! Load the friction coefficients corresponding to peptide groups
5254 gamvec(ind1)=gamp(mnum)
5256 ! Load the friction coefficients corresponding to side chains
5260 gamsc(ntyp1_molec(j),j)=1.0d0
5267 gamvec(ii)=gamsc(iabs(iti),mnum)
5269 if (surfarea) call sdarea(gamvec)
5271 ! write (iout,*) "Matrix A and vector gamma"
5273 ! write (iout,'(i2,$)') i
5275 ! write (iout,'(f4.1,$)') A(i,j)
5277 ! write (iout,'(f8.3)') gamvec(i)
5281 write (iout,*) "Vector gamvec"
5283 write (iout,'(i5,f10.5)') i, gamvec(i)
5287 ! The friction matrix
5292 dtdi=dtdi+A(j,k)*A(j,i)*gamvec(j)
5299 write (iout,'(//a)') "Matrix fricmat"
5300 call matout2(dimen,dimen,nres2,nres2,fricmat)
5302 if (lang.eq.2 .or. lang.eq.3) then
5303 ! Mass-scale the friction matrix if non-direct integration will be performed
5309 Ginvfric(i,j)=Ginvfric(i,j)+ &
5310 Gsqrm(i,k)*Gsqrm(l,j)*fricmat(k,l)
5315 ! Diagonalize the friction matrix
5320 Ghalf(ind)=Ginvfric(i,j)
5323 call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
5326 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
5327 " mass-scaled friction matrix"
5328 call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
5330 ! Precompute matrices for tinker stochastic integrator
5337 mt1(i,j)=mt1(i,j)+fricvec(k,i)*gsqrm(k,j)
5338 mt2(i,j)=mt2(i,j)+fricvec(k,i)*gsqrp(k,j)
5344 else if (lang.eq.4) then
5345 ! Diagonalize the friction matrix
5350 Ghalf(ind)=fricmat(i,j)
5353 call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
5356 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
5358 call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
5360 ! Determine the number of zero eigenvalues of the friction matrix
5361 nzero=max0(dimen-dimen1,0)
5362 ! do while (fricgam(nzero+1).le.1.0d-5 .and. nzero.lt.dimen)
5365 write (iout,*) "Number of zero eigenvalues:",nzero
5370 fricmat(i,j)=fricmat(i,j) &
5371 +fricvec(i,k)*fricvec(j,k)/fricgam(k)
5376 write (iout,'(//a)') "Generalized inverse of fricmat"
5377 call matout(dimen,dimen,nres6,nres6,fricmat)
5382 if (nfgtasks.gt.1) then
5383 if (fg_rank.eq.0) then
5384 ! The matching BROADCAST for fg processors is called in ERGASTULUM
5390 call MPI_Bcast(10,1,MPI_INTEGER,king,FG_COMM,IERROR)
5392 time_Bcast=time_Bcast+MPI_Wtime()-time00
5394 time_Bcast=time_Bcast+tcpu()-time00
5396 ! print *,"Processor",myrank,
5397 ! & " BROADCAST iorder in SETUP_FRICMAT"
5400 write (iout,*) "setup_fricmat licznik"!,licznik !sp
5406 ! Scatter the friction matrix
5407 call MPI_Scatterv(fricmat(1,1),nginv_counts(0),&
5408 nginv_start(0),MPI_DOUBLE_PRECISION,fcopy(1,1),&
5409 myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
5412 time_scatter=time_scatter+MPI_Wtime()-time00
5413 time_scatter_fmat=time_scatter_fmat+MPI_Wtime()-time00
5415 time_scatter=time_scatter+tcpu()-time00
5416 time_scatter_fmat=time_scatter_fmat+tcpu()-time00
5420 do j=1,2*my_ng_count
5421 fricmat(j,i)=fcopy(i,j)
5424 ! write (iout,*) "My chunk of fricmat"
5425 ! call MATOUT2(my_ng_count,dimen,maxres2,maxres2,fcopy)
5429 end subroutine setup_fricmat
5430 !-----------------------------------------------------------------------------
5431 subroutine stochastic_force(stochforcvec)
5434 use random, only:anorm_distr
5435 ! implicit real*8 (a-h,o-z)
5436 ! include 'DIMENSIONS'
5437 use control, only: tcpu
5442 ! include 'COMMON.VAR'
5443 ! include 'COMMON.CHAIN'
5444 ! include 'COMMON.DERIV'
5445 ! include 'COMMON.GEO'
5446 ! include 'COMMON.LOCAL'
5447 ! include 'COMMON.INTERACT'
5448 ! include 'COMMON.MD'
5449 ! include 'COMMON.TIME1'
5451 ! include 'COMMON.LANGEVIN'
5453 ! include 'COMMON.LANGEVIN.lang0'
5455 ! include 'COMMON.IOUNITS'
5457 real(kind=8) :: x,sig,lowb,highb
5458 real(kind=8) :: ff(3),force(3,0:2*nres),zeta2,lowb2
5459 real(kind=8) :: highb2,sig2,forcvec(6*nres),stochforcvec(6*nres)
5460 real(kind=8) :: time00
5461 logical :: lprn = .false.
5462 integer :: i,j,ind,mnum
5466 stochforc(j,i)=0.0d0
5476 ! Compute the stochastic forces acting on bodies. Store in force.
5482 force(j,i)=anorm_distr(x,sig,lowb,highb)
5490 force(j,i+nres)=anorm_distr(x,sig2,lowb2,highb2)
5494 time_fsample=time_fsample+MPI_Wtime()-time00
5496 time_fsample=time_fsample+tcpu()-time00
5498 ! Compute the stochastic forces acting on virtual-bond vectors.
5504 stochforc(j,i)=ff(j)+0.5d0*force(j,i)
5507 ff(j)=ff(j)+force(j,i)
5509 ! if (itype(i+1,1).ne.ntyp1) then
5511 if (itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
5513 stochforc(j,i)=stochforc(j,i)+force(j,i+nres+1)
5514 ff(j)=ff(j)+force(j,i+nres+1)
5519 stochforc(j,0)=ff(j)+force(j,nnt+nres)
5523 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5524 .and.(mnum.ne.5)) then
5525 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
5527 stochforc(j,i+nres)=force(j,i+nres)
5533 stochforcvec(j)=stochforc(j,0)
5538 stochforcvec(ind+j)=stochforc(j,i)
5544 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5545 .and.(mnum.ne.5)) then
5546 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
5548 stochforcvec(ind+j)=stochforc(j,i+nres)
5554 write (iout,*) "stochforcvec"
5556 write(iout,'(i5,e15.5)') i,stochforcvec(i)
5558 write(iout,*) "Stochastic forces backbone"
5560 write(iout,'(i5,3e15.5)') i,(stochforc(j,i),j=1,3)
5562 write(iout,*) "Stochastic forces side chain"
5564 write(iout,'(i5,3e15.5)') &
5565 i,(stochforc(j,i+nres),j=1,3)
5573 write (iout,*) i,ind
5575 forcvec(ind+j)=force(j,i)
5580 write (iout,*) i,ind
5582 forcvec(j+ind)=force(j,i+nres)
5587 write (iout,*) "forcvec"
5591 write (iout,'(2i3,2f10.5)') i,j,force(j,i),&
5598 write (iout,'(2i3,2f10.5)') i,j,force(j,i+nres),&
5607 end subroutine stochastic_force
5608 !-----------------------------------------------------------------------------
5609 subroutine sdarea(gamvec)
5611 ! Scale the friction coefficients according to solvent accessible surface areas
5612 ! Code adapted from TINKER
5616 ! implicit real*8 (a-h,o-z)
5617 ! include 'DIMENSIONS'
5618 ! include 'COMMON.CONTROL'
5619 ! include 'COMMON.VAR'
5620 ! include 'COMMON.MD'
5622 ! include 'COMMON.LANGEVIN'
5624 ! include 'COMMON.LANGEVIN.lang0'
5626 ! include 'COMMON.CHAIN'
5627 ! include 'COMMON.DERIV'
5628 ! include 'COMMON.GEO'
5629 ! include 'COMMON.LOCAL'
5630 ! include 'COMMON.INTERACT'
5631 ! include 'COMMON.IOUNITS'
5632 ! include 'COMMON.NAMES'
5633 real(kind=8),dimension(2*nres) :: radius,gamvec !(maxres2)
5634 real(kind=8),parameter :: twosix = 1.122462048309372981d0
5635 logical :: lprn = .false.
5636 real(kind=8) :: probe,area,ratio
5637 integer :: i,j,ind,iti,mnum
5639 ! determine new friction coefficients every few SD steps
5641 ! set the atomic radii to estimates of sigma values
5643 ! print *,"Entered sdarea"
5649 ! Load peptide group radii
5652 radius(i)=pstok(mnum)
5654 ! Load side chain radii
5658 radius(i+nres)=restok(iti,mnum)
5661 ! write (iout,*) "i",i," radius",radius(i)
5664 radius(i) = radius(i) / twosix
5665 if (radius(i) .ne. 0.0d0) radius(i) = radius(i) + probe
5668 ! scale atomic friction coefficients by accessible area
5670 if (lprn) write (iout,*) &
5671 "Original gammas, surface areas, scaling factors, new gammas, ",&
5672 "std's of stochastic forces"
5675 if (radius(i).gt.0.0d0) then
5676 call surfatom (i,area,radius)
5677 ratio = dmax1(area/(4.0d0*pi*radius(i)**2),1.0d-1)
5678 if (lprn) write (iout,'(i5,3f10.5,$)') &
5679 i,gamvec(ind+1),area,ratio
5682 gamvec(ind) = ratio * gamvec(ind)
5685 stdforcp(i)=stdfp(mnum)*dsqrt(gamvec(ind))
5686 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcp(i)
5690 if (radius(i+nres).gt.0.0d0) then
5691 call surfatom (i+nres,area,radius)
5692 ratio = dmax1(area/(4.0d0*pi*radius(i+nres)**2),1.0d-1)
5693 if (lprn) write (iout,'(i5,3f10.5,$)') &
5694 i,gamvec(ind+1),area,ratio
5697 gamvec(ind) = ratio * gamvec(ind)
5700 stdforcsc(i)=stdfsc(itype(i,mnum),mnum)*dsqrt(gamvec(ind))
5701 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcsc(i)
5706 end subroutine sdarea
5707 !-----------------------------------------------------------------------------
5709 !-----------------------------------------------------------------------------
5712 ! ###################################################
5713 ! ## COPYRIGHT (C) 1996 by Jay William Ponder ##
5714 ! ## All Rights Reserved ##
5715 ! ###################################################
5717 ! ################################################################
5719 ! ## subroutine surfatom -- exposed surface area of an atom ##
5721 ! ################################################################
5724 ! "surfatom" performs an analytical computation of the surface
5725 ! area of a specified atom; a simplified version of "surface"
5727 ! literature references:
5729 ! T. J. Richmond, "Solvent Accessible Surface Area and
5730 ! Excluded Volume in Proteins", Journal of Molecular Biology,
5733 ! L. Wesson and D. Eisenberg, "Atomic Solvation Parameters
5734 ! Applied to Molecular Dynamics of Proteins in Solution",
5735 ! Protein Science, 1, 227-235 (1992)
5737 ! variables and parameters:
5739 ! ir number of atom for which area is desired
5740 ! area accessible surface area of the atom
5741 ! radius radii of each of the individual atoms
5744 subroutine surfatom(ir,area,radius)
5746 ! implicit real*8 (a-h,o-z)
5747 ! include 'DIMENSIONS'
5749 ! include 'COMMON.GEO'
5750 ! include 'COMMON.IOUNITS'
5752 integer :: nsup,nstart_sup
5753 ! double precision c,dc,dc_old,d_c_work,xloc,xrot,dc_norm
5754 ! common /chain/ c(3,maxres2+2),dc(3,0:maxres2),dc_old(3,0:maxres2),
5755 ! & xloc(3,maxres),xrot(3,maxres),dc_norm(3,0:maxres2),
5756 ! & dc_work(MAXRES6),nres,nres0
5757 integer,parameter :: maxarc=300
5761 integer :: mi,ni,narc
5762 integer :: key(maxarc)
5763 integer :: intag(maxarc)
5764 integer :: intag1(maxarc)
5765 real(kind=8) :: area,arcsum
5766 real(kind=8) :: arclen,exang
5767 real(kind=8) :: delta,delta2
5768 real(kind=8) :: eps,rmove
5769 real(kind=8) :: xr,yr,zr
5770 real(kind=8) :: rr,rrsq
5771 real(kind=8) :: rplus,rminus
5772 real(kind=8) :: axx,axy,axz
5773 real(kind=8) :: ayx,ayy
5774 real(kind=8) :: azx,azy,azz
5775 real(kind=8) :: uxj,uyj,uzj
5776 real(kind=8) :: tx,ty,tz
5777 real(kind=8) :: txb,tyb,td
5778 real(kind=8) :: tr2,tr,txr,tyr
5779 real(kind=8) :: tk1,tk2
5780 real(kind=8) :: thec,the,t,tb
5781 real(kind=8) :: txk,tyk,tzk
5782 real(kind=8) :: t1,ti,tf,tt
5783 real(kind=8) :: txj,tyj,tzj
5784 real(kind=8) :: ccsq,cc,xysq
5785 real(kind=8) :: bsqk,bk,cosine
5786 real(kind=8) :: dsqj,gi,pix2
5787 real(kind=8) :: therk,dk,gk
5788 real(kind=8) :: risqk,rik
5789 real(kind=8) :: radius(2*nres) !(maxatm) (maxatm=maxres2)
5790 real(kind=8) :: ri(maxarc),risq(maxarc)
5791 real(kind=8) :: ux(maxarc),uy(maxarc),uz(maxarc)
5792 real(kind=8) :: xc(maxarc),yc(maxarc),zc(maxarc)
5793 real(kind=8) :: xc1(maxarc),yc1(maxarc),zc1(maxarc)
5794 real(kind=8) :: dsq(maxarc),bsq(maxarc)
5795 real(kind=8) :: dsq1(maxarc),bsq1(maxarc)
5796 real(kind=8) :: arci(maxarc),arcf(maxarc)
5797 real(kind=8) :: ex(maxarc),lt(maxarc),gr(maxarc)
5798 real(kind=8) :: b(maxarc),b1(maxarc),bg(maxarc)
5799 real(kind=8) :: kent(maxarc),kout(maxarc)
5800 real(kind=8) :: ther(maxarc)
5801 logical :: moved,top
5802 logical :: omit(maxarc)
5805 ! maxatm = 2*nres !maxres2 maxres2=2*maxres
5806 ! maxlight = 8*maxatm
5809 ! maxtors = 4*maxatm
5811 ! zero out the surface area for the sphere of interest
5814 ! write (2,*) "ir",ir," radius",radius(ir)
5815 if (radius(ir) .eq. 0.0d0) return
5817 ! set the overlap significance and connectivity shift
5821 delta2 = delta * delta
5826 ! store coordinates and radius of the sphere of interest
5834 ! initialize values of some counters and summations
5843 ! test each sphere to see if it overlaps the sphere of interest
5846 if (i.eq.ir .or. radius(i).eq.0.0d0) goto 30
5847 rplus = rr + radius(i)
5849 if (abs(tx) .ge. rplus) goto 30
5851 if (abs(ty) .ge. rplus) goto 30
5853 if (abs(tz) .ge. rplus) goto 30
5855 ! check for sphere overlap by testing distance against radii
5857 xysq = tx*tx + ty*ty
5858 if (xysq .lt. delta2) then
5865 if (rplus-cc .le. delta) goto 30
5866 rminus = rr - radius(i)
5868 ! check to see if sphere of interest is completely buried
5870 if (cc-abs(rminus) .le. delta) then
5871 if (rminus .le. 0.0d0) goto 170
5875 ! check for too many overlaps with sphere of interest
5877 if (io .ge. maxarc) then
5879 20 format (/,' SURFATOM -- Increase the Value of MAXARC')
5883 ! get overlap between current sphere and sphere of interest
5892 gr(io) = (ccsq+rplus*rminus) / (2.0d0*rr*b1(io))
5898 ! case where no other spheres overlap the sphere of interest
5901 area = 4.0d0 * pi * rrsq
5905 ! case where only one sphere overlaps the sphere of interest
5908 area = pix2 * (1.0d0 + gr(1))
5909 area = mod(area,4.0d0*pi) * rrsq
5913 ! case where many spheres intersect the sphere of interest;
5914 ! sort the intersecting spheres by their degree of overlap
5916 call sort2 (io,gr,key)
5919 intag(i) = intag1(k)
5928 ! get radius of each overlap circle on surface of the sphere
5933 risq(i) = rrsq - gi*gi
5934 ri(i) = sqrt(risq(i))
5935 ther(i) = 0.5d0*pi - asin(min(1.0d0,max(-1.0d0,gr(i))))
5938 ! find boundary of inaccessible area on sphere of interest
5941 if (.not. omit(k)) then
5948 ! check to see if J circle is intersecting K circle;
5949 ! get distance between circle centers and sum of radii
5952 if (omit(j)) goto 60
5953 cc = (txk*xc(j)+tyk*yc(j)+tzk*zc(j))/(bk*b(j))
5954 cc = acos(min(1.0d0,max(-1.0d0,cc)))
5955 td = therk + ther(j)
5957 ! check to see if circles enclose separate regions
5959 if (cc .ge. td) goto 60
5961 ! check for circle J completely inside circle K
5963 if (cc+ther(j) .lt. therk) goto 40
5965 ! check for circles that are essentially parallel
5967 if (cc .gt. delta) goto 50
5972 ! check to see if sphere of interest is completely buried
5975 if (pix2-cc .le. td) goto 170
5981 ! find T value of circle intersections
5984 if (omit(k)) goto 110
5999 ! rotation matrix elements
6011 if (.not. omit(j)) then
6016 ! rotate spheres so K vector colinear with z-axis
6018 uxj = txj*axx + tyj*axy - tzj*axz
6019 uyj = tyj*ayy - txj*ayx
6020 uzj = txj*azx + tyj*azy + tzj*azz
6021 cosine = min(1.0d0,max(-1.0d0,uzj/b(j)))
6022 if (acos(cosine) .lt. therk+ther(j)) then
6023 dsqj = uxj*uxj + uyj*uyj
6028 tr2 = risqk*dsqj - tb*tb
6034 ! get T values of intersection for K circle
6037 tb = min(1.0d0,max(-1.0d0,tb))
6039 if (tyb-txr .lt. 0.0d0) tk1 = pix2 - tk1
6041 tb = min(1.0d0,max(-1.0d0,tb))
6043 if (tyb+txr .lt. 0.0d0) tk2 = pix2 - tk2
6044 thec = (rrsq*uzj-gk*bg(j)) / (rik*ri(j)*b(j))
6045 if (abs(thec) .lt. 1.0d0) then
6047 else if (thec .ge. 1.0d0) then
6049 else if (thec .le. -1.0d0) then
6053 ! see if "tk1" is entry or exit point; check t=0 point;
6054 ! "ti" is exit point, "tf" is entry point
6056 cosine = min(1.0d0,max(-1.0d0, &
6057 (uzj*gk-uxj*rik)/(b(j)*rr)))
6058 if ((acos(cosine)-ther(j))*(tk2-tk1) .le. 0.0d0) then
6066 if (narc .ge. maxarc) then
6068 70 format (/,' SURFATOM -- Increase the Value',&
6072 if (tf .le. ti) then
6093 ! special case; K circle without intersections
6095 if (narc .le. 0) goto 90
6097 ! general case; sum up arclength and set connectivity code
6099 call sort2 (narc,arci,key)
6104 if (narc .gt. 1) then
6107 if (t .lt. arci(j)) then
6108 arcsum = arcsum + arci(j) - t
6109 exang = exang + ex(ni)
6111 if (jb .ge. maxarc) then
6113 80 format (/,' SURFATOM -- Increase the Value',&
6118 kent(jb) = maxarc*i + k
6120 kout(jb) = maxarc*k + i
6129 arcsum = arcsum + pix2 - t
6131 exang = exang + ex(ni)
6134 kent(jb) = maxarc*i + k
6136 kout(jb) = maxarc*k + i
6143 arclen = arclen + gr(k)*arcsum
6146 if (arclen .eq. 0.0d0) goto 170
6147 if (jb .eq. 0) goto 150
6149 ! find number of independent boundaries and check connectivity
6153 if (kout(k) .ne. 0) then
6160 if (m .eq. kent(ii)) then
6163 if (j .eq. jb) goto 150
6175 ! attempt to fix connectivity error by moving atom slightly
6179 140 format (/,' SURFATOM -- Connectivity Error at Atom',i6)
6188 ! compute the exposed surface area for the sphere of interest
6191 area = ib*pix2 + exang + arclen
6192 area = mod(area,4.0d0*pi) * rrsq
6194 ! attempt to fix negative area by moving atom slightly
6196 if (area .lt. 0.0d0) then
6199 160 format (/,' SURFATOM -- Negative Area at Atom',i6)
6210 end subroutine surfatom
6211 !----------------------------------------------------------------
6212 !----------------------------------------------------------------
6213 subroutine alloc_MD_arrays
6214 !EL Allocation of arrays used by MD module
6216 integer :: nres2,nres6
6219 !----------------------
6223 allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
6224 allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
6225 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
6226 allocate(fricvec(nres2,nres2))
6227 allocate(pfric_mat(nres2,nres2),vfric_mat(nres2,nres2))
6228 allocate(afric_mat(nres2,nres2),prand_mat(nres2,nres2))
6229 allocate(vrand_mat1(nres2,nres2),vrand_mat2(nres2,nres2)) !(MAXRES2,MAXRES2)
6230 allocate(pfric0_mat(nres2,nres2,0:maxflag_stoch))
6231 allocate(afric0_mat(nres2,nres2,0:maxflag_stoch))
6232 allocate(vfric0_mat(nres2,nres2,0:maxflag_stoch))
6233 allocate(prand0_mat(nres2,nres2,0:maxflag_stoch))
6234 allocate(vrand0_mat1(nres2,nres2,0:maxflag_stoch))
6235 allocate(vrand0_mat2(nres2,nres2,0:maxflag_stoch)) !(MAXRES2,MAXRES2,0:maxflag_stoch)
6236 allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
6238 allocate(mt1(nres2,nres2),mt2(nres2,nres2),mt3(nres2,nres2)) !(maxres2,maxres2)
6239 !----------------------
6241 ! commom.langevin.lang0
6243 allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
6244 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
6245 allocate(fricvec(nres2,nres2)) !(MAXRES2,MAXRES2)
6246 allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
6247 allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
6250 if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2))
6251 !----------------------
6252 ! commom.hairpin in CSA module
6253 !----------------------
6254 ! common.mce in MCM_MD module
6255 !----------------------
6257 ! common /mdgrad/ in module.energy
6258 ! common /back_constr/ in module.energy
6259 ! common /qmeas/ in module.energy
6262 allocate(potEcomp(0:n_ene+4)) !(0:n_ene+4)
6264 allocate(d_t(3,0:nres2),d_a(3,0:nres2),d_t_old(3,0:nres2)) !(3,0:MAXRES2)
6265 allocate(d_a_work(nres6)) !(6*MAXRES)
6267 allocate(DM(nres2),DU1(nres2),DU2(nres2))
6268 allocate(DMorig(nres2),DU1orig(nres2),DU2orig(nres2))
6270 write (iout,*) "Before A Allocation",nfgtasks-1
6272 allocate(Gmat(nres2,nres2),A(nres2,nres2))
6273 if(.not.allocated(Ginv)) allocate(Ginv(nres2,nres2)) !in control: ergastulum
6274 allocate(Gsqrp(nres2,nres2),Gsqrm(nres2,nres2),Gvec(nres2,nres2)) !(maxres2,maxres2)
6276 allocate(Geigen(nres2)) !(maxres2)
6277 if(.not.allocated(vtot)) allocate(vtot(nres2)) !(maxres2)
6278 ! common /inertia/ in io_conf: parmread
6279 ! real(kind=8),dimension(:),allocatable :: ISC,msc !(ntyp+1)
6280 ! common /langevin/in io read_MDpar
6281 ! real(kind=8),dimension(:),allocatable :: gamsc !(ntyp1)
6282 ! real(kind=8),dimension(:),allocatable :: stdfsc !(ntyp)
6283 ! in io_conf: parmread
6284 ! real(kind=8),dimension(:),allocatable :: restok !(ntyp+1)
6285 ! common /mdpmpi/ in control: ergastulum
6286 if(.not.allocated(ng_start)) allocate(ng_start(0:nfgtasks-1))
6287 if(.not.allocated(ng_counts)) allocate(ng_counts(0:nfgtasks-1))
6288 if(.not.allocated(nginv_counts)) allocate(nginv_counts(0:nfgtasks-1)) !(0:MaxProcs-1)
6289 if(.not.allocated(nginv_start)) allocate(nginv_start(0:nfgtasks)) !(0:MaxProcs)
6290 !----------------------
6291 ! common.muca in read_muca
6292 ! common /double_muca/
6293 ! real(kind=8) :: elow,ehigh,factor,hbin,factor_min
6294 ! real(kind=8),dimension(:),allocatable :: emuca,nemuca,&
6295 ! nemuca2,hist !(4*maxres)
6296 ! common /integer_muca/
6297 ! integer :: nmuca,imtime,muca_smooth
6299 ! real(kind=8),dimension(:),allocatable :: elowi,ehighi !(maxprocs)
6300 !----------------------
6302 ! common /mdgrad/ in module.energy
6303 ! common /back_constr/ in module.energy
6304 ! common /qmeas/ in module.energy
6308 allocate(d_t_work(nres6),d_t_work_new(nres6),d_af_work(nres6))
6309 allocate(d_as_work(nres6),kinetic_force(nres6)) !(MAXRES6)
6310 allocate(d_t_new(3,0:nres2),d_a_old(3,0:nres2),d_a_short(3,0:nres2)) !,d_a !(3,0:MAXRES2)
6311 allocate(stdforcp(nres),stdforcsc(nres)) !(MAXRES)
6312 !----------------------
6314 allocate(D_ban(nres6)) !(MAXRES6) maxres6=6*maxres
6315 ! common /stochcalc/ stochforcvec
6316 allocate(stochforcvec(nres6)) !(MAXRES6) maxres6=6*maxres
6317 !----------------------
6319 end subroutine alloc_MD_arrays
6320 !-----------------------------------------------------------------------------
6321 !-----------------------------------------------------------------------------