2 !-----------------------------------------------------------------------------
15 !-----------------------------------------------------------------------------
17 ! common /mdgrad/ in module.energy
18 ! common /back_constr/ in module.energy
19 ! common /qmeas/ in module.energy
23 real(kind=8),dimension(:),allocatable :: d_t_work,&
24 d_t_work_new,d_af_work,d_as_work,kinetic_force !(MAXRES6)
25 real(kind=8),dimension(:,:),allocatable :: d_t_new,&
26 d_a_old,d_a_short!,d_a !(3,0:MAXRES2)
27 ! real(kind=8),dimension(:),allocatable :: d_a_work !(6*MAXRES)
28 ! real(kind=8),dimension(:,:),allocatable :: Gmat,Ginv,A,&
29 ! Gsqrp,Gsqrm,Gvec !(maxres2,maxres2)
30 ! real(kind=8),dimension(:),allocatable :: Geigen !(maxres2)
31 ! integer :: dimen,dimen1,dimen3
32 ! integer :: lang,count_reset_moment,count_reset_vel
33 ! logical :: reset_moment,reset_vel,rattle,RESPA
36 ! real(kind=8) :: rwat,etawat,stdfp,cPoise
37 ! real(kind=8),dimension(:),allocatable :: gamsc !(ntyp1)
38 ! real(kind=8),dimension(:),allocatable :: stdfsc !(ntyp)
39 real(kind=8),dimension(:),allocatable :: stdforcp,stdforcsc !(MAXRES)
40 !-----------------------------------------------------------------------------
44 ! ###################################################
45 ! ## COPYRIGHT (C) 1992 by Jay William Ponder ##
46 ! ## All Rights Reserved ##
47 ! ###################################################
49 ! #############################################################
51 ! ## sizes.i -- parameter values to set array dimensions ##
53 ! #############################################################
56 ! "sizes.i" sets values for critical array dimensions used
57 ! throughout the software; these parameters will fix the size
58 ! of the largest systems that can be handled; values too large
59 ! for the computer's memory and/or swap space to accomodate
60 ! will result in poor performance or outright failure
62 ! parameter: maximum allowed number of:
64 ! maxatm atoms in the molecular system
65 ! maxval atoms directly bonded to an atom
66 ! maxgrp ! user-defined groups of atoms
67 ! maxtyp force field atom type definitions
68 ! maxclass force field atom class definitions
69 ! maxkey lines in the keyword file
70 ! maxrot bonds for torsional rotation
71 ! maxvar optimization variables (vector storage)
72 ! maxopt optimization variables (matrix storage)
73 ! maxhess off-diagonal Hessian elements
74 ! maxlight sites for method of lights neighbors
75 ! maxvib vibrational frequencies
76 ! maxgeo distance geometry points
77 ! maxcell unit cells in replicated crystal
78 ! maxring 3-, 4-, or 5-membered rings
79 ! maxfix geometric restraints
80 ! maxbio biopolymer atom definitions
81 ! maxres residues in the macromolecule
82 ! maxamino amino acid residue types
83 ! maxnuc nucleic acid residue types
84 ! maxbnd covalent bonds in molecular system
85 ! maxang bond angles in molecular system
86 ! maxtors torsional angles in molecular system
87 ! maxpi atoms in conjugated pisystem
88 ! maxpib covalent bonds involving pisystem
89 ! maxpit torsional angles involving pisystem
92 !el integer maxatm,maxval,maxgrp
93 !el integer maxtyp,maxclass,maxkey
94 !el integer maxrot,maxopt
95 !el integer maxhess,maxlight,maxvib
96 !el integer maxgeo,maxcell,maxring
97 !el integer maxfix,maxbio
98 !el integer maxamino,maxnuc,maxbnd
99 !el integer maxang,maxtors,maxpi
100 !el integer maxpib,maxpit
101 ! integer :: maxatm !=2*nres !maxres2 maxres2=2*maxres
102 ! integer,parameter :: maxval=8
103 ! integer,parameter :: maxgrp=1000
104 ! integer,parameter :: maxtyp=3000
105 ! integer,parameter :: maxclass=500
106 ! integer,parameter :: maxkey=10000
107 ! integer,parameter :: maxrot=1000
108 ! integer,parameter :: maxopt=1000
109 ! integer,parameter :: maxhess=1000000
110 ! integer :: maxlight !=8*maxatm
111 ! integer,parameter :: maxvib=1000
112 ! integer,parameter :: maxgeo=1000
113 ! integer,parameter :: maxcell=10000
114 ! integer,parameter :: maxring=10000
115 ! integer,parameter :: maxfix=10000
116 ! integer,parameter :: maxbio=10000
117 ! integer,parameter :: maxamino=31
118 ! integer,parameter :: maxnuc=12
119 ! integer :: maxbnd !=2*maxatm
120 ! integer :: maxang !=3*maxatm
121 ! integer :: maxtors !=4*maxatm
122 ! integer,parameter :: maxpi=100
123 ! integer,parameter :: maxpib=2*maxpi
124 ! integer,parameter :: maxpit=4*maxpi
125 !-----------------------------------------------------------------------------
126 ! Maximum number of seed
127 ! integer,parameter :: max_seed=1
128 !-----------------------------------------------------------------------------
129 real(kind=8),dimension(:),allocatable :: stochforcvec !(MAXRES6) maxres6=6*maxres
130 ! common /stochcalc/ stochforcvec
131 !-----------------------------------------------------------------------------
132 ! common /przechowalnia/ subroutines: rattle1,rattle2,rattle_brown
133 real(kind=8),dimension(:,:),allocatable :: GGinv !(2*nres,2*nres) maxres2=2*maxres
134 real(kind=8),dimension(:,:,:),allocatable :: gdc !(3,2*nres,2*nres) maxres2=2*maxres
135 real(kind=8),dimension(:,:),allocatable :: Cmat !(2*nres,2*nres) maxres2=2*maxres
136 !-----------------------------------------------------------------------------
137 ! common /syfek/ subroutines: friction_force,setup_fricmat
138 !el real(kind=8),dimension(:),allocatable :: gamvec !(MAXRES6) or (MAXRES2)
139 !-----------------------------------------------------------------------------
140 ! common /przechowalnia/ subroutines: friction_force,setup_fricmat
141 real(kind=8),dimension(:,:),allocatable :: ginvfric !(2*nres,2*nres) !maxres2=2*maxres
142 !-----------------------------------------------------------------------------
143 ! common /przechowalnia/ subroutine: setup_fricmat
144 real(kind=8),dimension(:,:),allocatable :: fcopy !(2*nres,2*nres)
145 !-----------------------------------------------------------------------------
148 !-----------------------------------------------------------------------------
150 !-----------------------------------------------------------------------------
152 !-----------------------------------------------------------------------------
153 subroutine brown_step(itime)
154 !------------------------------------------------
155 ! Perform a single Euler integration step of Brownian dynamics
156 !------------------------------------------------
157 ! implicit real*8 (a-h,o-z)
159 use control, only: tcpu
162 ! use io_conf, only:cartprint
163 ! include 'DIMENSIONS'
167 ! include 'COMMON.CONTROL'
168 ! include 'COMMON.VAR'
169 ! include 'COMMON.MD'
171 ! include 'COMMON.LANGEVIN'
173 ! include 'COMMON.LANGEVIN.lang0'
175 ! include 'COMMON.CHAIN'
176 ! include 'COMMON.DERIV'
177 ! include 'COMMON.GEO'
178 ! include 'COMMON.LOCAL'
179 ! include 'COMMON.INTERACT'
180 ! include 'COMMON.IOUNITS'
181 ! include 'COMMON.NAMES'
182 ! include 'COMMON.TIME1'
183 real(kind=8),dimension(6*nres) :: zapas !(MAXRES6) maxres6=6*maxres
184 integer :: rstcount !ilen,
186 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
187 real(kind=8),dimension(6*nres,2*nres) :: Bmat,GBmat,Tmat !(MAXRES6,MAXRES2) (maxres2=2*maxres,maxres6=6*maxres)
188 real(kind=8),dimension(2*nres,2*nres) :: Cmat_,Cinv !(maxres2,maxres2) maxres2=2*maxres
189 real(kind=8),dimension(6*nres,6*nres) :: Pmat !(maxres6,maxres6) maxres6=6*maxres
190 ! real(kind=8),dimension(:,:),allocatable :: Bmat,GBmat,Tmat !(MAXRES6,MAXRES2) (maxres2=2*maxres,maxres6=6*maxres)
191 ! real(kind=8),dimension(:,:),allocatable :: Cmat_,Cinv !(maxres2,maxres2) maxres2=2*maxres
192 ! real(kind=8),dimension(:,:),allocatable :: Pmat !(maxres6,maxres6) maxres6=6*maxres
193 real(kind=8),dimension(6*nres) :: Td !(maxres6) maxres6=6*maxres
194 real(kind=8),dimension(2*nres) :: ppvec !(maxres2) maxres2=2*maxres
195 !el common /stochcalc/ stochforcvec
196 !el real(kind=8),dimension(3) :: cm !el
197 !el common /gucio/ cm
199 logical :: lprn = .false.,lprn1 = .false.
200 integer :: maxiter = 5
201 real(kind=8) :: difftol = 1.0d-5
202 real(kind=8) :: xx,diffmax,blen2,diffbond,tt0
203 integer :: i,j,nbond,k,ind,ind1,iter
204 integer :: nres2,nres6
208 ! if (.not.allocated(Bmat)) allocate(Bmat(nres6,nres2))
209 ! if (.not.allocated(GBmat)) allocate (GBmat(nres6,nres2))
210 ! if (.not.allocated(Tmat)) allocate (Tmat(nres6,nres2))
211 ! if (.not.allocated(Cmat_)) allocate(Cmat_(nres2,nres2))
212 ! if (.not.allocated(Cinv)) allocate (Cinv(nres2,nres2))
213 ! if (.not.allocated(Pmat)) allocate(Pmat(6*nres,6*nres))
215 if (.not.allocated(stochforcvec)) allocate(stochforcvec(nres6)) !(MAXRES6) maxres6=6*maxres
219 if (itype(i,1).ne.10) nbond=nbond+1
223 write (iout,*) "Generalized inverse of fricmat"
224 call matout(dimen,dimen,nres6,nres6,fricmat)
236 Bmat(ind+j,ind1)=dC_norm(j,i)
241 if (itype(i,1).ne.10) then
244 Bmat(ind+j,ind1)=dC_norm(j,i+nres)
250 write (iout,*) "Matrix Bmat"
251 call MATOUT(nbond,dimen,nres6,nres6,Bmat)
257 GBmat(i,j)=GBmat(i,j)+fricmat(i,k)*Bmat(k,j)
262 write (iout,*) "Matrix GBmat"
263 call MATOUT(nbond,dimen,nres6,nres2,Gbmat)
269 Cmat_(i,j)=Cmat_(i,j)+Bmat(k,i)*GBmat(k,j)
274 write (iout,*) "Matrix Cmat"
275 call MATOUT(nbond,nbond,nres2,nres2,Cmat_)
277 call matinvert(nbond,nres2,Cmat_,Cinv,osob)
279 write (iout,*) "Matrix Cinv"
280 call MATOUT(nbond,nbond,nres2,nres2,Cinv)
286 Tmat(i,j)=Tmat(i,j)+GBmat(i,k)*Cinv(k,j)
291 write (iout,*) "Matrix Tmat"
292 call MATOUT(nbond,dimen,nres6,nres2,Tmat)
302 Pmat(i,j)=Pmat(i,j)-Tmat(i,k)*Bmat(j,k)
307 write (iout,*) "Matrix Pmat"
308 call MATOUT(dimen,dimen,nres6,nres6,Pmat)
315 Td(i)=Td(i)+vbl*Tmat(i,ind)
318 if (itype(k,1).ne.10) then
320 Td(i)=Td(i)+vbldsc0(1,itype(k,1))*Tmat(i,ind)
325 write (iout,*) "Vector Td"
327 write (iout,'(i5,f10.5)') i,Td(i)
330 call stochastic_force(stochforcvec)
332 write (iout,*) "stochforcvec"
334 write (iout,*) i,stochforcvec(i)
338 zapas(j)=-gcart(j,0)+stochforcvec(j)
340 dC_work(j)=dC_old(j,0)
346 zapas(ind)=-gcart(j,i)+stochforcvec(ind)
347 dC_work(ind)=dC_old(j,i)
351 if (itype(i,1).ne.10) then
354 zapas(ind)=-gxcart(j,i)+stochforcvec(ind)
355 dC_work(ind)=dC_old(j,i+nres)
361 write (iout,*) "Initial d_t_work"
363 write (iout,*) i,d_t_work(i)
370 d_t_work(i)=d_t_work(i)+fricmat(i,j)*zapas(j)
377 zapas(i)=zapas(i)+Pmat(i,j)*(dC_work(j)+d_t_work(j)*d_time)
381 write (iout,*) "Final d_t_work and zapas"
383 write (iout,*) i,d_t_work(i),zapas(i)
397 dc_work(ind+j)=dc(j,i)
403 d_t(j,i+nres)=d_t_work(ind+j)
404 dc(j,i+nres)=zapas(ind+j)
405 dc_work(ind+j)=dc(j,i+nres)
411 write (iout,*) "Before correction for rotational lengthening"
412 write (iout,*) "New coordinates",&
413 " and differences between actual and standard bond lengths"
418 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
422 if (itype(i,1).ne.10) then
424 xx=vbld(i+nres)-vbldsc0(1,itype(i,1))
425 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
426 i,(dC(j,i+nres),j=1,3),xx
430 ! Second correction (rotational lengthening)
436 blen2 = scalar(dc(1,i),dc(1,i))
437 ppvec(ind)=2*vbl**2-blen2
438 diffbond=dabs(vbl-dsqrt(blen2))
439 if (diffbond.gt.diffmax) diffmax=diffbond
440 if (ppvec(ind).gt.0.0d0) then
441 ppvec(ind)=dsqrt(ppvec(ind))
446 write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind)
450 if (itype(i,1).ne.10) then
452 blen2 = scalar(dc(1,i+nres),dc(1,i+nres))
453 ppvec(ind)=2*vbldsc0(1,itype(i,1))**2-blen2
454 diffbond=dabs(vbldsc0(1,itype(i,1))-dsqrt(blen2))
455 if (diffbond.gt.diffmax) diffmax=diffbond
456 if (ppvec(ind).gt.0.0d0) then
457 ppvec(ind)=dsqrt(ppvec(ind))
462 write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind)
466 if (lprn) write (iout,*) "iter",iter," diffmax",diffmax
467 if (diffmax.lt.difftol) goto 10
471 Td(i)=Td(i)+ppvec(j)*Tmat(i,j)
477 zapas(i)=zapas(i)+Pmat(i,j)*dc_work(j)
488 dc_work(ind+j)=zapas(ind+j)
493 if (itype(i,1).ne.10) then
495 dc(j,i+nres)=zapas(ind+j)
496 dc_work(ind+j)=zapas(ind+j)
501 ! Building the chain from the newly calculated coordinates
504 if (large.and. mod(itime,ntwe).eq.0) then
505 write (iout,*) "Cartesian and internal coordinates: step 1"
508 write (iout,'(a)') "Potential forces"
510 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(-gcart(j,i),j=1,3),&
513 write (iout,'(a)') "Stochastic forces"
515 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(stochforc(j,i),j=1,3),&
516 (stochforc(j,i+nres),j=1,3)
518 write (iout,'(a)') "Velocities"
520 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
521 (d_t(j,i+nres),j=1,3)
526 write (iout,*) "After correction for rotational lengthening"
527 write (iout,*) "New coordinates",&
528 " and differences between actual and standard bond lengths"
533 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
537 if (itype(i,1).ne.10) then
539 xx=vbld(i+nres)-vbldsc0(1,itype(i,1))
540 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
541 i,(dC(j,i+nres),j=1,3),xx
546 ! write (iout,*) "Too many attempts at correcting the bonds"
554 ! Calculate energy and forces
556 call etotal(potEcomp)
557 potE=potEcomp(0)-potEcomp(20)
561 ! Calculate the kinetic and total energy and the kinetic temperature
564 t_enegrad=t_enegrad+MPI_Wtime()-tt0
566 t_enegrad=t_enegrad+tcpu()-tt0
569 kinetic_T=2.0d0/(dimen*Rb)*EK
571 end subroutine brown_step
572 !-----------------------------------------------------------------------------
574 !-----------------------------------------------------------------------------
575 subroutine gauss(RO,AP,MT,M,N,*)
577 ! CALCULATES (RO**(-1))*AP BY GAUSS ELIMINATION
578 ! RO IS A SQUARE MATRIX
579 ! THE CALCULATED PRODUCT IS STORED IN AP
580 ! ABNORMAL EXIT IF RO IS SINGULAR
582 integer :: MT, M, N, M1,I,J,IM,&
584 real(kind=8) :: RO(MT,M),AP(MT,N),X,RM,PR,Y
590 if(dabs(X).le.1.0D-13) return 1
602 if(DABS(RO(J,I)).LE.RM) goto 2
616 if(dabs(X).le.1.0E-13) return 1
625 8 AP(J,K)=AP(J,K)-Y*AP(I,K)
627 9 RO(J,K)=RO(J,K)-Y*RO(I,K)
631 if(dabs(X).le.1.0E-13) return 1
641 15 X=X-AP(K,J)*RO(MI,K)
646 !-----------------------------------------------------------------------------
648 !-----------------------------------------------------------------------------
649 subroutine kinetic(KE_total)
650 !----------------------------------------------------------------
651 ! This subroutine calculates the total kinetic energy of the chain
652 !-----------------------------------------------------------------
654 ! implicit real*8 (a-h,o-z)
655 ! include 'DIMENSIONS'
656 ! include 'COMMON.VAR'
657 ! include 'COMMON.CHAIN'
658 ! include 'COMMON.DERIV'
659 ! include 'COMMON.GEO'
660 ! include 'COMMON.LOCAL'
661 ! include 'COMMON.INTERACT'
662 ! include 'COMMON.MD'
663 ! include 'COMMON.IOUNITS'
664 real(kind=8) :: KE_total,mscab
666 integer :: i,j,k,iti,mnum
667 real(kind=8) :: KEt_p,KEt_sc,KEr_p,KEr_sc,incr(3),&
670 write (iout,*) "Velocities, kietic"
672 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
673 (d_t(j,i+nres),j=1,3)
678 ! write (iout,*) "ISC",(isc(itype(i,1)),i=1,nres)
679 ! The translational part for peptide virtual bonds
685 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
686 ! write (iout,*) "Kinetic trp:",i,(incr(j),j=1,3)
688 v(j)=incr(j)+0.5d0*d_t(j,i)
690 vtot(i)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
691 KEt_p=KEt_p+mp(mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
693 incr(j)=incr(j)+d_t(j,i)
696 ! write(iout,*) 'KEt_p', KEt_p
697 ! The translational part for the side chain virtual bond
698 ! Only now we can initialize incr with zeros. It must be equal
699 ! to the velocities of the first Calpha.
705 iti=iabs(itype(i,mnum))
711 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
712 if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
713 .or.(mnum.eq.5)) then
719 v(j)=incr(j)+d_t(j,nres+i)
722 ! write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3)
723 ! write (iout,*) "i",i," msc",msc(iti)," v",(v(j),j=1,3)
724 KEt_sc=KEt_sc+mscab*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
725 vtot(i+nres)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
727 incr(j)=incr(j)+d_t(j,i)
731 ! write(iout,*) 'KEt_sc', KEt_sc
732 ! The part due to stretching and rotation of the peptide groups
736 ! write (iout,*) "i",i
737 ! write (iout,*) "i",i," mag1",mag1," mag2",mag2
741 ! write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3)
742 KEr_p=KEr_p+Ip(mnum)*(incr(1)*incr(1)+incr(2)*incr(2) &
746 ! write(iout,*) 'KEr_p', KEr_p
747 ! The rotational part of the side chain virtual bond
751 iti=iabs(itype(i,mnum))
752 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
753 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
754 .and.(mnum.ne.5)) then
756 incr(j)=d_t(j,nres+i)
758 ! write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3)
759 KEr_sc=KEr_sc+Isc(iti,mnum)*(incr(1)*incr(1)+incr(2)*incr(2)+ &
763 ! The total kinetic energy
765 ! write(iout,*) 'KEr_sc', KEr_sc
766 KE_total=0.5d0*(KEt_p+KEt_sc+0.25d0*KEr_p+KEr_sc)
767 ! write (iout,*) "KE_total",KE_total
769 end subroutine kinetic
770 !-----------------------------------------------------------------------------
772 !-----------------------------------------------------------------------------
774 !------------------------------------------------
775 ! The driver for molecular dynamics subroutines
776 !------------------------------------------------
779 use control, only:tcpu,ovrtim
780 ! use io_comm, only:ilen
782 use compare, only:secondary2,hairpin
783 use io, only:cartout,statout
784 ! implicit real*8 (a-h,o-z)
785 ! include 'DIMENSIONS'
788 integer :: IERROR,ERRCODE
790 ! include 'COMMON.SETUP'
791 ! include 'COMMON.CONTROL'
792 ! include 'COMMON.VAR'
793 ! include 'COMMON.MD'
795 ! include 'COMMON.LANGEVIN'
797 ! include 'COMMON.LANGEVIN.lang0'
799 ! include 'COMMON.CHAIN'
800 ! include 'COMMON.DERIV'
801 ! include 'COMMON.GEO'
802 ! include 'COMMON.LOCAL'
803 ! include 'COMMON.INTERACT'
804 ! include 'COMMON.IOUNITS'
805 ! include 'COMMON.NAMES'
806 ! include 'COMMON.TIME1'
807 ! include 'COMMON.HAIRPIN'
808 real(kind=8),dimension(3) :: L,vcm
810 real(kind=8),dimension(6*nres) :: v_work,v_transf !(maxres6) maxres6=6*maxres
812 integer :: rstcount !ilen,
814 character(len=50) :: tytul
815 !el common /gucio/ cm
816 integer :: itime,i,j,nharp
817 integer,dimension(4,nres/3) :: iharp !(4,nres/3)(4,maxres/3)
819 real(kind=8) :: tt0,scalfac
825 print *,"MY tmpdir",tmpdir,ilen(tmpdir)
826 if (ilen(tmpdir).gt.0) &
827 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_" &
828 //liczba(:ilen(liczba))//'.rst')
830 if (ilen(tmpdir).gt.0) &
831 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"//'.rst')
838 write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
844 print *,"just befor setup matix",nres
845 ! Determine the inverse of the inertia matrix.
846 call setup_MD_matrices
848 print *,"AFTER SETUP MATRICES"
850 print *,"AFTER INIT MD"
853 t_MDsetup = MPI_Wtime()-tt0
855 t_MDsetup = tcpu()-tt0
858 ! Entering the MD loop
864 if (lang.eq.2 .or. lang.eq.3) then
868 call sd_verlet_p_setup
870 call sd_verlet_ciccotti_setup
874 pfric0_mat(i,j,0)=pfric_mat(i,j)
875 afric0_mat(i,j,0)=afric_mat(i,j)
876 vfric0_mat(i,j,0)=vfric_mat(i,j)
877 prand0_mat(i,j,0)=prand_mat(i,j)
878 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
879 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
884 flag_stoch(i)=.false.
888 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
890 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
894 else if (lang.eq.1 .or. lang.eq.4) then
895 print *,"before setup_fricmat"
897 print *,"after setup_fricmat"
900 t_langsetup=MPI_Wtime()-tt0
903 t_langsetup=tcpu()-tt0
906 do itime=1,n_timestep
908 if (large.and. mod(itime,ntwe).eq.0) &
909 write (iout,*) "itime",itime
911 if (lang.gt.0 .and. surfarea .and. &
912 mod(itime,reset_fricmat).eq.0) then
913 if (lang.eq.2 .or. lang.eq.3) then
917 call sd_verlet_p_setup
919 call sd_verlet_ciccotti_setup
923 pfric0_mat(i,j,0)=pfric_mat(i,j)
924 afric0_mat(i,j,0)=afric_mat(i,j)
925 vfric0_mat(i,j,0)=vfric_mat(i,j)
926 prand0_mat(i,j,0)=prand_mat(i,j)
927 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
928 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
933 flag_stoch(i)=.false.
936 else if (lang.eq.1 .or. lang.eq.4) then
939 write (iout,'(a,i10)') &
940 "Friction matrix reset based on surface area, itime",itime
942 if (reset_vel .and. tbf .and. lang.eq.0 &
943 .and. mod(itime,count_reset_vel).eq.0) then
945 write(iout,'(a,f20.2)') &
946 "Velocities reset to random values, time",totT
949 d_t_old(j,i)=d_t(j,i)
953 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
957 d_t(j,0)=d_t(j,0)-vcm(j)
960 kinetic_T=2.0d0/(dimen3*Rb)*EK
961 scalfac=dsqrt(T_bath/kinetic_T)
962 write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
965 d_t_old(j,i)=scalfac*d_t(j,i)
971 ! Time-reversible RESPA algorithm
972 ! (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
973 call RESPA_step(itime)
975 ! Variable time step algorithm.
976 call velverlet_step(itime)
980 call brown_step(itime)
982 print *,"Brown dynamics not here!"
984 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
990 if (mod(itime,ntwe).eq.0) then
993 ! call check_ecartint
1003 v_work(ind)=d_t(j,i)
1008 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1.and.mnum.ne.5) then
1011 v_work(ind)=d_t(j,i+nres)
1016 write (66,'(80f10.5)') &
1017 ((d_t(j,i),j=1,3),i=0,nres-1),((d_t(j,i+nres),j=1,3),i=1,nres)
1021 v_transf(i)=v_transf(i)+gvec(j,i)*v_work(j)
1023 v_transf(i)= v_transf(i)*dsqrt(geigen(i))
1025 write (67,'(80f10.5)') (v_transf(i),i=1,ind)
1028 if (mod(itime,ntwx).eq.0) then
1030 write (tytul,'("time",f8.2)') totT
1032 call hairpin(.true.,nharp,iharp)
1033 call secondary2(.true.)
1034 call pdbout(potE,tytul,ipdb)
1039 if (rstcount.eq.1000.or.itime.eq.n_timestep) then
1040 open(irest2,file=rest2name,status='unknown')
1041 write(irest2,*) totT,EK,potE,totE,t_bath
1043 ! AL 4/17/17: Now writing d_t(0,:) too
1045 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
1047 ! AL 4/17/17: Now writing d_c(0,:) too
1049 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
1057 t_MD=MPI_Wtime()-tt0
1061 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))') &
1063 'MD calculations setup:',t_MDsetup,&
1064 'Energy & gradient evaluation:',t_enegrad,&
1065 'Stochastic MD setup:',t_langsetup,&
1066 'Stochastic MD step setup:',t_sdsetup,&
1068 write (iout,'(/28(1h=),a25,27(1h=))') &
1069 ' End of MD calculation '
1071 write (iout,*) "time for etotal",t_etotal," elong",t_elong,&
1073 write (iout,*) "time_fric",time_fric," time_stoch",time_stoch,&
1074 " time_fricmatmult",time_fricmatmult," time_fsample ",&
1079 !-----------------------------------------------------------------------------
1080 subroutine velverlet_step(itime)
1081 !-------------------------------------------------------------------------------
1082 ! Perform a single velocity Verlet step; the time step can be rescaled if
1083 ! increments in accelerations exceed the threshold
1084 !-------------------------------------------------------------------------------
1085 ! implicit real*8 (a-h,o-z)
1086 ! include 'DIMENSIONS'
1088 use control, only:tcpu
1092 integer :: ierror,ierrcode
1093 real(kind=8) :: errcode
1095 ! include 'COMMON.SETUP'
1096 ! include 'COMMON.VAR'
1097 ! include 'COMMON.MD'
1099 ! include 'COMMON.LANGEVIN'
1101 ! include 'COMMON.LANGEVIN.lang0'
1103 ! include 'COMMON.CHAIN'
1104 ! include 'COMMON.DERIV'
1105 ! include 'COMMON.GEO'
1106 ! include 'COMMON.LOCAL'
1107 ! include 'COMMON.INTERACT'
1108 ! include 'COMMON.IOUNITS'
1109 ! include 'COMMON.NAMES'
1110 ! include 'COMMON.TIME1'
1111 ! include 'COMMON.MUCA'
1112 real(kind=8),dimension(3) :: vcm,incr
1113 real(kind=8),dimension(3) :: L
1114 integer :: count,rstcount !ilen,
1116 character(len=50) :: tytul
1117 integer :: maxcount_scale = 20
1118 !el common /gucio/ cm
1119 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1120 !el common /stochcalc/ stochforcvec
1121 integer :: itime,icount_scale,itime_scal,i,j,ifac_time
1123 real(kind=8) :: epdrift,tt0,fac_time
1125 if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres)) !(MAXRES6) maxres6=6*maxres
1131 else if (lang.eq.2 .or. lang.eq.3) then
1133 call stochastic_force(stochforcvec)
1136 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
1138 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1145 icount_scale=icount_scale+1
1146 if (icount_scale.gt.maxcount_scale) then
1148 "ERROR: too many attempts at scaling down the time step. ",&
1149 "amax=",amax,"epdrift=",epdrift,&
1150 "damax=",damax,"edriftmax=",edriftmax,&
1154 call MPI_Abort(MPI_COMM_WORLD,IERROR,IERRCODE)
1158 ! First step of the velocity Verlet algorithm
1163 else if (lang.eq.3) then
1165 call sd_verlet1_ciccotti
1167 else if (lang.eq.1) then
1172 ! Build the chain from the newly calculated coordinates
1173 call chainbuild_cart
1174 if (rattle) call rattle1
1176 if (large.and. mod(itime,ntwe).eq.0) then
1177 write (iout,*) "Cartesian and internal coordinates: step 1"
1182 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),&
1183 (dc(j,i+nres),j=1,3)
1185 write (iout,*) "Accelerations"
1187 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1188 (d_a(j,i+nres),j=1,3)
1190 write (iout,*) "Velocities, step 1"
1192 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1193 (d_t(j,i+nres),j=1,3)
1202 ! Calculate energy and forces
1204 call etotal(potEcomp)
1205 ! AL 4/17/17: Reduce the steps if NaNs occurred.
1206 if (potEcomp(0).gt.0.99e18 .or. isnan(potEcomp(0)).gt.0) then
1211 if (large.and. mod(itime,ntwe).eq.0) &
1212 call enerprint(potEcomp)
1215 t_etotal=t_etotal+MPI_Wtime()-tt0
1217 t_etotal=t_etotal+tcpu()-tt0
1220 potE=potEcomp(0)-potEcomp(20)
1222 ! Get the new accelerations
1225 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1227 t_enegrad=t_enegrad+tcpu()-tt0
1229 ! Determine maximum acceleration and scale down the timestep if needed
1231 amax=amax/(itime_scal+1)**2
1232 call predict_edrift(epdrift)
1233 if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then
1234 ! Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
1236 ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax)) &
1238 itime_scal=itime_scal+ifac_time
1239 ! fac_time=dmin1(damax/amax,0.5d0)
1240 fac_time=0.5d0**ifac_time
1241 d_time=d_time*fac_time
1242 if (lang.eq.2 .or. lang.eq.3) then
1244 ! write (iout,*) "Calling sd_verlet_setup: 1"
1245 ! Rescale the stochastic forces and recalculate or restore
1246 ! the matrices of tinker integrator
1247 if (itime_scal.gt.maxflag_stoch) then
1248 if (large) write (iout,'(a,i5,a)') &
1249 "Calculate matrices for stochastic step;",&
1250 " itime_scal ",itime_scal
1252 call sd_verlet_p_setup
1254 call sd_verlet_ciccotti_setup
1256 write (iout,'(2a,i3,a,i3,1h.)') &
1257 "Warning: cannot store matrices for stochastic",&
1258 " integration because the index",itime_scal,&
1259 " is greater than",maxflag_stoch
1260 write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",&
1261 " integration Langevin algorithm for better efficiency."
1262 else if (flag_stoch(itime_scal)) then
1263 if (large) write (iout,'(a,i5,a,l1)') &
1264 "Restore matrices for stochastic step; itime_scal ",&
1265 itime_scal," flag ",flag_stoch(itime_scal)
1268 pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
1269 afric_mat(i,j)=afric0_mat(i,j,itime_scal)
1270 vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
1271 prand_mat(i,j)=prand0_mat(i,j,itime_scal)
1272 vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
1273 vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
1277 if (large) write (iout,'(2a,i5,a,l1)') &
1278 "Calculate & store matrices for stochastic step;",&
1279 " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
1281 call sd_verlet_p_setup
1283 call sd_verlet_ciccotti_setup
1285 flag_stoch(ifac_time)=.true.
1288 pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
1289 afric0_mat(i,j,itime_scal)=afric_mat(i,j)
1290 vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
1291 prand0_mat(i,j,itime_scal)=prand_mat(i,j)
1292 vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
1293 vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
1297 fac_time=1.0d0/dsqrt(fac_time)
1299 stochforcvec(i)=fac_time*stochforcvec(i)
1302 else if (lang.eq.1) then
1303 ! Rescale the accelerations due to stochastic forces
1304 fac_time=1.0d0/dsqrt(fac_time)
1306 d_as_work(i)=d_as_work(i)*fac_time
1309 if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)') &
1310 "itime",itime," Timestep scaled down to ",&
1311 d_time," ifac_time",ifac_time," itime_scal",itime_scal
1313 ! Second step of the velocity Verlet algorithm
1318 else if (lang.eq.3) then
1320 call sd_verlet2_ciccotti
1322 else if (lang.eq.1) then
1327 if (rattle) call rattle2
1330 if (d_time.ne.d_time0) then
1333 if (lang.eq.2 .or. lang.eq.3) then
1334 if (large) write (iout,'(a)') &
1335 "Restore original matrices for stochastic step"
1336 ! write (iout,*) "Calling sd_verlet_setup: 2"
1337 ! Restore the matrices of tinker integrator if the time step has been restored
1340 pfric_mat(i,j)=pfric0_mat(i,j,0)
1341 afric_mat(i,j)=afric0_mat(i,j,0)
1342 vfric_mat(i,j)=vfric0_mat(i,j,0)
1343 prand_mat(i,j)=prand0_mat(i,j,0)
1344 vrand_mat1(i,j)=vrand0_mat1(i,j,0)
1345 vrand_mat2(i,j)=vrand0_mat2(i,j,0)
1354 ! Calculate the kinetic and the total energy and the kinetic temperature
1358 ! call kinetic1(EK1)
1359 ! write (iout,*) "step",itime," EK",EK," EK1",EK1
1361 ! Couple the system to Berendsen bath if needed
1362 if (tbf .and. lang.eq.0) then
1365 kinetic_T=2.0d0/(dimen3*Rb)*EK
1366 ! Backup the coordinates, velocities, and accelerations
1370 d_t_old(j,i)=d_t(j,i)
1371 d_a_old(j,i)=d_a(j,i)
1375 if (mod(itime,ntwe).eq.0 .and. large) then
1376 write (iout,*) "Velocities, step 2"
1378 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1379 (d_t(j,i+nres),j=1,3)
1384 end subroutine velverlet_step
1385 !-----------------------------------------------------------------------------
1386 subroutine RESPA_step(itime)
1387 !-------------------------------------------------------------------------------
1388 ! Perform a single RESPA step.
1389 !-------------------------------------------------------------------------------
1390 ! implicit real*8 (a-h,o-z)
1391 ! include 'DIMENSIONS'
1395 use control, only:tcpu
1397 ! use io_conf, only:cartprint
1400 integer :: IERROR,ERRCODE
1402 ! include 'COMMON.SETUP'
1403 ! include 'COMMON.CONTROL'
1404 ! include 'COMMON.VAR'
1405 ! include 'COMMON.MD'
1407 ! include 'COMMON.LANGEVIN'
1409 ! include 'COMMON.LANGEVIN.lang0'
1411 ! include 'COMMON.CHAIN'
1412 ! include 'COMMON.DERIV'
1413 ! include 'COMMON.GEO'
1414 ! include 'COMMON.LOCAL'
1415 ! include 'COMMON.INTERACT'
1416 ! include 'COMMON.IOUNITS'
1417 ! include 'COMMON.NAMES'
1418 ! include 'COMMON.TIME1'
1419 real(kind=8),dimension(0:n_ene) :: energia_short,energia_long
1420 real(kind=8),dimension(3) :: L,vcm,incr
1421 real(kind=8),dimension(3,0:2*nres) :: dc_old0,d_t_old0,d_a_old0 !(3,0:maxres2) maxres2=2*maxres
1422 logical :: PRINT_AMTS_MSG = .false.
1423 integer :: count,rstcount !ilen,
1425 character(len=50) :: tytul
1426 integer :: maxcount_scale = 10
1427 !el common /gucio/ cm
1428 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1429 !el common /stochcalc/ stochforcvec
1430 integer :: itime,itt,i,j,itsplit
1432 !el common /cipiszcze/ itt
1434 real(kind=8) :: epdrift,tt0,epdriftmax
1437 if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres)) !(MAXRES6) maxres6=6*maxres
1441 if (large.and. mod(itime,ntwe).eq.0) then
1442 write (iout,*) "***************** RESPA itime",itime
1443 write (iout,*) "Cartesian and internal coordinates: step 0"
1445 call pdbout(0.0d0,"cipiszcze",iout)
1447 write (iout,*) "Accelerations from long-range forces"
1449 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1450 (d_a(j,i+nres),j=1,3)
1452 write (iout,*) "Velocities, step 0"
1454 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1455 (d_t(j,i+nres),j=1,3)
1460 ! Perform the initial RESPA step (increment velocities)
1461 ! write (iout,*) "*********************** RESPA ini"
1464 if (mod(itime,ntwe).eq.0 .and. large) then
1465 write (iout,*) "Velocities, end"
1467 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1468 (d_t(j,i+nres),j=1,3)
1472 ! Compute the short-range forces
1478 ! 7/2/2009 commented out
1480 ! call etotal_short(energia_short)
1483 ! 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
1486 d_a(j,i)=d_a_short(j,i)
1490 if (large.and. mod(itime,ntwe).eq.0) then
1491 write (iout,*) "energia_short",energia_short(0)
1492 write (iout,*) "Accelerations from short-range forces"
1494 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1495 (d_a(j,i+nres),j=1,3)
1500 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1502 t_enegrad=t_enegrad+tcpu()-tt0
1507 d_t_old(j,i)=d_t(j,i)
1508 d_a_old(j,i)=d_a(j,i)
1511 ! 6/30/08 A-MTS: attempt at increasing the split number
1514 dc_old0(j,i)=dc_old(j,i)
1515 d_t_old0(j,i)=d_t_old(j,i)
1516 d_a_old0(j,i)=d_a_old(j,i)
1519 if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
1520 if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
1527 ! write (iout,*) "itime",itime," ntime_split",ntime_split
1528 ! Split the time step
1529 d_time=d_time0/ntime_split
1530 ! Perform the short-range RESPA steps (velocity Verlet increments of
1531 ! positions and velocities using short-range forces)
1532 ! write (iout,*) "*********************** RESPA split"
1533 do itsplit=1,ntime_split
1536 else if (lang.eq.2 .or. lang.eq.3) then
1538 call stochastic_force(stochforcvec)
1541 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
1543 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1548 ! First step of the velocity Verlet algorithm
1553 else if (lang.eq.3) then
1555 call sd_verlet1_ciccotti
1557 else if (lang.eq.1) then
1562 ! Build the chain from the newly calculated coordinates
1563 call chainbuild_cart
1564 if (rattle) call rattle1
1566 if (large.and. mod(itime,ntwe).eq.0) then
1567 write (iout,*) "***** ITSPLIT",itsplit
1568 write (iout,*) "Cartesian and internal coordinates: step 1"
1569 call pdbout(0.0d0,"cipiszcze",iout)
1572 write (iout,*) "Velocities, step 1"
1574 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1575 (d_t(j,i+nres),j=1,3)
1584 ! Calculate energy and forces
1586 call etotal_short(energia_short)
1587 ! AL 4/17/17: Exit itime_split loop when energy goes infinite
1588 if (energia_short(0).gt.0.99e20 .or. isnan(energia_short(0)) ) then
1589 if (PRINT_AMTS_MSG) &
1590 write (iout,*) "Infinities/NaNs in energia_short",energia_short(0),"; increasing ntime_split to",ntime_split
1591 ntime_split=ntime_split*2
1592 if (ntime_split.gt.maxtime_split) then
1595 "Cannot rescue the run; aborting job. Retry with a smaller time step"
1597 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1600 "Cannot rescue the run; terminating. Retry with a smaller time step"
1606 if (large.and. mod(itime,ntwe).eq.0) &
1607 call enerprint(energia_short)
1610 t_eshort=t_eshort+MPI_Wtime()-tt0
1612 t_eshort=t_eshort+tcpu()-tt0
1616 ! Get the new accelerations
1618 ! 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
1621 d_a_short(j,i)=d_a(j,i)
1625 if (large.and. mod(itime,ntwe).eq.0) then
1626 write (iout,*)"energia_short",energia_short(0)
1627 write (iout,*) "Accelerations from short-range forces"
1629 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1630 (d_a(j,i+nres),j=1,3)
1635 ! Determine maximum acceleration and scale down the timestep if needed
1637 amax=amax/ntime_split**2
1638 call predict_edrift(epdrift)
1639 if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0) &
1640 write (iout,*) "amax",amax," damax",damax,&
1641 " epdrift",epdrift," epdriftmax",epdriftmax
1642 ! Exit loop and try with increased split number if the change of
1643 ! acceleration is too big
1644 if (amax.gt.damax .or. epdrift.gt.edriftmax) then
1645 if (ntime_split.lt.maxtime_split) then
1647 ntime_split=ntime_split*2
1648 ! AL 4/17/17: We should exit the itime_split loop when acceleration change is too big
1652 dc_old(j,i)=dc_old0(j,i)
1653 d_t_old(j,i)=d_t_old0(j,i)
1654 d_a_old(j,i)=d_a_old0(j,i)
1657 if (PRINT_AMTS_MSG) then
1658 write (iout,*) "acceleration/energy drift too large",amax,&
1659 epdrift," split increased to ",ntime_split," itime",itime,&
1665 "Uh-hu. Bumpy landscape. Maximum splitting number",&
1667 " already reached!!! Trying to carry on!"
1671 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1673 t_enegrad=t_enegrad+tcpu()-tt0
1675 ! Second step of the velocity Verlet algorithm
1680 else if (lang.eq.3) then
1682 call sd_verlet2_ciccotti
1684 else if (lang.eq.1) then
1689 if (rattle) call rattle2
1690 ! Backup the coordinates, velocities, and accelerations
1694 d_t_old(j,i)=d_t(j,i)
1695 d_a_old(j,i)=d_a(j,i)
1702 ! Restore the time step
1704 ! Compute long-range forces
1711 call etotal_long(energia_long)
1712 if (energia_long(0).gt.0.99e20 .or. isnan(energia_long(0))) then
1715 "Infinitied/NaNs in energia_long, Aborting MPI job."
1717 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1719 write (iout,*) "Infinitied/NaNs in energia_long, terminating."
1723 if (large.and. mod(itime,ntwe).eq.0) &
1724 call enerprint(energia_long)
1727 t_elong=t_elong+MPI_Wtime()-tt0
1729 t_elong=t_elong+tcpu()-tt0
1735 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1737 t_enegrad=t_enegrad+tcpu()-tt0
1739 ! Compute accelerations from long-range forces
1741 if (large.and. mod(itime,ntwe).eq.0) then
1742 write (iout,*) "energia_long",energia_long(0)
1743 write (iout,*) "Cartesian and internal coordinates: step 2"
1745 call pdbout(0.0d0,"cipiszcze",iout)
1747 write (iout,*) "Accelerations from long-range forces"
1749 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1750 (d_a(j,i+nres),j=1,3)
1752 write (iout,*) "Velocities, step 2"
1754 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1755 (d_t(j,i+nres),j=1,3)
1759 ! Compute the final RESPA step (increment velocities)
1760 ! write (iout,*) "*********************** RESPA fin"
1762 ! Compute the complete potential energy
1764 potEcomp(i)=energia_short(i)+energia_long(i)
1766 potE=potEcomp(0)-potEcomp(20)
1767 ! potE=energia_short(0)+energia_long(0)
1770 ! Calculate the kinetic and the total energy and the kinetic temperature
1773 ! Couple the system to Berendsen bath if needed
1774 if (tbf .and. lang.eq.0) then
1777 kinetic_T=2.0d0/(dimen3*Rb)*EK
1778 ! Backup the coordinates, velocities, and accelerations
1780 if (mod(itime,ntwe).eq.0 .and. large) then
1781 write (iout,*) "Velocities, end"
1783 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1784 (d_t(j,i+nres),j=1,3)
1789 end subroutine RESPA_step
1790 !-----------------------------------------------------------------------------
1791 subroutine RESPA_vel
1792 ! First and last RESPA step (incrementing velocities using long-range
1795 ! implicit real*8 (a-h,o-z)
1796 ! include 'DIMENSIONS'
1797 ! include 'COMMON.CONTROL'
1798 ! include 'COMMON.VAR'
1799 ! include 'COMMON.MD'
1800 ! include 'COMMON.CHAIN'
1801 ! include 'COMMON.DERIV'
1802 ! include 'COMMON.GEO'
1803 ! include 'COMMON.LOCAL'
1804 ! include 'COMMON.INTERACT'
1805 ! include 'COMMON.IOUNITS'
1806 ! include 'COMMON.NAMES'
1807 integer :: i,j,inres,mnum
1810 d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
1814 d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
1819 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
1820 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
1821 .and.(mnum.ne.5)) then
1824 d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
1829 end subroutine RESPA_vel
1830 !-----------------------------------------------------------------------------
1832 ! Applying velocity Verlet algorithm - step 1 to coordinates
1834 ! implicit real*8 (a-h,o-z)
1835 ! include 'DIMENSIONS'
1836 ! include 'COMMON.CONTROL'
1837 ! include 'COMMON.VAR'
1838 ! include 'COMMON.MD'
1839 ! include 'COMMON.CHAIN'
1840 ! include 'COMMON.DERIV'
1841 ! include 'COMMON.GEO'
1842 ! include 'COMMON.LOCAL'
1843 ! include 'COMMON.INTERACT'
1844 ! include 'COMMON.IOUNITS'
1845 ! include 'COMMON.NAMES'
1846 real(kind=8) :: adt,adt2
1847 integer :: i,j,inres,mnum
1850 write (iout,*) "VELVERLET1 START: DC"
1852 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),&
1853 (dc(j,i+nres),j=1,3)
1857 adt=d_a_old(j,0)*d_time
1859 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1860 d_t_new(j,0)=d_t_old(j,0)+adt2
1861 d_t(j,0)=d_t_old(j,0)+adt
1865 adt=d_a_old(j,i)*d_time
1867 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1868 d_t_new(j,i)=d_t_old(j,i)+adt2
1869 d_t(j,i)=d_t_old(j,i)+adt
1874 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
1875 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
1876 .and.(mnum.ne.5)) then
1879 adt=d_a_old(j,inres)*d_time
1881 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1882 d_t_new(j,inres)=d_t_old(j,inres)+adt2
1883 d_t(j,inres)=d_t_old(j,inres)+adt
1888 write (iout,*) "VELVERLET1 END: DC"
1890 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),&
1891 (dc(j,i+nres),j=1,3)
1895 end subroutine verlet1
1896 !-----------------------------------------------------------------------------
1898 ! Step 2 of the velocity Verlet algorithm: update velocities
1900 ! implicit real*8 (a-h,o-z)
1901 ! include 'DIMENSIONS'
1902 ! include 'COMMON.CONTROL'
1903 ! include 'COMMON.VAR'
1904 ! include 'COMMON.MD'
1905 ! include 'COMMON.CHAIN'
1906 ! include 'COMMON.DERIV'
1907 ! include 'COMMON.GEO'
1908 ! include 'COMMON.LOCAL'
1909 ! include 'COMMON.INTERACT'
1910 ! include 'COMMON.IOUNITS'
1911 ! include 'COMMON.NAMES'
1912 integer :: i,j,inres,mnum
1915 d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
1919 d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
1924 ! iti=iabs(itype(i,mnum))
1925 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
1926 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
1927 .and.(mnum.ne.5)) then
1930 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
1935 end subroutine verlet2
1936 !-----------------------------------------------------------------------------
1937 subroutine sddir_precalc
1938 ! Applying velocity Verlet algorithm - step 1 to coordinates
1939 ! implicit real*8 (a-h,o-z)
1940 ! include 'DIMENSIONS'
1946 ! include 'COMMON.CONTROL'
1947 ! include 'COMMON.VAR'
1948 ! include 'COMMON.MD'
1950 ! include 'COMMON.LANGEVIN'
1952 ! include 'COMMON.LANGEVIN.lang0'
1954 ! include 'COMMON.CHAIN'
1955 ! include 'COMMON.DERIV'
1956 ! include 'COMMON.GEO'
1957 ! include 'COMMON.LOCAL'
1958 ! include 'COMMON.INTERACT'
1959 ! include 'COMMON.IOUNITS'
1960 ! include 'COMMON.NAMES'
1961 ! include 'COMMON.TIME1'
1962 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1963 !el common /stochcalc/ stochforcvec
1964 real(kind=8) :: time00
1966 ! Compute friction and stochastic forces
1971 time_fric=time_fric+MPI_Wtime()-time00
1973 call stochastic_force(stochforcvec)
1974 time_stoch=time_stoch+MPI_Wtime()-time00
1977 ! Compute the acceleration due to friction forces (d_af_work) and stochastic
1978 ! forces (d_as_work)
1980 call ginv_mult(fric_work, d_af_work)
1981 call ginv_mult(stochforcvec, d_as_work)
1983 end subroutine sddir_precalc
1984 !-----------------------------------------------------------------------------
1985 subroutine sddir_verlet1
1986 ! Applying velocity Verlet algorithm - step 1 to velocities
1989 ! implicit real*8 (a-h,o-z)
1990 ! include 'DIMENSIONS'
1991 ! include 'COMMON.CONTROL'
1992 ! include 'COMMON.VAR'
1993 ! include 'COMMON.MD'
1995 ! include 'COMMON.LANGEVIN'
1997 ! include 'COMMON.LANGEVIN.lang0'
1999 ! include 'COMMON.CHAIN'
2000 ! include 'COMMON.DERIV'
2001 ! include 'COMMON.GEO'
2002 ! include 'COMMON.LOCAL'
2003 ! include 'COMMON.INTERACT'
2004 ! include 'COMMON.IOUNITS'
2005 ! include 'COMMON.NAMES'
2006 ! Revised 3/31/05 AL: correlation between random contributions to
2007 ! position and velocity increments included.
2008 real(kind=8) :: sqrt13 = 0.57735026918962576451d0 ! 1/sqrt(3)
2009 real(kind=8) :: adt,adt2
2010 integer :: i,j,ind,inres,mnum
2012 ! Add the contribution from BOTH friction and stochastic force to the
2013 ! coordinates, but ONLY the contribution from the friction forces to velocities
2016 adt=(d_a_old(j,0)+d_af_work(j))*d_time
2017 adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
2018 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
2019 d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
2020 d_t(j,0)=d_t_old(j,0)+adt
2025 adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
2026 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
2027 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
2028 d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
2029 d_t(j,i)=d_t_old(j,i)+adt
2035 ! iti=iabs(itype(i,mnum))
2036 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2037 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2038 .and.(mnum.ne.5)) then
2041 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
2042 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
2043 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
2044 d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
2045 d_t(j,inres)=d_t_old(j,inres)+adt
2051 end subroutine sddir_verlet1
2052 !-----------------------------------------------------------------------------
2053 subroutine sddir_verlet2
2054 ! Calculating the adjusted velocities for accelerations
2057 ! implicit real*8 (a-h,o-z)
2058 ! include 'DIMENSIONS'
2059 ! include 'COMMON.CONTROL'
2060 ! include 'COMMON.VAR'
2061 ! include 'COMMON.MD'
2063 ! include 'COMMON.LANGEVIN'
2065 ! include 'COMMON.LANGEVIN.lang0'
2067 ! include 'COMMON.CHAIN'
2068 ! include 'COMMON.DERIV'
2069 ! include 'COMMON.GEO'
2070 ! include 'COMMON.LOCAL'
2071 ! include 'COMMON.INTERACT'
2072 ! include 'COMMON.IOUNITS'
2073 ! include 'COMMON.NAMES'
2074 real(kind=8),dimension(6*nres) :: stochforcvec,d_as_work1 !(MAXRES6) maxres6=6*maxres
2075 real(kind=8) :: cos60 = 0.5d0, sin60 = 0.86602540378443864676d0
2076 integer :: i,j,ind,inres,mnum
2077 ! Revised 3/31/05 AL: correlation between random contributions to
2078 ! position and velocity increments included.
2079 ! The correlation coefficients are calculated at low-friction limit.
2080 ! Also, friction forces are now not calculated with new velocities.
2082 ! call friction_force
2083 call stochastic_force(stochforcvec)
2085 ! Compute the acceleration due to friction forces (d_af_work) and stochastic
2086 ! forces (d_as_work)
2088 call ginv_mult(stochforcvec, d_as_work1)
2094 d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j)) &
2095 +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
2100 d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j)) &
2101 +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
2107 ! iti=iabs(itype(i,mnum))
2108 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2109 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2110 .and.(mnum.ne.5)) then
2113 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres) &
2114 +d_af_work(ind+j))+sin60*d_as_work(ind+j) &
2115 +cos60*d_as_work1(ind+j))*d_time
2121 end subroutine sddir_verlet2
2122 !-----------------------------------------------------------------------------
2123 subroutine max_accel
2125 ! Find the maximum difference in the accelerations of the the sites
2126 ! at the beginning and the end of the time step.
2129 ! implicit real*8 (a-h,o-z)
2130 ! include 'DIMENSIONS'
2131 ! include 'COMMON.CONTROL'
2132 ! include 'COMMON.VAR'
2133 ! include 'COMMON.MD'
2134 ! include 'COMMON.CHAIN'
2135 ! include 'COMMON.DERIV'
2136 ! include 'COMMON.GEO'
2137 ! include 'COMMON.LOCAL'
2138 ! include 'COMMON.INTERACT'
2139 ! include 'COMMON.IOUNITS'
2140 real(kind=8),dimension(3) :: aux,accel,accel_old
2141 real(kind=8) :: dacc
2145 ! aux(j)=d_a(j,0)-d_a_old(j,0)
2146 accel_old(j)=d_a_old(j,0)
2153 ! 7/3/08 changed to asymmetric difference
2155 ! accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
2156 accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
2157 accel(j)=accel(j)+0.5d0*d_a(j,i)
2158 ! if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
2159 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
2160 dacc=dabs(accel(j)-accel_old(j))
2161 ! write (iout,*) i,dacc
2162 if (dacc.gt.amax) amax=dacc
2170 accel_old(j)=d_a_old(j,0)
2175 accel_old(j)=accel_old(j)+d_a_old(j,1)
2176 accel(j)=accel(j)+d_a(j,1)
2181 ! iti=iabs(itype(i,mnum))
2182 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2183 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2184 .and.(mnum.ne.5)) then
2186 ! accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
2187 accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
2188 accel(j)=accel(j)+d_a(j,i+nres)
2192 ! if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
2193 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
2194 dacc=dabs(accel(j)-accel_old(j))
2195 ! write (iout,*) "side-chain",i,dacc
2196 if (dacc.gt.amax) amax=dacc
2200 accel_old(j)=accel_old(j)+d_a_old(j,i)
2201 accel(j)=accel(j)+d_a(j,i)
2202 ! aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
2206 end subroutine max_accel
2207 !-----------------------------------------------------------------------------
2208 subroutine predict_edrift(epdrift)
2210 ! Predict the drift of the potential energy
2213 use control_data, only: lmuca
2214 ! implicit real*8 (a-h,o-z)
2215 ! include 'DIMENSIONS'
2216 ! include 'COMMON.CONTROL'
2217 ! include 'COMMON.VAR'
2218 ! include 'COMMON.MD'
2219 ! include 'COMMON.CHAIN'
2220 ! include 'COMMON.DERIV'
2221 ! include 'COMMON.GEO'
2222 ! include 'COMMON.LOCAL'
2223 ! include 'COMMON.INTERACT'
2224 ! include 'COMMON.IOUNITS'
2225 ! include 'COMMON.MUCA'
2226 real(kind=8) :: epdrift,epdriftij
2228 ! Drift of the potential energy
2234 epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
2235 if (lmuca) epdriftij=epdriftij*factor
2236 ! write (iout,*) "back",i,j,epdriftij
2237 if (epdriftij.gt.epdrift) epdrift=epdriftij
2241 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1.and.&
2242 molnum(i).ne.5) then
2245 dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
2246 if (lmuca) epdriftij=epdriftij*factor
2247 ! write (iout,*) "side",i,j,epdriftij
2248 if (epdriftij.gt.epdrift) epdrift=epdriftij
2252 epdrift=0.5d0*epdrift*d_time*d_time
2253 ! write (iout,*) "epdrift",epdrift
2255 end subroutine predict_edrift
2256 !-----------------------------------------------------------------------------
2257 subroutine verlet_bath
2259 ! Coupling to the thermostat by using the Berendsen algorithm
2262 ! implicit real*8 (a-h,o-z)
2263 ! include 'DIMENSIONS'
2264 ! include 'COMMON.CONTROL'
2265 ! include 'COMMON.VAR'
2266 ! include 'COMMON.MD'
2267 ! include 'COMMON.CHAIN'
2268 ! include 'COMMON.DERIV'
2269 ! include 'COMMON.GEO'
2270 ! include 'COMMON.LOCAL'
2271 ! include 'COMMON.INTERACT'
2272 ! include 'COMMON.IOUNITS'
2273 ! include 'COMMON.NAMES'
2274 real(kind=8) :: T_half,fact
2275 integer :: i,j,inres,mnum
2277 T_half=2.0d0/(dimen3*Rb)*EK
2278 fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
2279 ! write(iout,*) "T_half", T_half
2280 ! write(iout,*) "EK", EK
2281 ! write(iout,*) "fact", fact
2283 d_t(j,0)=fact*d_t(j,0)
2287 d_t(j,i)=fact*d_t(j,i)
2292 ! iti=iabs(itype(i,mnum))
2293 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2294 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2295 .and.(mnum.ne.5)) then
2298 d_t(j,inres)=fact*d_t(j,inres)
2303 end subroutine verlet_bath
2304 !-----------------------------------------------------------------------------
2306 ! Set up the initial conditions of a MD simulation
2309 use control, only:tcpu
2310 !el use io_basic, only:ilen
2313 use minimm, only:minim_dc,minimize,sc_move
2314 use io_config, only:readrst
2315 use io, only:statout
2316 ! implicit real*8 (a-h,o-z)
2317 ! include 'DIMENSIONS'
2320 character(len=16) :: form
2321 integer :: IERROR,ERRCODE
2323 integer :: iranmin,itrial,itmp
2324 ! include 'COMMON.SETUP'
2325 ! include 'COMMON.CONTROL'
2326 ! include 'COMMON.VAR'
2327 ! include 'COMMON.MD'
2329 ! include 'COMMON.LANGEVIN'
2331 ! include 'COMMON.LANGEVIN.lang0'
2333 ! include 'COMMON.CHAIN'
2334 ! include 'COMMON.DERIV'
2335 ! include 'COMMON.GEO'
2336 ! include 'COMMON.LOCAL'
2337 ! include 'COMMON.INTERACT'
2338 ! include 'COMMON.IOUNITS'
2339 ! include 'COMMON.NAMES'
2340 ! include 'COMMON.REMD'
2341 real(kind=8),dimension(0:n_ene) :: energia_long,energia_short
2342 real(kind=8),dimension(3) :: vcm,incr,L
2343 real(kind=8) :: xv,sigv,lowb,highb
2344 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
2345 character(len=256) :: qstr
2348 character(len=50) :: tytul
2349 logical :: file_exist
2350 !el common /gucio/ cm
2351 integer :: i,j,ipos,iq,iw,nft_sc,iretcode,nfun,itime,ierr,mnum
2352 real(kind=8) :: etot,tt0
2356 ! write(iout,*) "d_time", d_time
2357 ! Compute the standard deviations of stochastic forces for Langevin dynamics
2358 ! if the friction coefficients do not depend on surface area
2359 if (lang.gt.0 .and. .not.surfarea) then
2362 stdforcp(i)=stdfp(mnum)*dsqrt(gamp(mnum))
2366 stdforcsc(i)=stdfsc(iabs(itype(i,mnum)),mnum) &
2367 *dsqrt(gamsc(iabs(itype(i,mnum)),mnum))
2370 ! Open the pdb file for snapshotshots
2373 if (ilen(tmpdir).gt.0) &
2374 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2375 liczba(:ilen(liczba))//".pdb")
2377 file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2381 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) &
2382 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2383 liczba(:ilen(liczba))//".x")
2384 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2387 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) &
2388 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2389 liczba(:ilen(liczba))//".cx")
2390 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2396 if (ilen(tmpdir).gt.0) &
2397 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
2398 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
2400 if (ilen(tmpdir).gt.0) &
2401 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
2402 cartname=prefix(:ilen(prefix))//"_MD.cx"
2406 write (qstr,'(256(1h ))')
2409 iq = qinfrag(i,iset)*10
2410 iw = wfrag(i,iset)/100
2412 if(me.eq.king.or..not.out1file) &
2413 write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
2414 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
2419 iq = qinpair(i,iset)*10
2420 iw = wpair(i,iset)/100
2422 if(me.eq.king.or..not.out1file) &
2423 write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
2424 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
2428 ! pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
2430 ! cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
2432 ! cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
2434 ! statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
2438 if (restart1file) then
2440 inquire(file=mremd_rst_name,exist=file_exist)
2441 write (*,*) me," Before broadcast: file_exist",file_exist
2443 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,&
2446 write (*,*) me," After broadcast: file_exist",file_exist
2447 ! inquire(file=mremd_rst_name,exist=file_exist)
2448 if(me.eq.king.or..not.out1file) &
2449 write(iout,*) "Initial state read by master and distributed"
2451 if (ilen(tmpdir).gt.0) &
2452 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_' &
2453 //liczba(:ilen(liczba))//'.rst')
2454 inquire(file=rest2name,exist=file_exist)
2457 if(.not.restart1file) then
2458 if(me.eq.king.or..not.out1file) &
2459 write(iout,*) "Initial state will be read from file ",&
2460 rest2name(:ilen(rest2name))
2463 call rescale_weights(t_bath)
2465 if(me.eq.king.or..not.out1file)then
2466 if (restart1file) then
2467 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),&
2470 write(iout,*) "File ",rest2name(:ilen(rest2name)),&
2473 write(iout,*) "Initial velocities randomly generated"
2480 ! Generate initial velocities
2481 if(me.eq.king.or..not.out1file) &
2482 write(iout,*) "Initial velocities randomly generated"
2487 ! rest2name = prefix(:ilen(prefix))//'.rst'
2488 if(me.eq.king.or..not.out1file)then
2489 write (iout,*) "Initial velocities"
2491 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2492 (d_t(j,i+nres),j=1,3)
2494 ! Zeroing the total angular momentum of the system
2495 write(iout,*) "Calling the zero-angular momentum subroutine"
2498 ! Getting the potential energy and forces and velocities and accelerations
2500 ! write (iout,*) "velocity of the center of the mass:"
2501 ! write (iout,*) (vcm(j),j=1,3)
2503 d_t(j,0)=d_t(j,0)-vcm(j)
2505 ! Removing the velocity of the center of mass
2507 if(me.eq.king.or..not.out1file)then
2508 write (iout,*) "vcm right after adjustment:"
2509 write (iout,*) (vcm(j),j=1,3)
2513 if(iranconf.ne.0 .or.indpdb.gt.0.and..not.unres_pdb .or.preminim) then
2515 print *, 'Calling OVERLAP_SC'
2516 call overlap_sc(fail)
2517 print *,'after OVERLAP'
2520 print *,'call SC_MOVE'
2521 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
2522 print *,'SC_move',nft_sc,etot
2523 if(me.eq.king.or..not.out1file) &
2524 write(iout,*) 'SC_move',nft_sc,etot
2528 print *, 'Calling MINIM_DC'
2529 call minim_dc(etot,iretcode,nfun)
2531 call geom_to_var(nvar,varia)
2532 print *,'Calling MINIMIZE.'
2533 call minimize(etot,varia,iretcode,nfun)
2534 call var_to_geom(nvar,varia)
2536 if(me.eq.king.or..not.out1file) &
2537 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2539 if(iranconf.ne.0) then
2540 !c 8/22/17 AL Loop to produce a low-energy random conformation
2543 if(me.eq.king.or..not.out1file) &
2544 write (iout,*) 'Calling OVERLAP_SC'
2545 call overlap_sc(fail)
2546 endif !endif overlap
2549 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
2550 print *,'SC_move',nft_sc,etot
2551 if(me.eq.king.or..not.out1file) &
2552 write(iout,*) 'SC_move',nft_sc,etot
2556 print *, 'Calling MINIM_DC'
2557 call minim_dc(etot,iretcode,nfun)
2559 call geom_to_var(nvar,varia)
2560 print *,'Calling MINIMIZE.'
2561 call minimize(etot,varia,iretcode,nfun)
2562 call var_to_geom(nvar,varia)
2564 if(me.eq.king.or..not.out1file) &
2565 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2567 if (isnan(etot) .or. etot.gt.4.0d6) then
2568 write (iout,*) "Energy too large",etot, &
2569 " trying another random conformation"
2572 call gen_rand_conf(itmp,*30)
2574 30 write (iout,*) 'Failed to generate random conformation', &
2576 write (*,*) 'Processor:',me, &
2577 ' Failed to generate random conformation',&
2586 write (iout,'(a,i3,a)') 'Processor:',me, &
2587 ' error in generating random conformation.'
2588 write (*,'(a,i3,a)') 'Processor:',me, &
2589 ' error in generating random conformation.'
2592 ! call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
2593 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
2603 write (iout,'(a,i3,a)') 'Processor:',me, &
2604 ' failed to generate a low-energy random conformation.'
2605 write (*,'(a,i3,a,f10.3)') 'Processor:',me, &
2606 ' failed to generate a low-energy random conformation.',etot
2610 ! call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
2611 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
2618 call chainbuild_cart
2623 kinetic_T=2.0d0/(dimen3*Rb)*EK
2624 if(me.eq.king.or..not.out1file)then
2634 write(iout,*) "before ETOTAL"
2635 call etotal(potEcomp)
2636 if (large) call enerprint(potEcomp)
2639 t_etotal=t_etotal+MPI_Wtime()-tt0
2641 t_etotal=t_etotal+tcpu()-tt0
2648 if (amax*d_time .gt. dvmax) then
2649 d_time=d_time*dvmax/amax
2650 if(me.eq.king.or..not.out1file) write (iout,*) &
2651 "Time step reduced to",d_time,&
2652 " because of too large initial acceleration."
2654 if(me.eq.king.or..not.out1file)then
2655 write(iout,*) "Potential energy and its components"
2656 call enerprint(potEcomp)
2657 ! write(iout,*) (potEcomp(i),i=0,n_ene)
2659 potE=potEcomp(0)-potEcomp(20)
2662 if (ntwe.ne.0) call statout(itime)
2663 if(me.eq.king.or..not.out1file) &
2664 write (iout,'(/a/3(a25,1pe14.5/))') "Initial:", &
2665 " Kinetic energy",EK," Potential energy",potE, &
2666 " Total energy",totE," Maximum acceleration ", &
2669 write (iout,*) "Initial coordinates"
2671 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),&
2674 write (iout,*) "Initial dC"
2676 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),&
2677 (dc(j,i+nres),j=1,3)
2679 write (iout,*) "Initial velocities"
2680 write (iout,"(13x,' backbone ',23x,' side chain')")
2682 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2683 (d_t(j,i+nres),j=1,3)
2685 write (iout,*) "Initial accelerations"
2687 ! write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2688 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),&
2689 (d_a(j,i+nres),j=1,3)
2695 d_t_old(j,i)=d_t(j,i)
2696 d_a_old(j,i)=d_a(j,i)
2698 ! write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
2707 call etotal_short(energia_short)
2708 if (large) call enerprint(potEcomp)
2711 t_eshort=t_eshort+MPI_Wtime()-tt0
2713 t_eshort=t_eshort+tcpu()-tt0
2718 if(.not.out1file .and. large) then
2719 write (iout,*) "energia_long",energia_long(0),&
2720 " energia_short",energia_short(0),&
2721 " total",energia_long(0)+energia_short(0)
2722 write (iout,*) "Initial fast-force accelerations"
2724 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
2725 (d_a(j,i+nres),j=1,3)
2728 ! 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
2731 d_a_short(j,i)=d_a(j,i)
2740 call etotal_long(energia_long)
2741 if (large) call enerprint(potEcomp)
2744 t_elong=t_elong+MPI_Wtime()-tt0
2746 t_elong=t_elong+tcpu()-tt0
2751 if(.not.out1file .and. large) then
2752 write (iout,*) "energia_long",energia_long(0)
2753 write (iout,*) "Initial slow-force accelerations"
2755 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
2756 (d_a(j,i+nres),j=1,3)
2760 t_enegrad=t_enegrad+MPI_Wtime()-tt0
2762 t_enegrad=t_enegrad+tcpu()-tt0
2766 end subroutine init_MD
2767 !-----------------------------------------------------------------------------
2768 subroutine random_vel
2770 ! implicit real*8 (a-h,o-z)
2772 use random, only:anorm_distr
2774 ! include 'DIMENSIONS'
2775 ! include 'COMMON.CONTROL'
2776 ! include 'COMMON.VAR'
2777 ! include 'COMMON.MD'
2779 ! include 'COMMON.LANGEVIN'
2781 ! include 'COMMON.LANGEVIN.lang0'
2783 ! include 'COMMON.CHAIN'
2784 ! include 'COMMON.DERIV'
2785 ! include 'COMMON.GEO'
2786 ! include 'COMMON.LOCAL'
2787 ! include 'COMMON.INTERACT'
2788 ! include 'COMMON.IOUNITS'
2789 ! include 'COMMON.NAMES'
2790 ! include 'COMMON.TIME1'
2791 real(kind=8) :: xv,sigv,lowb,highb ,Ek1
2794 real(kind=8) ,allocatable, dimension(:) :: DDU1,DDU2,DL2,DL1,xsolv,DML,rs
2795 real(kind=8) :: sumx
2797 real(kind=8) ,allocatable, dimension(:) :: rsold
2798 real (kind=8),allocatable,dimension(:,:) :: matold
2802 integer :: i,j,ii,k,ind,mark,imark,mnum
2803 ! Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
2804 ! First generate velocities in the eigenspace of the G matrix
2805 ! write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
2812 sigv=dsqrt((Rb*t_bath)/geigen(i))
2815 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
2817 write (iout,*) "i",i," ii",ii," geigen",geigen(i),&
2818 " d_t_work_new",d_t_work_new(ii)
2829 Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
2832 write (iout,*) "Ek from eigenvectors",Ek1
2833 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
2837 ! Transform velocities to UNRES coordinate space
2838 allocate (DL1(2*nres))
2839 allocate (DDU1(2*nres))
2840 allocate (DL2(2*nres))
2841 allocate (DDU2(2*nres))
2842 allocate (xsolv(2*nres))
2843 allocate (DML(2*nres))
2844 allocate (rs(2*nres))
2846 allocate (rsold(2*nres))
2847 allocate (matold(2*nres,2*nres))
2849 matold(1,1)=DMorig(1)
2850 matold(1,2)=DU1orig(1)
2851 matold(1,3)=DU2orig(1)
2852 write (*,*) DMorig(1),DU1orig(1),DU2orig(1)
2857 matold(i,j)=DMorig(i)
2858 matold(i,j-1)=DU1orig(i-1)
2860 matold(i,j-2)=DU2orig(i-2)
2868 matold(i,j+1)=DU1orig(i)
2874 matold(i,j+2)=DU2orig(i)
2878 matold(dimen,dimen)=DMorig(dimen)
2879 matold(dimen,dimen-1)=DU1orig(dimen-1)
2880 matold(dimen,dimen-2)=DU2orig(dimen-2)
2881 write(iout,*) "old gmatrix"
2882 call matout(dimen,dimen,2*nres,2*nres,matold)
2886 ! Find the ith eigenvector of the pentadiagonal inertiq matrix
2890 DML(j)=DMorig(j)-geigen(i)
2893 DML(j-1)=DMorig(j)-geigen(i)
2898 DDU1(imark-1)=DU2orig(imark-1)
2899 do j=imark+1,dimen-1
2900 DDU1(j-1)=DU1orig(j)
2908 DDU2(j)=DU2orig(j+1)
2917 write (iout,*) "DL2,DL1,DML,DDU1,DDU2"
2918 write(iout,'(10f10.5)') (DL2(k),k=3,dimen-1)
2919 write(iout,'(10f10.5)') (DL1(k),k=2,dimen-1)
2920 write(iout,'(10f10.5)') (DML(k),k=1,dimen-1)
2921 write(iout,'(10f10.5)') (DDU1(k),k=1,dimen-2)
2922 write(iout,'(10f10.5)') (DDU2(k),k=1,dimen-3)
2925 if (imark.gt.2) rs(imark-2)=-DU2orig(imark-2)
2926 if (imark.gt.1) rs(imark-1)=-DU1orig(imark-1)
2927 if (imark.lt.dimen) rs(imark)=-DU1orig(imark)
2928 if (imark.lt.dimen-1) rs(imark+1)=-DU2orig(imark)
2932 ! write (iout,*) "Vector rs"
2934 ! write (iout,*) j,rs(j)
2937 call FDIAG(dimen-1,DL2,DL1,DML,DDU1,DDU2,rs,xsolv,mark)
2944 sumx=-geigen(i)*xsolv(j)
2946 sumx=sumx+matold(j,k)*xsolv(k)
2949 sumx=sumx+matold(j,k)*xsolv(k-1)
2951 write(iout,'(i5,3f10.5)') j,sumx,rsold(j),sumx-rsold(j)
2954 sumx=-geigen(i)*xsolv(j-1)
2956 sumx=sumx+matold(j,k)*xsolv(k)
2959 sumx=sumx+matold(j,k)*xsolv(k-1)
2961 write(iout,'(i5,3f10.5)') j-1,sumx,rsold(j-1),sumx-rsold(j-1)
2965 "Solution of equations system",i," for eigenvalue",geigen(i)
2967 write(iout,'(i5,f10.5)') j,xsolv(j)
2970 do j=dimen-1,imark,-1
2975 write (iout,*) "Un-normalized eigenvector",i," for eigenvalue",geigen(i)
2977 write(iout,'(i5,f10.5)') j,xsolv(j)
2980 ! Normalize ith eigenvector
2983 sumx=sumx+xsolv(j)**2
2987 xsolv(j)=xsolv(j)/sumx
2990 write (iout,*) "Eigenvector",i," for eigenvalue",geigen(i)
2992 write(iout,'(i5,3f10.5)') j,xsolv(j)
2995 ! All done at this point for eigenvector i, exit loop
3003 write (iout,*) "Unable to find eigenvector",i
3006 ! write (iout,*) "i=",i
3008 ! write (iout,*) "k=",k
3011 ! write(iout,*) "ind",ind," ind1",3*(i-1)+k
3012 d_t_work(ind)=d_t_work(ind) &
3013 +xsolv(j)*d_t_work_new(3*(i-1)+k)
3016 enddo ! i (loop over the eigenvectors)
3019 write (iout,*) "d_t_work"
3021 write (iout,"(i5,f10.5)") i,d_t_work(i)
3026 ! if (itype(i,1).eq.10) then
3028 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
3029 iti=iabs(itype(i,mnum))
3030 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3031 if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
3032 .or.(mnum.eq.5)) then
3039 ! write (iout,*) "k",k," ii+k",ii+k," ii+j+k",ii+j+k,"EK1",Ek1
3040 Ek1=Ek1+0.5d0*mp(mnum)*((d_t_work(ii+j+k)+d_t_work_new(ii+k))/2)**2+&
3041 0.5d0*0.25d0*IP(mnum)*(d_t_work(ii+j+k)-d_t_work_new(ii+k))**2
3044 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3045 .and.(mnum.ne.5)) ii=ii+3
3046 write (iout,*) "i",i," itype",itype(i,mnum)," mass",msc(itype(i,mnum),mnum)
3047 write (iout,*) "ii",ii
3050 write (iout,*) "k",k," ii",ii,"EK1",EK1
3051 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3053 Ek1=Ek1+0.5d0*Isc(iabs(itype(i,mnum)),mnum)*(d_t_work(ii)-d_t_work(ii-3))**2
3054 Ek1=Ek1+0.5d0*msc(iabs(itype(i,mnum)),mnum)*d_t_work(ii)**2
3056 write (iout,*) "i",i," ii",ii
3058 write (iout,*) "Ek from d_t_work",Ek1
3059 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
3061 if(allocated(DDU1)) deallocate(DDU1)
3062 if(allocated(DDU2)) deallocate(DDU2)
3063 if(allocated(DL2)) deallocate(DL2)
3064 if(allocated(DL1)) deallocate(DL1)
3065 if(allocated(xsolv)) deallocate(xsolv)
3066 if(allocated(DML)) deallocate(DML)
3067 if(allocated(rs)) deallocate(rs)
3069 if(allocated(matold)) deallocate(matold)
3070 if(allocated(rsold)) deallocate(rsold)
3075 d_t(k,j)=d_t_work(ind)
3079 if (itype(j,1).ne.10 .and. itype(j,mnum).ne.ntyp1_molec(mnum)&
3080 .and.(mnum.ne.5)) then
3082 d_t(k,j+nres)=d_t_work(ind)
3088 write (iout,*) "Random velocities in the Calpha,SC space"
3090 write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3)
3093 write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3)
3100 ! if (itype(i,1).eq.10) then
3102 if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
3103 .or.(mnum.eq.5)) then
3105 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
3109 d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
3110 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
3115 write (iout,*)"Random velocities in the virtual-bond-vector space"
3117 write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3)
3120 write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3)
3123 write (iout,*) "Ek from d_t_work",Ek1
3124 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
3132 d_t_work(ind)=d_t_work(ind) &
3133 +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
3135 ! write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
3139 ! Transfer to the d_t vector
3141 d_t(j,0)=d_t_work(j)
3147 d_t(j,i)=d_t_work(ind)
3152 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3153 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3154 .and.(mnum.ne.5)) then
3157 d_t(j,i+nres)=d_t_work(ind)
3163 ! write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",&
3164 ! 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
3166 ! write(iout,*) "end init MD"
3168 end subroutine random_vel
3169 !-----------------------------------------------------------------------------
3171 subroutine sd_verlet_p_setup
3172 ! Sets up the parameters of stochastic Verlet algorithm
3173 ! implicit real*8 (a-h,o-z)
3174 ! include 'DIMENSIONS'
3175 use control, only: tcpu
3180 ! include 'COMMON.CONTROL'
3181 ! include 'COMMON.VAR'
3182 ! include 'COMMON.MD'
3184 ! include 'COMMON.LANGEVIN'
3186 ! include 'COMMON.LANGEVIN.lang0'
3188 ! include 'COMMON.CHAIN'
3189 ! include 'COMMON.DERIV'
3190 ! include 'COMMON.GEO'
3191 ! include 'COMMON.LOCAL'
3192 ! include 'COMMON.INTERACT'
3193 ! include 'COMMON.IOUNITS'
3194 ! include 'COMMON.NAMES'
3195 ! include 'COMMON.TIME1'
3196 real(kind=8),dimension(6*nres) :: emgdt !(MAXRES6) maxres6=6*maxres
3197 real(kind=8) :: pterm,vterm,rho,rhoc,vsig
3198 real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
3199 prand_vec,vrand_vec1,vrand_vec2 !(MAXRES6) maxres6=6*maxres
3200 logical :: lprn = .false.
3201 real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
3202 real(kind=8) :: ktm,gdt,egdt,gdt2,gdt3,gdt4,gdt5,gdt6,gdt7,gdt8,&
3204 integer :: i,maxres2
3211 ! AL 8/17/04 Code adapted from tinker
3213 ! Get the frictional and random terms for stochastic dynamics in the
3214 ! eigenspace of mass-scaled UNRES friction matrix
3218 gdt = fricgam(i) * d_time
3220 ! Stochastic dynamics reduces to simple MD for zero friction
3222 if (gdt .le. zero) then
3223 pfric_vec(i) = 1.0d0
3224 vfric_vec(i) = d_time
3225 afric_vec(i) = 0.5d0 * d_time * d_time
3226 prand_vec(i) = 0.0d0
3227 vrand_vec1(i) = 0.0d0
3228 vrand_vec2(i) = 0.0d0
3230 ! Analytical expressions when friction coefficient is large
3233 if (gdt .ge. gdt_radius) then
3236 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
3237 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
3238 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
3239 vterm = 1.0d0 - egdt**2
3240 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
3242 ! Use series expansions when friction coefficient is small
3253 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0 &
3254 - gdt5/120.0d0 + gdt6/720.0d0 &
3255 - gdt7/5040.0d0 + gdt8/40320.0d0 &
3256 - gdt9/362880.0d0) / fricgam(i)**2
3257 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
3258 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
3259 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0 &
3260 + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0 &
3261 + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0 &
3262 + 127.0d0*gdt9/90720.0d0
3263 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0 &
3264 - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0 &
3265 - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0 &
3266 - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
3267 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0 &
3268 - 17.0d0*gdt2/1280.0d0 &
3269 + 17.0d0*gdt3/6144.0d0 &
3270 + 40967.0d0*gdt4/34406400.0d0 &
3271 - 57203.0d0*gdt5/275251200.0d0 &
3272 - 1429487.0d0*gdt6/13212057600.0d0)
3275 ! Compute the scaling factors of random terms for the nonzero friction case
3277 ktm = 0.5d0*d_time/fricgam(i)
3278 psig = dsqrt(ktm*pterm) / fricgam(i)
3279 vsig = dsqrt(ktm*vterm)
3280 rhoc = dsqrt(1.0d0 - rho*rho)
3282 vrand_vec1(i) = vsig * rho
3283 vrand_vec2(i) = vsig * rhoc
3288 "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",&
3291 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),&
3292 afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
3296 ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
3299 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
3300 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
3301 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
3302 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
3303 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
3304 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
3307 t_sdsetup=t_sdsetup+MPI_Wtime()
3309 t_sdsetup=t_sdsetup+tcpu()-tt0
3312 end subroutine sd_verlet_p_setup
3313 !-----------------------------------------------------------------------------
3314 subroutine eigtransf1(n,ndim,ab,d,c)
3318 real(kind=8) :: ab(ndim,ndim,n),c(ndim,n),d(ndim)
3324 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
3329 end subroutine eigtransf1
3330 !-----------------------------------------------------------------------------
3331 subroutine eigtransf(n,ndim,a,b,d,c)
3335 real(kind=8) :: a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
3341 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
3346 end subroutine eigtransf
3347 !-----------------------------------------------------------------------------
3348 subroutine sd_verlet1
3350 ! Applying stochastic velocity Verlet algorithm - step 1 to velocities
3352 ! implicit real*8 (a-h,o-z)
3353 ! include 'DIMENSIONS'
3354 ! include 'COMMON.CONTROL'
3355 ! include 'COMMON.VAR'
3356 ! include 'COMMON.MD'
3358 ! include 'COMMON.LANGEVIN'
3360 ! include 'COMMON.LANGEVIN.lang0'
3362 ! include 'COMMON.CHAIN'
3363 ! include 'COMMON.DERIV'
3364 ! include 'COMMON.GEO'
3365 ! include 'COMMON.LOCAL'
3366 ! include 'COMMON.INTERACT'
3367 ! include 'COMMON.IOUNITS'
3368 ! include 'COMMON.NAMES'
3369 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
3370 !el common /stochcalc/ stochforcvec
3371 logical :: lprn = .false.
3372 real(kind=8) :: ddt1,ddt2
3373 integer :: i,j,ind,inres
3375 ! write (iout,*) "dc_old"
3377 ! write (iout,'(i5,3f10.5,5x,3f10.5)')
3378 ! & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3381 dc_work(j)=dc_old(j,0)
3382 d_t_work(j)=d_t_old(j,0)
3383 d_a_work(j)=d_a_old(j,0)
3388 dc_work(ind+j)=dc_old(j,i)
3389 d_t_work(ind+j)=d_t_old(j,i)
3390 d_a_work(ind+j)=d_a_old(j,i)
3396 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3397 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3398 .and.(mnum.ne.5)) then
3400 dc_work(ind+j)=dc_old(j,i+nres)
3401 d_t_work(ind+j)=d_t_old(j,i+nres)
3402 d_a_work(ind+j)=d_a_old(j,i+nres)
3410 "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",&
3414 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
3415 vfric_mat(i,j),afric_mat(i,j),&
3416 prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3424 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) &
3425 +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3426 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3427 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3429 d_t_work_new(i)=ddt1+0.5d0*ddt2
3430 d_t_work(i)=ddt1+ddt2
3435 d_t(j,0)=d_t_work(j)
3440 dc(j,i)=dc_work(ind+j)
3441 d_t(j,i)=d_t_work(ind+j)
3447 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3448 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3449 .and.(mnum.ne.5)) then
3452 dc(j,inres)=dc_work(ind+j)
3453 d_t(j,inres)=d_t_work(ind+j)
3459 end subroutine sd_verlet1
3460 !-----------------------------------------------------------------------------
3461 subroutine sd_verlet2
3463 ! Calculating the adjusted velocities for accelerations
3465 ! implicit real*8 (a-h,o-z)
3466 ! include 'DIMENSIONS'
3467 ! include 'COMMON.CONTROL'
3468 ! include 'COMMON.VAR'
3469 ! include 'COMMON.MD'
3471 ! include 'COMMON.LANGEVIN'
3473 ! include 'COMMON.LANGEVIN.lang0'
3475 ! include 'COMMON.CHAIN'
3476 ! include 'COMMON.DERIV'
3477 ! include 'COMMON.GEO'
3478 ! include 'COMMON.LOCAL'
3479 ! include 'COMMON.INTERACT'
3480 ! include 'COMMON.IOUNITS'
3481 ! include 'COMMON.NAMES'
3482 !el real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV !(MAXRES6) maxres6=6*maxres
3483 real(kind=8),dimension(6*nres) :: stochforcvecV !(MAXRES6) maxres6=6*maxres
3484 !el common /stochcalc/ stochforcvec
3486 real(kind=8) :: ddt1,ddt2
3487 integer :: i,j,ind,inres
3488 ! Compute the stochastic forces which contribute to velocity change
3490 call stochastic_force(stochforcvecV)
3497 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3498 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+ &
3499 vrand_mat2(i,j)*stochforcvecV(j)
3501 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
3505 d_t(j,0)=d_t_work(j)
3510 d_t(j,i)=d_t_work(ind+j)
3516 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3517 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3518 .and.(mnum.ne.5)) then
3521 d_t(j,inres)=d_t_work(ind+j)
3527 end subroutine sd_verlet2
3528 !-----------------------------------------------------------------------------
3529 subroutine sd_verlet_ciccotti_setup
3531 ! Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
3533 ! implicit real*8 (a-h,o-z)
3534 ! include 'DIMENSIONS'
3535 use control, only: tcpu
3540 ! include 'COMMON.CONTROL'
3541 ! include 'COMMON.VAR'
3542 ! include 'COMMON.MD'
3544 ! include 'COMMON.LANGEVIN'
3546 ! include 'COMMON.LANGEVIN.lang0'
3548 ! include 'COMMON.CHAIN'
3549 ! include 'COMMON.DERIV'
3550 ! include 'COMMON.GEO'
3551 ! include 'COMMON.LOCAL'
3552 ! include 'COMMON.INTERACT'
3553 ! include 'COMMON.IOUNITS'
3554 ! include 'COMMON.NAMES'
3555 ! include 'COMMON.TIME1'
3556 real(kind=8),dimension(6*nres) :: emgdt !(MAXRES6) maxres6=6*maxres
3557 real(kind=8) :: pterm,vterm,rho,rhoc,vsig
3558 real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
3559 prand_vec,vrand_vec1,vrand_vec2 !(MAXRES6) maxres6=6*maxres
3560 logical :: lprn = .false.
3561 real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
3562 real(kind=8) :: ktm,gdt,egdt,tt0
3563 integer :: i,maxres2
3570 ! AL 8/17/04 Code adapted from tinker
3572 ! Get the frictional and random terms for stochastic dynamics in the
3573 ! eigenspace of mass-scaled UNRES friction matrix
3577 write (iout,*) "i",i," fricgam",fricgam(i)
3578 gdt = fricgam(i) * d_time
3580 ! Stochastic dynamics reduces to simple MD for zero friction
3582 if (gdt .le. zero) then
3583 pfric_vec(i) = 1.0d0
3584 vfric_vec(i) = d_time
3585 afric_vec(i) = 0.5d0*d_time*d_time
3586 prand_vec(i) = afric_vec(i)
3587 vrand_vec2(i) = vfric_vec(i)
3589 ! Analytical expressions when friction coefficient is large
3594 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
3595 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
3596 prand_vec(i) = afric_vec(i)
3597 vrand_vec2(i) = vfric_vec(i)
3599 ! Compute the scaling factors of random terms for the nonzero friction case
3601 ! ktm = 0.5d0*d_time/fricgam(i)
3602 ! psig = dsqrt(ktm*pterm) / fricgam(i)
3603 ! vsig = dsqrt(ktm*vterm)
3604 ! prand_vec(i) = psig*afric_vec(i)
3605 ! vrand_vec2(i) = vsig*vfric_vec(i)
3610 "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",&
3613 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),&
3614 afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
3618 ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
3620 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
3621 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
3622 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
3623 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
3624 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
3626 t_sdsetup=t_sdsetup+MPI_Wtime()
3628 t_sdsetup=t_sdsetup+tcpu()-tt0
3631 end subroutine sd_verlet_ciccotti_setup
3632 !-----------------------------------------------------------------------------
3633 subroutine sd_verlet1_ciccotti
3635 ! Applying stochastic velocity Verlet algorithm - step 1 to velocities
3636 ! implicit real*8 (a-h,o-z)
3638 ! include 'DIMENSIONS'
3642 ! include 'COMMON.CONTROL'
3643 ! include 'COMMON.VAR'
3644 ! include 'COMMON.MD'
3646 ! include 'COMMON.LANGEVIN'
3648 ! include 'COMMON.LANGEVIN.lang0'
3650 ! include 'COMMON.CHAIN'
3651 ! include 'COMMON.DERIV'
3652 ! include 'COMMON.GEO'
3653 ! include 'COMMON.LOCAL'
3654 ! include 'COMMON.INTERACT'
3655 ! include 'COMMON.IOUNITS'
3656 ! include 'COMMON.NAMES'
3657 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
3658 !el common /stochcalc/ stochforcvec
3659 logical :: lprn = .false.
3660 real(kind=8) :: ddt1,ddt2
3661 integer :: i,j,ind,inres
3662 ! write (iout,*) "dc_old"
3664 ! write (iout,'(i5,3f10.5,5x,3f10.5)')
3665 ! & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3668 dc_work(j)=dc_old(j,0)
3669 d_t_work(j)=d_t_old(j,0)
3670 d_a_work(j)=d_a_old(j,0)
3675 dc_work(ind+j)=dc_old(j,i)
3676 d_t_work(ind+j)=d_t_old(j,i)
3677 d_a_work(ind+j)=d_a_old(j,i)
3682 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3684 dc_work(ind+j)=dc_old(j,i+nres)
3685 d_t_work(ind+j)=d_t_old(j,i+nres)
3686 d_a_work(ind+j)=d_a_old(j,i+nres)
3695 "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",&
3699 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
3700 vfric_mat(i,j),afric_mat(i,j),&
3701 prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3709 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) &
3710 +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3711 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3712 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3714 d_t_work_new(i)=ddt1+0.5d0*ddt2
3715 d_t_work(i)=ddt1+ddt2
3720 d_t(j,0)=d_t_work(j)
3725 dc(j,i)=dc_work(ind+j)
3726 d_t(j,i)=d_t_work(ind+j)
3732 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3733 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3734 .and.(mnum.ne.5)) then
3737 dc(j,inres)=dc_work(ind+j)
3738 d_t(j,inres)=d_t_work(ind+j)
3744 end subroutine sd_verlet1_ciccotti
3745 !-----------------------------------------------------------------------------
3746 subroutine sd_verlet2_ciccotti
3748 ! Calculating the adjusted velocities for accelerations
3750 ! implicit real*8 (a-h,o-z)
3751 ! include 'DIMENSIONS'
3752 ! include 'COMMON.CONTROL'
3753 ! include 'COMMON.VAR'
3754 ! include 'COMMON.MD'
3756 ! include 'COMMON.LANGEVIN'
3758 ! include 'COMMON.LANGEVIN.lang0'
3760 ! include 'COMMON.CHAIN'
3761 ! include 'COMMON.DERIV'
3762 ! include 'COMMON.GEO'
3763 ! include 'COMMON.LOCAL'
3764 ! include 'COMMON.INTERACT'
3765 ! include 'COMMON.IOUNITS'
3766 ! include 'COMMON.NAMES'
3767 !el real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV !(MAXRES6) maxres6=6*maxres
3768 real(kind=8),dimension(6*nres) :: stochforcvecV !(MAXRES6) maxres6=6*maxres
3769 !el common /stochcalc/ stochforcvec
3770 real(kind=8) :: ddt1,ddt2
3771 integer :: i,j,ind,inres
3773 ! Compute the stochastic forces which contribute to velocity change
3775 call stochastic_force(stochforcvecV)
3782 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3783 ! ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
3784 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
3786 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
3790 d_t(j,0)=d_t_work(j)
3795 d_t(j,i)=d_t_work(ind+j)
3801 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3803 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3806 d_t(j,inres)=d_t_work(ind+j)
3812 end subroutine sd_verlet2_ciccotti
3814 !-----------------------------------------------------------------------------
3816 !-----------------------------------------------------------------------------
3817 subroutine inertia_tensor
3819 ! Calculating the intertia tensor for the entire protein in order to
3820 ! remove the perpendicular components of velocity matrix which cause
3821 ! the molecule to rotate.
3824 ! implicit real*8 (a-h,o-z)
3825 ! include 'DIMENSIONS'
3826 ! include 'COMMON.CONTROL'
3827 ! include 'COMMON.VAR'
3828 ! include 'COMMON.MD'
3829 ! include 'COMMON.CHAIN'
3830 ! include 'COMMON.DERIV'
3831 ! include 'COMMON.GEO'
3832 ! include 'COMMON.LOCAL'
3833 ! include 'COMMON.INTERACT'
3834 ! include 'COMMON.IOUNITS'
3835 ! include 'COMMON.NAMES'
3837 real(kind=8),dimension(3,3) :: Im,Imcp,eigvec,Id
3838 real(kind=8),dimension(3) :: pr,eigval,L,vp,vrot
3839 real(kind=8) :: M_SC,mag,mag2,M_PEP
3840 real(kind=8),dimension(3,0:nres) :: vpp !(3,0:MAXRES)
3841 real(kind=8),dimension(3) :: vs_p,pp,incr,v
3842 real(kind=8),dimension(3,3) :: pr1,pr2
3844 !el common /gucio/ cm
3845 integer :: iti,inres,i,j,k,mnum
3856 ! calculating the center of the mass of the protein
3860 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
3861 if (itype(i,mnum).eq.ntyp1_molec(mnum)) cycle
3862 M_PEP=M_PEP+mp(mnum)
3864 cm(j)=cm(j)+(c(j,i)+0.5d0*dc(j,i))*mp(mnum)
3873 if (mnum.eq.5) cycle
3874 iti=iabs(itype(i,mnum))
3875 M_SC=M_SC+msc(iabs(iti),mnum)
3878 cm(j)=cm(j)+msc(iabs(iti),mnum)*c(j,inres)
3882 cm(j)=cm(j)/(M_SC+M_PEP)
3887 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
3889 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
3891 Im(1,1)=Im(1,1)+mp(mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
3892 Im(1,2)=Im(1,2)-mp(mnum)*pr(1)*pr(2)
3893 Im(1,3)=Im(1,3)-mp(mnum)*pr(1)*pr(3)
3894 Im(2,3)=Im(2,3)-mp(mnum)*pr(2)*pr(3)
3895 Im(2,2)=Im(2,2)+mp(mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
3896 Im(3,3)=Im(3,3)+mp(mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
3901 iti=iabs(itype(i,mnum))
3902 if (mnum.eq.5) cycle
3905 pr(j)=c(j,inres)-cm(j)
3907 Im(1,1)=Im(1,1)+msc(iabs(iti),mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
3908 Im(1,2)=Im(1,2)-msc(iabs(iti),mnum)*pr(1)*pr(2)
3909 Im(1,3)=Im(1,3)-msc(iabs(iti),mnum)*pr(1)*pr(3)
3910 Im(2,3)=Im(2,3)-msc(iabs(iti),mnum)*pr(2)*pr(3)
3911 Im(2,2)=Im(2,2)+msc(iabs(iti),mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
3912 Im(3,3)=Im(3,3)+msc(iabs(iti),mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
3917 Im(1,1)=Im(1,1)+Ip(mnum)*(1-dc_norm(1,i)*dc_norm(1,i))* &
3918 vbld(i+1)*vbld(i+1)*0.25d0
3919 Im(1,2)=Im(1,2)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(2,i))* &
3920 vbld(i+1)*vbld(i+1)*0.25d0
3921 Im(1,3)=Im(1,3)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(3,i))* &
3922 vbld(i+1)*vbld(i+1)*0.25d0
3923 Im(2,3)=Im(2,3)+Ip(mnum)*(-dc_norm(2,i)*dc_norm(3,i))* &
3924 vbld(i+1)*vbld(i+1)*0.25d0
3925 Im(2,2)=Im(2,2)+Ip(mnum)*(1-dc_norm(2,i)*dc_norm(2,i))* &
3926 vbld(i+1)*vbld(i+1)*0.25d0
3927 Im(3,3)=Im(3,3)+Ip(mnum)*(1-dc_norm(3,i)*dc_norm(3,i))* &
3928 vbld(i+1)*vbld(i+1)*0.25d0
3934 ! if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
3935 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3936 .and.(mnum.ne.5)) then
3937 iti=iabs(itype(i,mnum))
3939 Im(1,1)=Im(1,1)+Isc(iti,mnum)*(1-dc_norm(1,inres)* &
3940 dc_norm(1,inres))*vbld(inres)*vbld(inres)
3941 Im(1,2)=Im(1,2)-Isc(iti,mnum)*(dc_norm(1,inres)* &
3942 dc_norm(2,inres))*vbld(inres)*vbld(inres)
3943 Im(1,3)=Im(1,3)-Isc(iti,mnum)*(dc_norm(1,inres)* &
3944 dc_norm(3,inres))*vbld(inres)*vbld(inres)
3945 Im(2,3)=Im(2,3)-Isc(iti,mnum)*(dc_norm(2,inres)* &
3946 dc_norm(3,inres))*vbld(inres)*vbld(inres)
3947 Im(2,2)=Im(2,2)+Isc(iti,mnum)*(1-dc_norm(2,inres)* &
3948 dc_norm(2,inres))*vbld(inres)*vbld(inres)
3949 Im(3,3)=Im(3,3)+Isc(iti,mnum)*(1-dc_norm(3,inres)* &
3950 dc_norm(3,inres))*vbld(inres)*vbld(inres)
3955 ! write(iout,*) "The angular momentum before adjustment:"
3956 ! write(iout,*) (L(j),j=1,3)
3962 ! Copying the Im matrix for the djacob subroutine
3970 ! Finding the eigenvectors and eignvalues of the inertia tensor
3971 call djacob(3,3,10000,1.0d-10,Imcp,eigvec,eigval)
3972 ! write (iout,*) "Eigenvalues & Eigenvectors"
3973 ! write (iout,'(5x,3f10.5)') (eigval(i),i=1,3)
3976 ! write (iout,'(i5,3f10.5)') i,(eigvec(i,j),j=1,3)
3978 ! Constructing the diagonalized matrix
3980 if (dabs(eigval(i)).gt.1.0d-15) then
3981 Id(i,i)=1.0d0/eigval(i)
3988 Imcp(i,j)=eigvec(j,i)
3994 pr1(i,j)=pr1(i,j)+Id(i,k)*Imcp(k,j)
4001 pr2(i,j)=pr2(i,j)+eigvec(i,k)*pr1(k,j)
4005 ! Calculating the total rotational velocity of the molecule
4008 vrot(i)=vrot(i)+pr2(i,j)*L(j)
4011 ! Resetting the velocities
4013 call vecpr(vrot(1),dc(1,i),vp)
4015 d_t(j,i)=d_t(j,i)-vp(j)
4020 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4021 .and.(mnum.ne.5)) then
4022 ! if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
4023 ! if(itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
4025 call vecpr(vrot(1),dc(1,inres),vp)
4027 d_t(j,inres)=d_t(j,inres)-vp(j)
4032 ! write(iout,*) "The angular momentum after adjustment:"
4033 ! write(iout,*) (L(j),j=1,3)
4036 end subroutine inertia_tensor
4037 !-----------------------------------------------------------------------------
4038 subroutine angmom(cm,L)
4041 ! implicit real*8 (a-h,o-z)
4042 ! include 'DIMENSIONS'
4043 ! include 'COMMON.CONTROL'
4044 ! include 'COMMON.VAR'
4045 ! include 'COMMON.MD'
4046 ! include 'COMMON.CHAIN'
4047 ! include 'COMMON.DERIV'
4048 ! include 'COMMON.GEO'
4049 ! include 'COMMON.LOCAL'
4050 ! include 'COMMON.INTERACT'
4051 ! include 'COMMON.IOUNITS'
4052 ! include 'COMMON.NAMES'
4053 real(kind=8) :: mscab
4054 real(kind=8),dimension(3) :: L,cm,pr,vp,vrot,incr,v,pp
4055 integer :: iti,inres,i,j,mnum
4056 ! Calculate the angular momentum
4065 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
4067 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
4070 v(j)=incr(j)+0.5d0*d_t(j,i)
4073 incr(j)=incr(j)+d_t(j,i)
4075 call vecpr(pr(1),v(1),vp)
4077 L(j)=L(j)+mp(mnum)*vp(j)
4081 pp(j)=0.5d0*d_t(j,i)
4083 call vecpr(pr(1),pp(1),vp)
4085 L(j)=L(j)+Ip(mnum)*vp(j)
4093 iti=iabs(itype(i,mnum))
4101 pr(j)=c(j,inres)-cm(j)
4103 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4104 .and.(mnum.ne.5)) then
4106 v(j)=incr(j)+d_t(j,inres)
4113 call vecpr(pr(1),v(1),vp)
4114 ! write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),&
4115 ! " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
4117 L(j)=L(j)+mscab*vp(j)
4119 ! write (iout,*) "L",(l(j),j=1,3)
4120 if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4121 .and.(mnum.ne.5)) then
4123 v(j)=incr(j)+d_t(j,inres)
4125 call vecpr(dc(1,inres),d_t(1,inres),vp)
4127 L(j)=L(j)+Isc(iti,mnum)*vp(j)
4131 incr(j)=incr(j)+d_t(j,i)
4135 end subroutine angmom
4136 !-----------------------------------------------------------------------------
4137 subroutine vcm_vel(vcm)
4140 ! implicit real*8 (a-h,o-z)
4141 ! include 'DIMENSIONS'
4142 ! include 'COMMON.VAR'
4143 ! include 'COMMON.MD'
4144 ! include 'COMMON.CHAIN'
4145 ! include 'COMMON.DERIV'
4146 ! include 'COMMON.GEO'
4147 ! include 'COMMON.LOCAL'
4148 ! include 'COMMON.INTERACT'
4149 ! include 'COMMON.IOUNITS'
4150 real(kind=8),dimension(3) :: vcm,vv
4151 real(kind=8) :: summas,amas
4161 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
4163 summas=summas+mp(mnum)
4165 vcm(j)=vcm(j)+mp(mnum)*(vv(j)+0.5d0*d_t(j,i))
4169 amas=msc(iabs(itype(i,mnum)),mnum)
4174 if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4175 .and.(mnum.ne.5)) then
4177 vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres))
4181 vcm(j)=vcm(j)+amas*vv(j)
4185 vv(j)=vv(j)+d_t(j,i)
4188 ! write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas
4190 vcm(j)=vcm(j)/summas
4193 end subroutine vcm_vel
4194 !-----------------------------------------------------------------------------
4196 !-----------------------------------------------------------------------------
4198 ! RATTLE algorithm for velocity Verlet - step 1, UNRES
4202 ! implicit real*8 (a-h,o-z)
4203 ! include 'DIMENSIONS'
4205 ! include 'COMMON.CONTROL'
4206 ! include 'COMMON.VAR'
4207 ! include 'COMMON.MD'
4209 ! include 'COMMON.LANGEVIN'
4211 ! include 'COMMON.LANGEVIN.lang0'
4213 ! include 'COMMON.CHAIN'
4214 ! include 'COMMON.DERIV'
4215 ! include 'COMMON.GEO'
4216 ! include 'COMMON.LOCAL'
4217 ! include 'COMMON.INTERACT'
4218 ! include 'COMMON.IOUNITS'
4219 ! include 'COMMON.NAMES'
4220 ! include 'COMMON.TIME1'
4221 !el real(kind=8) :: gginv(2*nres,2*nres),&
4222 !el gdc(3,2*nres,2*nres)
4223 real(kind=8) :: dC_uncor(3,2*nres) !,&
4224 !el real(kind=8) :: Cmat(2*nres,2*nres)
4225 real(kind=8) :: x(2*nres),xcorr(3,2*nres) !maxres2=2*maxres
4226 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
4227 !el common /przechowalnia/ nbond
4228 integer :: max_rattle = 5
4229 logical :: lprn = .false., lprn1 = .false., not_done
4230 real(kind=8) :: tol_rattle = 1.0d-5
4232 integer :: ii,i,j,jj,l,ind,ind1,nres2
4235 !el /common/ przechowalnia
4237 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4238 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4239 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4241 if (lprn) write (iout,*) "RATTLE1"
4245 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4246 .and.(mnum.ne.5)) nbond=nbond+1
4248 ! Make a folded form of the Ginv-matrix
4261 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4266 if (itype(k,1).ne.10 .and. itype(k,mnum).ne.ntyp1_molec(mnum)&
4267 .and.(mnum.ne.5)) then
4271 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4279 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4290 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4294 if (itype(k,1).ne.10) then
4298 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4306 write (iout,*) "Matrix GGinv"
4307 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
4313 if (iter.gt.max_rattle) then
4314 write (iout,*) "Error - too many iterations in RATTLE."
4317 ! Calculate the matrix C = GG**(-1) dC_old o dC
4322 dC_uncor(j,ind1)=dC(j,i)
4326 if (itype(i,1).ne.10) then
4329 dC_uncor(j,ind1)=dC(j,i+nres)
4338 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
4342 if (itype(k,1).ne.10) then
4345 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
4350 ! Calculate deviations from standard virtual-bond lengths
4354 x(ind)=vbld(i+1)**2-vbl**2
4357 if (itype(i,1).ne.10) then
4359 x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4363 write (iout,*) "Coordinates and violations"
4365 write(iout,'(i5,3f10.5,5x,e15.5)') &
4366 i,(dC_uncor(j,i),j=1,3),x(i)
4368 write (iout,*) "Velocities and violations"
4372 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4373 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4377 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4378 .and.(mnum.ne.5)) then
4381 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4382 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4383 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4386 ! write (iout,*) "gdc"
4388 ! write (iout,*) "i",i
4390 ! write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4396 if (dabs(x(i)).gt.xmax) then
4400 if (xmax.lt.tol_rattle) then
4404 ! Calculate the matrix of the system of equations
4409 Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
4414 write (iout,*) "Matrix Cmat"
4415 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4417 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4418 ! Add constraint term to positions
4425 xx = xx+x(ii)*gdc(j,ind,ii)
4429 d_t_new(j,i)=d_t_new(j,i)-xx/d_time
4433 if (itype(i,1).ne.10) then
4438 xx = xx+x(ii)*gdc(j,ind,ii)
4441 dC(j,i+nres)=dC(j,i+nres)-xx
4442 d_t_new(j,i+nres)=d_t_new(j,i+nres)-xx/d_time
4446 ! Rebuild the chain using the new coordinates
4447 call chainbuild_cart
4449 write (iout,*) "New coordinates, Lagrange multipliers,",&
4450 " and differences between actual and standard bond lengths"
4454 xx=vbld(i+1)**2-vbl**2
4455 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4456 i,(dC(j,i),j=1,3),x(ind),xx
4460 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4463 xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4464 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4465 i,(dC(j,i+nres),j=1,3),x(ind),xx
4468 write (iout,*) "Velocities and violations"
4472 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4473 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4476 if (itype(i,1).ne.10) then
4478 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4479 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4480 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4487 10 write (iout,*) "Error - singularity in solving the system",&
4488 " of equations for Lagrange multipliers."
4492 "RATTLE inactive; use -DRATTLE switch at compile time."
4495 end subroutine rattle1
4496 !-----------------------------------------------------------------------------
4498 ! RATTLE algorithm for velocity Verlet - step 2, UNRES
4502 ! implicit real*8 (a-h,o-z)
4503 ! include 'DIMENSIONS'
4505 ! include 'COMMON.CONTROL'
4506 ! include 'COMMON.VAR'
4507 ! include 'COMMON.MD'
4509 ! include 'COMMON.LANGEVIN'
4511 ! include 'COMMON.LANGEVIN.lang0'
4513 ! include 'COMMON.CHAIN'
4514 ! include 'COMMON.DERIV'
4515 ! include 'COMMON.GEO'
4516 ! include 'COMMON.LOCAL'
4517 ! include 'COMMON.INTERACT'
4518 ! include 'COMMON.IOUNITS'
4519 ! include 'COMMON.NAMES'
4520 ! include 'COMMON.TIME1'
4521 !el real(kind=8) :: gginv(2*nres,2*nres),&
4522 !el gdc(3,2*nres,2*nres)
4523 real(kind=8) :: dC_uncor(3,2*nres) !,&
4524 !el Cmat(2*nres,2*nres)
4525 real(kind=8) :: x(2*nres) !maxres2=2*maxres
4526 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
4527 !el common /przechowalnia/ nbond
4528 integer :: max_rattle = 5
4529 logical :: lprn = .false., lprn1 = .false., not_done
4530 real(kind=8) :: tol_rattle = 1.0d-5
4534 !el /common/ przechowalnia
4535 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4536 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4537 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4539 if (lprn) write (iout,*) "RATTLE2"
4540 if (lprn) write (iout,*) "Velocity correction"
4541 ! Calculate the matrix G dC
4547 gdc(j,i,ind)=GGinv(i,ind)*dC(j,k)
4552 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4553 .and.(mnum.ne.5)) then
4556 gdc(j,i,ind)=GGinv(i,ind)*dC(j,k+nres)
4562 ! write (iout,*) "gdc"
4564 ! write (iout,*) "i",i
4566 ! write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4570 ! Calculate the matrix of the system of equations
4577 Cmat(ind,j)=Cmat(ind,j)+dC(k,i)*gdc(k,ind,j)
4583 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4584 .and.(mnum.ne.5)) then
4589 Cmat(ind,j)=Cmat(ind,j)+dC(k,i+nres)*gdc(k,ind,j)
4594 ! Calculate the scalar product dC o d_t_new
4598 x(ind)=scalar(d_t(1,i),dC(1,i))
4602 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4603 .and.(mnum.ne.5)) then
4605 x(ind)=scalar(d_t(1,i+nres),dC(1,i+nres))
4609 write (iout,*) "Velocities and violations"
4613 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4614 i,ind,(d_t(j,i),j=1,3),x(ind)
4618 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4619 .and.(mnum.ne.5)) then
4621 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4622 i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind)
4628 if (dabs(x(i)).gt.xmax) then
4632 if (xmax.lt.tol_rattle) then
4637 write (iout,*) "Matrix Cmat"
4638 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4640 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4641 ! Add constraint term to velocities
4648 xx = xx+x(ii)*gdc(j,ind,ii)
4650 d_t(j,i)=d_t(j,i)-xx
4655 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4656 .and.(mnum.ne.5)) then
4661 xx = xx+x(ii)*gdc(j,ind,ii)
4663 d_t(j,i+nres)=d_t(j,i+nres)-xx
4669 "New velocities, Lagrange multipliers violations"
4673 if (lprn) write (iout,'(2i5,3f10.5,5x,2e15.5)') &
4674 i,ind,(d_t(j,i),j=1,3),x(ind),scalar(d_t(1,i),dC(1,i))
4678 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4681 write (iout,'(2i5,3f10.5,5x,2e15.5)') &
4682 i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind),&
4683 scalar(d_t(1,i+nres),dC(1,i+nres))
4689 10 write (iout,*) "Error - singularity in solving the system",&
4690 " of equations for Lagrange multipliers."
4694 "RATTLE inactive; use -DRATTLE option at compile time."
4697 end subroutine rattle2
4698 !-----------------------------------------------------------------------------
4699 subroutine rattle_brown
4700 ! RATTLE/LINCS algorithm for Brownian dynamics, UNRES
4704 ! implicit real*8 (a-h,o-z)
4705 ! include 'DIMENSIONS'
4707 ! include 'COMMON.CONTROL'
4708 ! include 'COMMON.VAR'
4709 ! include 'COMMON.MD'
4711 ! include 'COMMON.LANGEVIN'
4713 ! include 'COMMON.LANGEVIN.lang0'
4715 ! include 'COMMON.CHAIN'
4716 ! include 'COMMON.DERIV'
4717 ! include 'COMMON.GEO'
4718 ! include 'COMMON.LOCAL'
4719 ! include 'COMMON.INTERACT'
4720 ! include 'COMMON.IOUNITS'
4721 ! include 'COMMON.NAMES'
4722 ! include 'COMMON.TIME1'
4723 !el real(kind=8) :: gginv(2*nres,2*nres),&
4724 !el gdc(3,2*nres,2*nres)
4725 real(kind=8) :: dC_uncor(3,2*nres) !,&
4726 !el real(kind=8) :: Cmat(2*nres,2*nres)
4727 real(kind=8) :: x(2*nres) !maxres2=2*maxres
4728 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
4729 !el common /przechowalnia/ nbond
4730 integer :: max_rattle = 5
4731 logical :: lprn = .false., lprn1 = .false., not_done
4732 real(kind=8) :: tol_rattle = 1.0d-5
4736 !el /common/ przechowalnia
4737 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4738 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4739 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4742 if (lprn) write (iout,*) "RATTLE_BROWN"
4745 if (itype(i,1).ne.10) nbond=nbond+1
4747 ! Make a folded form of the Ginv-matrix
4760 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4764 if (itype(k,1).ne.10) then
4768 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4775 if (itype(i,1).ne.10) then
4785 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4789 if (itype(k,1).ne.10) then
4793 if (j.eq.1 .and. l.eq.1)GGinv(ii,jj)=fricmat(ind,ind1)
4801 write (iout,*) "Matrix GGinv"
4802 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
4808 if (iter.gt.max_rattle) then
4809 write (iout,*) "Error - too many iterations in RATTLE."
4812 ! Calculate the matrix C = GG**(-1) dC_old o dC
4817 dC_uncor(j,ind1)=dC(j,i)
4821 if (itype(i,1).ne.10) then
4824 dC_uncor(j,ind1)=dC(j,i+nres)
4833 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
4837 if (itype(k,1).ne.10) then
4840 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
4845 ! Calculate deviations from standard virtual-bond lengths
4849 x(ind)=vbld(i+1)**2-vbl**2
4852 if (itype(i,1).ne.10) then
4854 x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4858 write (iout,*) "Coordinates and violations"
4860 write(iout,'(i5,3f10.5,5x,e15.5)') &
4861 i,(dC_uncor(j,i),j=1,3),x(i)
4863 write (iout,*) "Velocities and violations"
4867 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4868 i,ind,(d_t(j,i),j=1,3),scalar(d_t(1,i),dC_old(1,i))
4871 if (itype(i,1).ne.10) then
4873 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4874 i+nres,ind,(d_t(j,i+nres),j=1,3),&
4875 scalar(d_t(1,i+nres),dC_old(1,i+nres))
4878 write (iout,*) "gdc"
4880 write (iout,*) "i",i
4882 write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4888 if (dabs(x(i)).gt.xmax) then
4892 if (xmax.lt.tol_rattle) then
4896 ! Calculate the matrix of the system of equations
4901 Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
4906 write (iout,*) "Matrix Cmat"
4907 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4909 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4910 ! Add constraint term to positions
4917 xx = xx+x(ii)*gdc(j,ind,ii)
4920 d_t(j,i)=d_t(j,i)+xx/d_time
4925 if (itype(i,1).ne.10) then
4930 xx = xx+x(ii)*gdc(j,ind,ii)
4933 d_t(j,i+nres)=d_t(j,i+nres)+xx/d_time
4934 dC(j,i+nres)=dC(j,i+nres)+xx
4938 ! Rebuild the chain using the new coordinates
4939 call chainbuild_cart
4941 write (iout,*) "New coordinates, Lagrange multipliers,",&
4942 " and differences between actual and standard bond lengths"
4946 xx=vbld(i+1)**2-vbl**2
4947 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4948 i,(dC(j,i),j=1,3),x(ind),xx
4951 if (itype(i,1).ne.10) then
4953 xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4954 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4955 i,(dC(j,i+nres),j=1,3),x(ind),xx
4958 write (iout,*) "Velocities and violations"
4962 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4963 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4966 if (itype(i,1).ne.10) then
4968 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4969 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4970 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4977 10 write (iout,*) "Error - singularity in solving the system",&
4978 " of equations for Lagrange multipliers."
4982 "RATTLE inactive; use -DRATTLE option at compile time"
4985 end subroutine rattle_brown
4986 !-----------------------------------------------------------------------------
4988 !-----------------------------------------------------------------------------
4989 subroutine friction_force
4994 ! implicit real*8 (a-h,o-z)
4995 ! include 'DIMENSIONS'
4996 ! include 'COMMON.VAR'
4997 ! include 'COMMON.CHAIN'
4998 ! include 'COMMON.DERIV'
4999 ! include 'COMMON.GEO'
5000 ! include 'COMMON.LOCAL'
5001 ! include 'COMMON.INTERACT'
5002 ! include 'COMMON.MD'
5004 ! include 'COMMON.LANGEVIN'
5006 ! include 'COMMON.LANGEVIN.lang0'
5008 ! include 'COMMON.IOUNITS'
5009 !el real(kind=8),dimension(6*nres) :: gamvec !(MAXRES6) maxres6=6*maxres
5010 !el common /syfek/ gamvec
5011 real(kind=8) :: vv(3),vvtot(3,nres),v_work(6*nres) !,&
5012 !el ginvfric(2*nres,2*nres) !maxres2=2*maxres
5013 !el common /przechowalnia/ ginvfric
5015 logical :: lprn = .false., checkmode = .false.
5016 integer :: i,j,ind,k,nres2,nres6,mnum
5020 if(.not.allocated(gamvec)) allocate(gamvec(nres6)) !(MAXRES6)
5021 if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
5029 d_t_work(j)=d_t(j,0)
5034 d_t_work(ind+j)=d_t(j,i)
5040 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5041 .and.(mnum.ne.5)) then
5043 d_t_work(ind+j)=d_t(j,i+nres)
5049 call fricmat_mult(d_t_work,fric_work)
5051 if (.not.checkmode) return
5054 write (iout,*) "d_t_work and fric_work"
5056 write (iout,'(i3,2e15.5)') i,d_t_work(i),fric_work(i)
5060 friction(j,0)=fric_work(j)
5065 friction(j,i)=fric_work(ind+j)
5071 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5072 .and.(mnum.ne.5)) then
5073 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
5075 friction(j,i+nres)=fric_work(ind+j)
5081 write(iout,*) "Friction backbone"
5083 write(iout,'(i5,3e15.5,5x,3e15.5)') &
5084 i,(friction(j,i),j=1,3),(d_t(j,i),j=1,3)
5086 write(iout,*) "Friction side chain"
5088 write(iout,'(i5,3e15.5,5x,3e15.5)') &
5089 i,(friction(j,i+nres),j=1,3),(d_t(j,i+nres),j=1,3)
5098 vvtot(j,i)=vv(j)+0.5d0*d_t(j,i)
5099 vvtot(j,i+nres)=vv(j)+d_t(j,i+nres)
5100 vv(j)=vv(j)+d_t(j,i)
5103 write (iout,*) "vvtot backbone and sidechain"
5105 write (iout,'(i5,3e15.5,5x,3e15.5)') i,(vvtot(j,i),j=1,3),&
5106 (vvtot(j,i+nres),j=1,3)
5111 v_work(ind+j)=vvtot(j,i)
5117 v_work(ind+j)=vvtot(j,i+nres)
5121 write (iout,*) "v_work gamvec and site-based friction forces"
5123 write (iout,'(i5,3e15.5)') i,v_work(i),gamvec(i),&
5127 ! fric_work1(i)=0.0d0
5129 ! fric_work1(i)=fric_work1(i)-A(j,i)*gamvec(j)*v_work(j)
5132 ! write (iout,*) "fric_work and fric_work1"
5134 ! write (iout,'(i5,2e15.5)') i,fric_work(i),fric_work1(i)
5140 ginvfric(i,j)=ginvfric(i,j)+ginv(i,k)*fricmat(k,j)
5144 write (iout,*) "ginvfric"
5146 write (iout,'(i5,100f8.3)') i,(ginvfric(i,j),j=1,dimen)
5148 write (iout,*) "symmetry check"
5151 write (iout,*) i,j,ginvfric(i,j)-ginvfric(j,i)
5156 end subroutine friction_force
5157 !-----------------------------------------------------------------------------
5158 subroutine setup_fricmat
5162 use control_data, only:time_Bcast
5163 use control, only:tcpu
5165 ! implicit real*8 (a-h,o-z)
5169 real(kind=8) :: time00
5171 ! include 'DIMENSIONS'
5172 ! include 'COMMON.VAR'
5173 ! include 'COMMON.CHAIN'
5174 ! include 'COMMON.DERIV'
5175 ! include 'COMMON.GEO'
5176 ! include 'COMMON.LOCAL'
5177 ! include 'COMMON.INTERACT'
5178 ! include 'COMMON.MD'
5179 ! include 'COMMON.SETUP'
5180 ! include 'COMMON.TIME1'
5181 ! integer licznik /0/
5184 ! include 'COMMON.LANGEVIN'
5186 ! include 'COMMON.LANGEVIN.lang0'
5188 ! include 'COMMON.IOUNITS'
5190 integer :: i,j,ind,ind1,m
5191 logical :: lprn = .false.
5192 real(kind=8) :: dtdi !el ,gamvec(2*nres)
5193 !el real(kind=8),dimension(2*nres,2*nres) :: ginvfric,fcopy
5194 ! real(kind=8),allocatable,dimension(:,:) :: fcopy
5195 !el real(kind=8),dimension(2*nres*(2*nres+1)/2) :: Ghalf !(mmaxres2) (mmaxres2=(maxres2*(maxres2+1)/2))
5196 !el common /syfek/ gamvec
5197 real(kind=8) :: work(8*2*nres)
5198 integer :: iwork(2*nres)
5199 !el common /przechowalnia/ ginvfric,Ghalf,fcopy
5200 integer :: ii,iti,k,l,nzero,nres2,nres6,ierr,mnum
5204 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
5205 if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
5206 if (fg_rank.ne.king) goto 10
5211 if(.not.allocated(gamvec)) allocate(gamvec(nres2)) !(MAXRES2)
5212 if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
5213 if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
5214 !el allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
5215 if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !maxres2=2*maxres
5217 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
5218 ! Zeroing out fricmat
5224 ! Load the friction coefficients corresponding to peptide groups
5229 gamvec(ind1)=gamp(mnum)
5231 ! Load the friction coefficients corresponding to side chains
5235 gamsc(ntyp1_molec(j),j)=1.0d0
5242 gamvec(ii)=gamsc(iabs(iti),mnum)
5244 if (surfarea) call sdarea(gamvec)
5246 ! write (iout,*) "Matrix A and vector gamma"
5248 ! write (iout,'(i2,$)') i
5250 ! write (iout,'(f4.1,$)') A(i,j)
5252 ! write (iout,'(f8.3)') gamvec(i)
5256 write (iout,*) "Vector gamvec"
5258 write (iout,'(i5,f10.5)') i, gamvec(i)
5262 ! The friction matrix
5267 dtdi=dtdi+A(j,k)*A(j,i)*gamvec(j)
5274 write (iout,'(//a)') "Matrix fricmat"
5275 call matout2(dimen,dimen,nres2,nres2,fricmat)
5277 if (lang.eq.2 .or. lang.eq.3) then
5278 ! Mass-scale the friction matrix if non-direct integration will be performed
5284 Ginvfric(i,j)=Ginvfric(i,j)+ &
5285 Gsqrm(i,k)*Gsqrm(l,j)*fricmat(k,l)
5290 ! Diagonalize the friction matrix
5295 Ghalf(ind)=Ginvfric(i,j)
5298 call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
5301 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
5302 " mass-scaled friction matrix"
5303 call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
5305 ! Precompute matrices for tinker stochastic integrator
5312 mt1(i,j)=mt1(i,j)+fricvec(k,i)*gsqrm(k,j)
5313 mt2(i,j)=mt2(i,j)+fricvec(k,i)*gsqrp(k,j)
5319 else if (lang.eq.4) then
5320 ! Diagonalize the friction matrix
5325 Ghalf(ind)=fricmat(i,j)
5328 call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
5331 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
5333 call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
5335 ! Determine the number of zero eigenvalues of the friction matrix
5336 nzero=max0(dimen-dimen1,0)
5337 ! do while (fricgam(nzero+1).le.1.0d-5 .and. nzero.lt.dimen)
5340 write (iout,*) "Number of zero eigenvalues:",nzero
5345 fricmat(i,j)=fricmat(i,j) &
5346 +fricvec(i,k)*fricvec(j,k)/fricgam(k)
5351 write (iout,'(//a)') "Generalized inverse of fricmat"
5352 call matout(dimen,dimen,nres6,nres6,fricmat)
5357 if (nfgtasks.gt.1) then
5358 if (fg_rank.eq.0) then
5359 ! The matching BROADCAST for fg processors is called in ERGASTULUM
5365 call MPI_Bcast(10,1,MPI_INTEGER,king,FG_COMM,IERROR)
5367 time_Bcast=time_Bcast+MPI_Wtime()-time00
5369 time_Bcast=time_Bcast+tcpu()-time00
5371 ! print *,"Processor",myrank,
5372 ! & " BROADCAST iorder in SETUP_FRICMAT"
5375 write (iout,*) "setup_fricmat licznik"!,licznik !sp
5381 ! Scatter the friction matrix
5382 call MPI_Scatterv(fricmat(1,1),nginv_counts(0),&
5383 nginv_start(0),MPI_DOUBLE_PRECISION,fcopy(1,1),&
5384 myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
5387 time_scatter=time_scatter+MPI_Wtime()-time00
5388 time_scatter_fmat=time_scatter_fmat+MPI_Wtime()-time00
5390 time_scatter=time_scatter+tcpu()-time00
5391 time_scatter_fmat=time_scatter_fmat+tcpu()-time00
5395 do j=1,2*my_ng_count
5396 fricmat(j,i)=fcopy(i,j)
5399 ! write (iout,*) "My chunk of fricmat"
5400 ! call MATOUT2(my_ng_count,dimen,maxres2,maxres2,fcopy)
5404 end subroutine setup_fricmat
5405 !-----------------------------------------------------------------------------
5406 subroutine stochastic_force(stochforcvec)
5409 use random, only:anorm_distr
5410 ! implicit real*8 (a-h,o-z)
5411 ! include 'DIMENSIONS'
5412 use control, only: tcpu
5417 ! include 'COMMON.VAR'
5418 ! include 'COMMON.CHAIN'
5419 ! include 'COMMON.DERIV'
5420 ! include 'COMMON.GEO'
5421 ! include 'COMMON.LOCAL'
5422 ! include 'COMMON.INTERACT'
5423 ! include 'COMMON.MD'
5424 ! include 'COMMON.TIME1'
5426 ! include 'COMMON.LANGEVIN'
5428 ! include 'COMMON.LANGEVIN.lang0'
5430 ! include 'COMMON.IOUNITS'
5432 real(kind=8) :: x,sig,lowb,highb
5433 real(kind=8) :: ff(3),force(3,0:2*nres),zeta2,lowb2
5434 real(kind=8) :: highb2,sig2,forcvec(6*nres),stochforcvec(6*nres)
5435 real(kind=8) :: time00
5436 logical :: lprn = .false.
5437 integer :: i,j,ind,mnum
5441 stochforc(j,i)=0.0d0
5451 ! Compute the stochastic forces acting on bodies. Store in force.
5457 force(j,i)=anorm_distr(x,sig,lowb,highb)
5465 force(j,i+nres)=anorm_distr(x,sig2,lowb2,highb2)
5469 time_fsample=time_fsample+MPI_Wtime()-time00
5471 time_fsample=time_fsample+tcpu()-time00
5473 ! Compute the stochastic forces acting on virtual-bond vectors.
5479 stochforc(j,i)=ff(j)+0.5d0*force(j,i)
5482 ff(j)=ff(j)+force(j,i)
5484 ! if (itype(i+1,1).ne.ntyp1) then
5486 if (itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
5488 stochforc(j,i)=stochforc(j,i)+force(j,i+nres+1)
5489 ff(j)=ff(j)+force(j,i+nres+1)
5494 stochforc(j,0)=ff(j)+force(j,nnt+nres)
5498 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5499 .and.(mnum.ne.5)) then
5500 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
5502 stochforc(j,i+nres)=force(j,i+nres)
5508 stochforcvec(j)=stochforc(j,0)
5513 stochforcvec(ind+j)=stochforc(j,i)
5519 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5520 .and.(mnum.ne.5)) then
5521 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
5523 stochforcvec(ind+j)=stochforc(j,i+nres)
5529 write (iout,*) "stochforcvec"
5531 write(iout,'(i5,e15.5)') i,stochforcvec(i)
5533 write(iout,*) "Stochastic forces backbone"
5535 write(iout,'(i5,3e15.5)') i,(stochforc(j,i),j=1,3)
5537 write(iout,*) "Stochastic forces side chain"
5539 write(iout,'(i5,3e15.5)') &
5540 i,(stochforc(j,i+nres),j=1,3)
5548 write (iout,*) i,ind
5550 forcvec(ind+j)=force(j,i)
5555 write (iout,*) i,ind
5557 forcvec(j+ind)=force(j,i+nres)
5562 write (iout,*) "forcvec"
5566 write (iout,'(2i3,2f10.5)') i,j,force(j,i),&
5573 write (iout,'(2i3,2f10.5)') i,j,force(j,i+nres),&
5582 end subroutine stochastic_force
5583 !-----------------------------------------------------------------------------
5584 subroutine sdarea(gamvec)
5586 ! Scale the friction coefficients according to solvent accessible surface areas
5587 ! Code adapted from TINKER
5591 ! implicit real*8 (a-h,o-z)
5592 ! include 'DIMENSIONS'
5593 ! include 'COMMON.CONTROL'
5594 ! include 'COMMON.VAR'
5595 ! include 'COMMON.MD'
5597 ! include 'COMMON.LANGEVIN'
5599 ! include 'COMMON.LANGEVIN.lang0'
5601 ! include 'COMMON.CHAIN'
5602 ! include 'COMMON.DERIV'
5603 ! include 'COMMON.GEO'
5604 ! include 'COMMON.LOCAL'
5605 ! include 'COMMON.INTERACT'
5606 ! include 'COMMON.IOUNITS'
5607 ! include 'COMMON.NAMES'
5608 real(kind=8),dimension(2*nres) :: radius,gamvec !(maxres2)
5609 real(kind=8),parameter :: twosix = 1.122462048309372981d0
5610 logical :: lprn = .false.
5611 real(kind=8) :: probe,area,ratio
5612 integer :: i,j,ind,iti,mnum
5614 ! determine new friction coefficients every few SD steps
5616 ! set the atomic radii to estimates of sigma values
5618 ! print *,"Entered sdarea"
5624 ! Load peptide group radii
5627 radius(i)=pstok(mnum)
5629 ! Load side chain radii
5633 radius(i+nres)=restok(iti,mnum)
5636 ! write (iout,*) "i",i," radius",radius(i)
5639 radius(i) = radius(i) / twosix
5640 if (radius(i) .ne. 0.0d0) radius(i) = radius(i) + probe
5643 ! scale atomic friction coefficients by accessible area
5645 if (lprn) write (iout,*) &
5646 "Original gammas, surface areas, scaling factors, new gammas, ",&
5647 "std's of stochastic forces"
5650 if (radius(i).gt.0.0d0) then
5651 call surfatom (i,area,radius)
5652 ratio = dmax1(area/(4.0d0*pi*radius(i)**2),1.0d-1)
5653 if (lprn) write (iout,'(i5,3f10.5,$)') &
5654 i,gamvec(ind+1),area,ratio
5657 gamvec(ind) = ratio * gamvec(ind)
5660 stdforcp(i)=stdfp(mnum)*dsqrt(gamvec(ind))
5661 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcp(i)
5665 if (radius(i+nres).gt.0.0d0) then
5666 call surfatom (i+nres,area,radius)
5667 ratio = dmax1(area/(4.0d0*pi*radius(i+nres)**2),1.0d-1)
5668 if (lprn) write (iout,'(i5,3f10.5,$)') &
5669 i,gamvec(ind+1),area,ratio
5672 gamvec(ind) = ratio * gamvec(ind)
5675 stdforcsc(i)=stdfsc(itype(i,mnum),mnum)*dsqrt(gamvec(ind))
5676 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcsc(i)
5681 end subroutine sdarea
5682 !-----------------------------------------------------------------------------
5684 !-----------------------------------------------------------------------------
5687 ! ###################################################
5688 ! ## COPYRIGHT (C) 1996 by Jay William Ponder ##
5689 ! ## All Rights Reserved ##
5690 ! ###################################################
5692 ! ################################################################
5694 ! ## subroutine surfatom -- exposed surface area of an atom ##
5696 ! ################################################################
5699 ! "surfatom" performs an analytical computation of the surface
5700 ! area of a specified atom; a simplified version of "surface"
5702 ! literature references:
5704 ! T. J. Richmond, "Solvent Accessible Surface Area and
5705 ! Excluded Volume in Proteins", Journal of Molecular Biology,
5708 ! L. Wesson and D. Eisenberg, "Atomic Solvation Parameters
5709 ! Applied to Molecular Dynamics of Proteins in Solution",
5710 ! Protein Science, 1, 227-235 (1992)
5712 ! variables and parameters:
5714 ! ir number of atom for which area is desired
5715 ! area accessible surface area of the atom
5716 ! radius radii of each of the individual atoms
5719 subroutine surfatom(ir,area,radius)
5721 ! implicit real*8 (a-h,o-z)
5722 ! include 'DIMENSIONS'
5724 ! include 'COMMON.GEO'
5725 ! include 'COMMON.IOUNITS'
5727 integer :: nsup,nstart_sup
5728 ! double precision c,dc,dc_old,d_c_work,xloc,xrot,dc_norm
5729 ! common /chain/ c(3,maxres2+2),dc(3,0:maxres2),dc_old(3,0:maxres2),
5730 ! & xloc(3,maxres),xrot(3,maxres),dc_norm(3,0:maxres2),
5731 ! & dc_work(MAXRES6),nres,nres0
5732 integer,parameter :: maxarc=300
5736 integer :: mi,ni,narc
5737 integer :: key(maxarc)
5738 integer :: intag(maxarc)
5739 integer :: intag1(maxarc)
5740 real(kind=8) :: area,arcsum
5741 real(kind=8) :: arclen,exang
5742 real(kind=8) :: delta,delta2
5743 real(kind=8) :: eps,rmove
5744 real(kind=8) :: xr,yr,zr
5745 real(kind=8) :: rr,rrsq
5746 real(kind=8) :: rplus,rminus
5747 real(kind=8) :: axx,axy,axz
5748 real(kind=8) :: ayx,ayy
5749 real(kind=8) :: azx,azy,azz
5750 real(kind=8) :: uxj,uyj,uzj
5751 real(kind=8) :: tx,ty,tz
5752 real(kind=8) :: txb,tyb,td
5753 real(kind=8) :: tr2,tr,txr,tyr
5754 real(kind=8) :: tk1,tk2
5755 real(kind=8) :: thec,the,t,tb
5756 real(kind=8) :: txk,tyk,tzk
5757 real(kind=8) :: t1,ti,tf,tt
5758 real(kind=8) :: txj,tyj,tzj
5759 real(kind=8) :: ccsq,cc,xysq
5760 real(kind=8) :: bsqk,bk,cosine
5761 real(kind=8) :: dsqj,gi,pix2
5762 real(kind=8) :: therk,dk,gk
5763 real(kind=8) :: risqk,rik
5764 real(kind=8) :: radius(2*nres) !(maxatm) (maxatm=maxres2)
5765 real(kind=8) :: ri(maxarc),risq(maxarc)
5766 real(kind=8) :: ux(maxarc),uy(maxarc),uz(maxarc)
5767 real(kind=8) :: xc(maxarc),yc(maxarc),zc(maxarc)
5768 real(kind=8) :: xc1(maxarc),yc1(maxarc),zc1(maxarc)
5769 real(kind=8) :: dsq(maxarc),bsq(maxarc)
5770 real(kind=8) :: dsq1(maxarc),bsq1(maxarc)
5771 real(kind=8) :: arci(maxarc),arcf(maxarc)
5772 real(kind=8) :: ex(maxarc),lt(maxarc),gr(maxarc)
5773 real(kind=8) :: b(maxarc),b1(maxarc),bg(maxarc)
5774 real(kind=8) :: kent(maxarc),kout(maxarc)
5775 real(kind=8) :: ther(maxarc)
5776 logical :: moved,top
5777 logical :: omit(maxarc)
5780 ! maxatm = 2*nres !maxres2 maxres2=2*maxres
5781 ! maxlight = 8*maxatm
5784 ! maxtors = 4*maxatm
5786 ! zero out the surface area for the sphere of interest
5789 ! write (2,*) "ir",ir," radius",radius(ir)
5790 if (radius(ir) .eq. 0.0d0) return
5792 ! set the overlap significance and connectivity shift
5796 delta2 = delta * delta
5801 ! store coordinates and radius of the sphere of interest
5809 ! initialize values of some counters and summations
5818 ! test each sphere to see if it overlaps the sphere of interest
5821 if (i.eq.ir .or. radius(i).eq.0.0d0) goto 30
5822 rplus = rr + radius(i)
5824 if (abs(tx) .ge. rplus) goto 30
5826 if (abs(ty) .ge. rplus) goto 30
5828 if (abs(tz) .ge. rplus) goto 30
5830 ! check for sphere overlap by testing distance against radii
5832 xysq = tx*tx + ty*ty
5833 if (xysq .lt. delta2) then
5840 if (rplus-cc .le. delta) goto 30
5841 rminus = rr - radius(i)
5843 ! check to see if sphere of interest is completely buried
5845 if (cc-abs(rminus) .le. delta) then
5846 if (rminus .le. 0.0d0) goto 170
5850 ! check for too many overlaps with sphere of interest
5852 if (io .ge. maxarc) then
5854 20 format (/,' SURFATOM -- Increase the Value of MAXARC')
5858 ! get overlap between current sphere and sphere of interest
5867 gr(io) = (ccsq+rplus*rminus) / (2.0d0*rr*b1(io))
5873 ! case where no other spheres overlap the sphere of interest
5876 area = 4.0d0 * pi * rrsq
5880 ! case where only one sphere overlaps the sphere of interest
5883 area = pix2 * (1.0d0 + gr(1))
5884 area = mod(area,4.0d0*pi) * rrsq
5888 ! case where many spheres intersect the sphere of interest;
5889 ! sort the intersecting spheres by their degree of overlap
5891 call sort2 (io,gr,key)
5894 intag(i) = intag1(k)
5903 ! get radius of each overlap circle on surface of the sphere
5908 risq(i) = rrsq - gi*gi
5909 ri(i) = sqrt(risq(i))
5910 ther(i) = 0.5d0*pi - asin(min(1.0d0,max(-1.0d0,gr(i))))
5913 ! find boundary of inaccessible area on sphere of interest
5916 if (.not. omit(k)) then
5923 ! check to see if J circle is intersecting K circle;
5924 ! get distance between circle centers and sum of radii
5927 if (omit(j)) goto 60
5928 cc = (txk*xc(j)+tyk*yc(j)+tzk*zc(j))/(bk*b(j))
5929 cc = acos(min(1.0d0,max(-1.0d0,cc)))
5930 td = therk + ther(j)
5932 ! check to see if circles enclose separate regions
5934 if (cc .ge. td) goto 60
5936 ! check for circle J completely inside circle K
5938 if (cc+ther(j) .lt. therk) goto 40
5940 ! check for circles that are essentially parallel
5942 if (cc .gt. delta) goto 50
5947 ! check to see if sphere of interest is completely buried
5950 if (pix2-cc .le. td) goto 170
5956 ! find T value of circle intersections
5959 if (omit(k)) goto 110
5974 ! rotation matrix elements
5986 if (.not. omit(j)) then
5991 ! rotate spheres so K vector colinear with z-axis
5993 uxj = txj*axx + tyj*axy - tzj*axz
5994 uyj = tyj*ayy - txj*ayx
5995 uzj = txj*azx + tyj*azy + tzj*azz
5996 cosine = min(1.0d0,max(-1.0d0,uzj/b(j)))
5997 if (acos(cosine) .lt. therk+ther(j)) then
5998 dsqj = uxj*uxj + uyj*uyj
6003 tr2 = risqk*dsqj - tb*tb
6009 ! get T values of intersection for K circle
6012 tb = min(1.0d0,max(-1.0d0,tb))
6014 if (tyb-txr .lt. 0.0d0) tk1 = pix2 - tk1
6016 tb = min(1.0d0,max(-1.0d0,tb))
6018 if (tyb+txr .lt. 0.0d0) tk2 = pix2 - tk2
6019 thec = (rrsq*uzj-gk*bg(j)) / (rik*ri(j)*b(j))
6020 if (abs(thec) .lt. 1.0d0) then
6022 else if (thec .ge. 1.0d0) then
6024 else if (thec .le. -1.0d0) then
6028 ! see if "tk1" is entry or exit point; check t=0 point;
6029 ! "ti" is exit point, "tf" is entry point
6031 cosine = min(1.0d0,max(-1.0d0, &
6032 (uzj*gk-uxj*rik)/(b(j)*rr)))
6033 if ((acos(cosine)-ther(j))*(tk2-tk1) .le. 0.0d0) then
6041 if (narc .ge. maxarc) then
6043 70 format (/,' SURFATOM -- Increase the Value',&
6047 if (tf .le. ti) then
6068 ! special case; K circle without intersections
6070 if (narc .le. 0) goto 90
6072 ! general case; sum up arclength and set connectivity code
6074 call sort2 (narc,arci,key)
6079 if (narc .gt. 1) then
6082 if (t .lt. arci(j)) then
6083 arcsum = arcsum + arci(j) - t
6084 exang = exang + ex(ni)
6086 if (jb .ge. maxarc) then
6088 80 format (/,' SURFATOM -- Increase the Value',&
6093 kent(jb) = maxarc*i + k
6095 kout(jb) = maxarc*k + i
6104 arcsum = arcsum + pix2 - t
6106 exang = exang + ex(ni)
6109 kent(jb) = maxarc*i + k
6111 kout(jb) = maxarc*k + i
6118 arclen = arclen + gr(k)*arcsum
6121 if (arclen .eq. 0.0d0) goto 170
6122 if (jb .eq. 0) goto 150
6124 ! find number of independent boundaries and check connectivity
6128 if (kout(k) .ne. 0) then
6135 if (m .eq. kent(ii)) then
6138 if (j .eq. jb) goto 150
6150 ! attempt to fix connectivity error by moving atom slightly
6154 140 format (/,' SURFATOM -- Connectivity Error at Atom',i6)
6163 ! compute the exposed surface area for the sphere of interest
6166 area = ib*pix2 + exang + arclen
6167 area = mod(area,4.0d0*pi) * rrsq
6169 ! attempt to fix negative area by moving atom slightly
6171 if (area .lt. 0.0d0) then
6174 160 format (/,' SURFATOM -- Negative Area at Atom',i6)
6185 end subroutine surfatom
6186 !----------------------------------------------------------------
6187 !----------------------------------------------------------------
6188 subroutine alloc_MD_arrays
6189 !EL Allocation of arrays used by MD module
6191 integer :: nres2,nres6
6194 !----------------------
6198 allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
6199 allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
6200 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
6201 allocate(fricvec(nres2,nres2))
6202 allocate(pfric_mat(nres2,nres2),vfric_mat(nres2,nres2))
6203 allocate(afric_mat(nres2,nres2),prand_mat(nres2,nres2))
6204 allocate(vrand_mat1(nres2,nres2),vrand_mat2(nres2,nres2)) !(MAXRES2,MAXRES2)
6205 allocate(pfric0_mat(nres2,nres2,0:maxflag_stoch))
6206 allocate(afric0_mat(nres2,nres2,0:maxflag_stoch))
6207 allocate(vfric0_mat(nres2,nres2,0:maxflag_stoch))
6208 allocate(prand0_mat(nres2,nres2,0:maxflag_stoch))
6209 allocate(vrand0_mat1(nres2,nres2,0:maxflag_stoch))
6210 allocate(vrand0_mat2(nres2,nres2,0:maxflag_stoch)) !(MAXRES2,MAXRES2,0:maxflag_stoch)
6211 allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
6213 allocate(mt1(nres2,nres2),mt2(nres2,nres2),mt3(nres2,nres2)) !(maxres2,maxres2)
6214 !----------------------
6216 ! commom.langevin.lang0
6218 allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
6219 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
6220 allocate(fricvec(nres2,nres2)) !(MAXRES2,MAXRES2)
6221 allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
6222 allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
6225 if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2))
6226 !----------------------
6227 ! commom.hairpin in CSA module
6228 !----------------------
6229 ! common.mce in MCM_MD module
6230 !----------------------
6232 ! common /mdgrad/ in module.energy
6233 ! common /back_constr/ in module.energy
6234 ! common /qmeas/ in module.energy
6237 allocate(potEcomp(0:n_ene+4)) !(0:n_ene+4)
6239 allocate(d_t(3,0:nres2),d_a(3,0:nres2),d_t_old(3,0:nres2)) !(3,0:MAXRES2)
6240 allocate(d_a_work(nres6)) !(6*MAXRES)
6242 allocate(DM(nres2),DU1(nres2),DU2(nres2))
6243 allocate(DMorig(nres2),DU1orig(nres2),DU2orig(nres2))
6245 write (iout,*) "Before A Allocation",nfgtasks-1
6247 allocate(Gmat(nres2,nres2),A(nres2,nres2))
6248 if(.not.allocated(Ginv)) allocate(Ginv(nres2,nres2)) !in control: ergastulum
6249 allocate(Gsqrp(nres2,nres2),Gsqrm(nres2,nres2),Gvec(nres2,nres2)) !(maxres2,maxres2)
6251 allocate(Geigen(nres2)) !(maxres2)
6252 if(.not.allocated(vtot)) allocate(vtot(nres2)) !(maxres2)
6253 ! common /inertia/ in io_conf: parmread
6254 ! real(kind=8),dimension(:),allocatable :: ISC,msc !(ntyp+1)
6255 ! common /langevin/in io read_MDpar
6256 ! real(kind=8),dimension(:),allocatable :: gamsc !(ntyp1)
6257 ! real(kind=8),dimension(:),allocatable :: stdfsc !(ntyp)
6258 ! in io_conf: parmread
6259 ! real(kind=8),dimension(:),allocatable :: restok !(ntyp+1)
6260 ! common /mdpmpi/ in control: ergastulum
6261 if(.not.allocated(ng_start)) allocate(ng_start(0:nfgtasks-1))
6262 if(.not.allocated(ng_counts)) allocate(ng_counts(0:nfgtasks-1))
6263 if(.not.allocated(nginv_counts)) allocate(nginv_counts(0:nfgtasks-1)) !(0:MaxProcs-1)
6264 if(.not.allocated(nginv_start)) allocate(nginv_start(0:nfgtasks)) !(0:MaxProcs)
6265 !----------------------
6266 ! common.muca in read_muca
6267 ! common /double_muca/
6268 ! real(kind=8) :: elow,ehigh,factor,hbin,factor_min
6269 ! real(kind=8),dimension(:),allocatable :: emuca,nemuca,&
6270 ! nemuca2,hist !(4*maxres)
6271 ! common /integer_muca/
6272 ! integer :: nmuca,imtime,muca_smooth
6274 ! real(kind=8),dimension(:),allocatable :: elowi,ehighi !(maxprocs)
6275 !----------------------
6277 ! common /mdgrad/ in module.energy
6278 ! common /back_constr/ in module.energy
6279 ! common /qmeas/ in module.energy
6283 allocate(d_t_work(nres6),d_t_work_new(nres6),d_af_work(nres6))
6284 allocate(d_as_work(nres6),kinetic_force(nres6)) !(MAXRES6)
6285 allocate(d_t_new(3,0:nres2),d_a_old(3,0:nres2),d_a_short(3,0:nres2)) !,d_a !(3,0:MAXRES2)
6286 allocate(stdforcp(nres),stdforcsc(nres)) !(MAXRES)
6287 !----------------------
6289 allocate(D_ban(nres6)) !(MAXRES6) maxres6=6*maxres
6290 ! common /stochcalc/ stochforcvec
6291 allocate(stochforcvec(nres6)) !(MAXRES6) maxres6=6*maxres
6292 !----------------------
6294 end subroutine alloc_MD_arrays
6295 !-----------------------------------------------------------------------------
6296 !-----------------------------------------------------------------------------