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 ! write(iout,*) "icount_scale",icount_scale,scalel
1147 if (icount_scale.gt.maxcount_scale) then
1149 "ERROR: too many attempts at scaling down the time step. ",&
1150 "amax=",amax,"epdrift=",epdrift,&
1151 "damax=",damax,"edriftmax=",edriftmax,&
1155 call MPI_Abort(MPI_COMM_WORLD,IERROR,IERRCODE)
1159 ! First step of the velocity Verlet algorithm
1164 else if (lang.eq.3) then
1166 call sd_verlet1_ciccotti
1168 else if (lang.eq.1) then
1173 ! Build the chain from the newly calculated coordinates
1174 call chainbuild_cart
1175 if (rattle) call rattle1
1177 if (large.and. mod(itime,ntwe).eq.0) then
1178 write (iout,*) "Cartesian and internal coordinates: step 1"
1183 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),&
1184 (dc(j,i+nres),j=1,3)
1186 write (iout,*) "Accelerations"
1188 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1189 (d_a(j,i+nres),j=1,3)
1191 write (iout,*) "Velocities, step 1"
1193 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1194 (d_t(j,i+nres),j=1,3)
1203 ! Calculate energy and forces
1205 call etotal(potEcomp)
1206 ! AL 4/17/17: Reduce the steps if NaNs occurred.
1207 if (potEcomp(0).gt.0.99e18 .or. isnan(potEcomp(0)).gt.0) then
1209 ! write (iout,*) "Tu jest problem",potEcomp(0),d_time
1213 if (large.and. mod(itime,ntwe).eq.0) &
1214 call enerprint(potEcomp)
1217 t_etotal=t_etotal+MPI_Wtime()-tt0
1219 t_etotal=t_etotal+tcpu()-tt0
1222 potE=potEcomp(0)-potEcomp(20)
1224 ! Get the new accelerations
1227 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1229 t_enegrad=t_enegrad+tcpu()-tt0
1231 ! Determine maximum acceleration and scale down the timestep if needed
1233 amax=amax/(itime_scal+1)**2
1234 call predict_edrift(epdrift)
1235 ! write(iout,*) "amax=",amax,damax,epdrift,edriftmax,amax/(itime_scal+1)
1237 ! write (iout,*) "before enter if",scalel,icount_scale
1238 if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then
1239 ! write(iout,*) "I enter if"
1240 ! Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
1242 ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax)) &
1244 itime_scal=itime_scal+ifac_time
1245 ! fac_time=dmin1(damax/amax,0.5d0)
1246 fac_time=0.5d0**ifac_time
1247 d_time=d_time*fac_time
1248 if (lang.eq.2 .or. lang.eq.3) then
1250 ! write (iout,*) "Calling sd_verlet_setup: 1"
1251 ! Rescale the stochastic forces and recalculate or restore
1252 ! the matrices of tinker integrator
1253 if (itime_scal.gt.maxflag_stoch) then
1254 if (large) write (iout,'(a,i5,a)') &
1255 "Calculate matrices for stochastic step;",&
1256 " itime_scal ",itime_scal
1258 call sd_verlet_p_setup
1260 call sd_verlet_ciccotti_setup
1262 write (iout,'(2a,i3,a,i3,1h.)') &
1263 "Warning: cannot store matrices for stochastic",&
1264 " integration because the index",itime_scal,&
1265 " is greater than",maxflag_stoch
1266 write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",&
1267 " integration Langevin algorithm for better efficiency."
1268 else if (flag_stoch(itime_scal)) then
1269 if (large) write (iout,'(a,i5,a,l1)') &
1270 "Restore matrices for stochastic step; itime_scal ",&
1271 itime_scal," flag ",flag_stoch(itime_scal)
1274 pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
1275 afric_mat(i,j)=afric0_mat(i,j,itime_scal)
1276 vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
1277 prand_mat(i,j)=prand0_mat(i,j,itime_scal)
1278 vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
1279 vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
1283 if (large) write (iout,'(2a,i5,a,l1)') &
1284 "Calculate & store matrices for stochastic step;",&
1285 " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
1287 call sd_verlet_p_setup
1289 call sd_verlet_ciccotti_setup
1291 flag_stoch(ifac_time)=.true.
1294 pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
1295 afric0_mat(i,j,itime_scal)=afric_mat(i,j)
1296 vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
1297 prand0_mat(i,j,itime_scal)=prand_mat(i,j)
1298 vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
1299 vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
1303 fac_time=1.0d0/dsqrt(fac_time)
1305 stochforcvec(i)=fac_time*stochforcvec(i)
1308 else if (lang.eq.1) then
1309 ! Rescale the accelerations due to stochastic forces
1310 fac_time=1.0d0/dsqrt(fac_time)
1312 d_as_work(i)=d_as_work(i)*fac_time
1315 if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)') &
1316 "itime",itime," Timestep scaled down to ",&
1317 d_time," ifac_time",ifac_time," itime_scal",itime_scal
1319 ! Second step of the velocity Verlet algorithm
1324 else if (lang.eq.3) then
1326 call sd_verlet2_ciccotti
1328 else if (lang.eq.1) then
1333 if (rattle) call rattle2
1336 if (d_time.ne.d_time0) then
1339 if (lang.eq.2 .or. lang.eq.3) then
1340 if (large) write (iout,'(a)') &
1341 "Restore original matrices for stochastic step"
1342 ! write (iout,*) "Calling sd_verlet_setup: 2"
1343 ! Restore the matrices of tinker integrator if the time step has been restored
1346 pfric_mat(i,j)=pfric0_mat(i,j,0)
1347 afric_mat(i,j)=afric0_mat(i,j,0)
1348 vfric_mat(i,j)=vfric0_mat(i,j,0)
1349 prand_mat(i,j)=prand0_mat(i,j,0)
1350 vrand_mat1(i,j)=vrand0_mat1(i,j,0)
1351 vrand_mat2(i,j)=vrand0_mat2(i,j,0)
1359 ! Calculate the kinetic and the total energy and the kinetic temperature
1363 ! call kinetic1(EK1)
1364 ! write (iout,*) "step",itime," EK",EK," EK1",EK1
1366 ! Couple the system to Berendsen bath if needed
1367 if (tbf .and. lang.eq.0) then
1370 kinetic_T=2.0d0/(dimen3*Rb)*EK
1371 ! Backup the coordinates, velocities, and accelerations
1375 d_t_old(j,i)=d_t(j,i)
1376 d_a_old(j,i)=d_a(j,i)
1380 if (mod(itime,ntwe).eq.0 .and. large) then
1381 write (iout,*) "Velocities, step 2"
1383 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1384 (d_t(j,i+nres),j=1,3)
1389 end subroutine velverlet_step
1390 !-----------------------------------------------------------------------------
1391 subroutine RESPA_step(itime)
1392 !-------------------------------------------------------------------------------
1393 ! Perform a single RESPA step.
1394 !-------------------------------------------------------------------------------
1395 ! implicit real*8 (a-h,o-z)
1396 ! include 'DIMENSIONS'
1400 use control, only:tcpu
1402 ! use io_conf, only:cartprint
1405 integer :: IERROR,ERRCODE
1407 ! include 'COMMON.SETUP'
1408 ! include 'COMMON.CONTROL'
1409 ! include 'COMMON.VAR'
1410 ! include 'COMMON.MD'
1412 ! include 'COMMON.LANGEVIN'
1414 ! include 'COMMON.LANGEVIN.lang0'
1416 ! include 'COMMON.CHAIN'
1417 ! include 'COMMON.DERIV'
1418 ! include 'COMMON.GEO'
1419 ! include 'COMMON.LOCAL'
1420 ! include 'COMMON.INTERACT'
1421 ! include 'COMMON.IOUNITS'
1422 ! include 'COMMON.NAMES'
1423 ! include 'COMMON.TIME1'
1424 real(kind=8),dimension(0:n_ene) :: energia_short,energia_long
1425 real(kind=8),dimension(3) :: L,vcm,incr
1426 real(kind=8),dimension(3,0:2*nres) :: dc_old0,d_t_old0,d_a_old0 !(3,0:maxres2) maxres2=2*maxres
1427 logical :: PRINT_AMTS_MSG = .false.
1428 integer :: count,rstcount !ilen,
1430 character(len=50) :: tytul
1431 integer :: maxcount_scale = 10
1432 !el common /gucio/ cm
1433 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1434 !el common /stochcalc/ stochforcvec
1435 integer :: itime,itt,i,j,itsplit
1437 !el common /cipiszcze/ itt
1439 real(kind=8) :: epdrift,tt0,epdriftmax
1442 if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres)) !(MAXRES6) maxres6=6*maxres
1446 if (large.and. mod(itime,ntwe).eq.0) then
1447 write (iout,*) "***************** RESPA itime",itime
1448 write (iout,*) "Cartesian and internal coordinates: step 0"
1450 call pdbout(0.0d0,"cipiszcze",iout)
1452 write (iout,*) "Accelerations from long-range forces"
1454 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1455 (d_a(j,i+nres),j=1,3)
1457 write (iout,*) "Velocities, step 0"
1459 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1460 (d_t(j,i+nres),j=1,3)
1465 ! Perform the initial RESPA step (increment velocities)
1466 ! write (iout,*) "*********************** RESPA ini"
1469 if (mod(itime,ntwe).eq.0 .and. large) then
1470 write (iout,*) "Velocities, end"
1472 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1473 (d_t(j,i+nres),j=1,3)
1477 ! Compute the short-range forces
1483 ! 7/2/2009 commented out
1485 ! call etotal_short(energia_short)
1488 ! 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
1491 d_a(j,i)=d_a_short(j,i)
1495 if (large.and. mod(itime,ntwe).eq.0) then
1496 write (iout,*) "energia_short",energia_short(0)
1497 write (iout,*) "Accelerations from short-range forces"
1499 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1500 (d_a(j,i+nres),j=1,3)
1505 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1507 t_enegrad=t_enegrad+tcpu()-tt0
1512 d_t_old(j,i)=d_t(j,i)
1513 d_a_old(j,i)=d_a(j,i)
1516 ! 6/30/08 A-MTS: attempt at increasing the split number
1519 dc_old0(j,i)=dc_old(j,i)
1520 d_t_old0(j,i)=d_t_old(j,i)
1521 d_a_old0(j,i)=d_a_old(j,i)
1524 if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
1525 if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
1532 ! write (iout,*) "itime",itime," ntime_split",ntime_split
1533 ! Split the time step
1534 d_time=d_time0/ntime_split
1535 ! Perform the short-range RESPA steps (velocity Verlet increments of
1536 ! positions and velocities using short-range forces)
1537 ! write (iout,*) "*********************** RESPA split"
1538 do itsplit=1,ntime_split
1541 else if (lang.eq.2 .or. lang.eq.3) then
1543 call stochastic_force(stochforcvec)
1546 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
1548 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1553 ! First step of the velocity Verlet algorithm
1558 else if (lang.eq.3) then
1560 call sd_verlet1_ciccotti
1562 else if (lang.eq.1) then
1567 ! Build the chain from the newly calculated coordinates
1568 call chainbuild_cart
1569 if (rattle) call rattle1
1571 if (large.and. mod(itime,ntwe).eq.0) then
1572 write (iout,*) "***** ITSPLIT",itsplit
1573 write (iout,*) "Cartesian and internal coordinates: step 1"
1574 call pdbout(0.0d0,"cipiszcze",iout)
1577 write (iout,*) "Velocities, step 1"
1579 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1580 (d_t(j,i+nres),j=1,3)
1589 ! Calculate energy and forces
1591 call etotal_short(energia_short)
1592 ! AL 4/17/17: Exit itime_split loop when energy goes infinite
1593 if (energia_short(0).gt.0.99e20 .or. isnan(energia_short(0)) ) then
1594 if (PRINT_AMTS_MSG) &
1595 write (iout,*) "Infinities/NaNs in energia_short",energia_short(0),"; increasing ntime_split to",ntime_split
1596 ntime_split=ntime_split*2
1597 if (ntime_split.gt.maxtime_split) then
1600 "Cannot rescue the run; aborting job. Retry with a smaller time step"
1602 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1605 "Cannot rescue the run; terminating. Retry with a smaller time step"
1611 if (large.and. mod(itime,ntwe).eq.0) &
1612 call enerprint(energia_short)
1615 t_eshort=t_eshort+MPI_Wtime()-tt0
1617 t_eshort=t_eshort+tcpu()-tt0
1621 ! Get the new accelerations
1623 ! 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
1626 d_a_short(j,i)=d_a(j,i)
1630 if (large.and. mod(itime,ntwe).eq.0) then
1631 write (iout,*)"energia_short",energia_short(0)
1632 write (iout,*) "Accelerations from short-range forces"
1634 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1635 (d_a(j,i+nres),j=1,3)
1640 ! Determine maximum acceleration and scale down the timestep if needed
1642 amax=amax/ntime_split**2
1643 call predict_edrift(epdrift)
1644 if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0) &
1645 write (iout,*) "amax",amax," damax",damax,&
1646 " epdrift",epdrift," epdriftmax",epdriftmax
1647 ! Exit loop and try with increased split number if the change of
1648 ! acceleration is too big
1649 if (amax.gt.damax .or. epdrift.gt.edriftmax) then
1650 if (ntime_split.lt.maxtime_split) then
1652 ntime_split=ntime_split*2
1653 ! AL 4/17/17: We should exit the itime_split loop when acceleration change is too big
1657 dc_old(j,i)=dc_old0(j,i)
1658 d_t_old(j,i)=d_t_old0(j,i)
1659 d_a_old(j,i)=d_a_old0(j,i)
1662 if (PRINT_AMTS_MSG) then
1663 write (iout,*) "acceleration/energy drift too large",amax,&
1664 epdrift," split increased to ",ntime_split," itime",itime,&
1670 "Uh-hu. Bumpy landscape. Maximum splitting number",&
1672 " already reached!!! Trying to carry on!"
1676 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1678 t_enegrad=t_enegrad+tcpu()-tt0
1680 ! Second step of the velocity Verlet algorithm
1685 else if (lang.eq.3) then
1687 call sd_verlet2_ciccotti
1689 else if (lang.eq.1) then
1694 if (rattle) call rattle2
1695 ! Backup the coordinates, velocities, and accelerations
1699 d_t_old(j,i)=d_t(j,i)
1700 d_a_old(j,i)=d_a(j,i)
1707 ! Restore the time step
1709 ! Compute long-range forces
1716 call etotal_long(energia_long)
1717 if (energia_long(0).gt.0.99e20 .or. isnan(energia_long(0))) then
1720 "Infinitied/NaNs in energia_long, Aborting MPI job."
1722 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1724 write (iout,*) "Infinitied/NaNs in energia_long, terminating."
1728 if (large.and. mod(itime,ntwe).eq.0) &
1729 call enerprint(energia_long)
1732 t_elong=t_elong+MPI_Wtime()-tt0
1734 t_elong=t_elong+tcpu()-tt0
1740 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1742 t_enegrad=t_enegrad+tcpu()-tt0
1744 ! Compute accelerations from long-range forces
1746 if (large.and. mod(itime,ntwe).eq.0) then
1747 write (iout,*) "energia_long",energia_long(0)
1748 write (iout,*) "Cartesian and internal coordinates: step 2"
1750 call pdbout(0.0d0,"cipiszcze",iout)
1752 write (iout,*) "Accelerations from long-range forces"
1754 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1755 (d_a(j,i+nres),j=1,3)
1757 write (iout,*) "Velocities, step 2"
1759 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1760 (d_t(j,i+nres),j=1,3)
1764 ! Compute the final RESPA step (increment velocities)
1765 ! write (iout,*) "*********************** RESPA fin"
1767 ! Compute the complete potential energy
1769 potEcomp(i)=energia_short(i)+energia_long(i)
1771 potE=potEcomp(0)-potEcomp(20)
1772 ! potE=energia_short(0)+energia_long(0)
1775 ! Calculate the kinetic and the total energy and the kinetic temperature
1778 ! Couple the system to Berendsen bath if needed
1779 if (tbf .and. lang.eq.0) then
1782 kinetic_T=2.0d0/(dimen3*Rb)*EK
1783 ! Backup the coordinates, velocities, and accelerations
1785 if (mod(itime,ntwe).eq.0 .and. large) then
1786 write (iout,*) "Velocities, end"
1788 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1789 (d_t(j,i+nres),j=1,3)
1794 end subroutine RESPA_step
1795 !-----------------------------------------------------------------------------
1796 subroutine RESPA_vel
1797 ! First and last RESPA step (incrementing velocities using long-range
1800 ! implicit real*8 (a-h,o-z)
1801 ! include 'DIMENSIONS'
1802 ! include 'COMMON.CONTROL'
1803 ! include 'COMMON.VAR'
1804 ! include 'COMMON.MD'
1805 ! include 'COMMON.CHAIN'
1806 ! include 'COMMON.DERIV'
1807 ! include 'COMMON.GEO'
1808 ! include 'COMMON.LOCAL'
1809 ! include 'COMMON.INTERACT'
1810 ! include 'COMMON.IOUNITS'
1811 ! include 'COMMON.NAMES'
1812 integer :: i,j,inres,mnum
1815 d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
1819 d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
1824 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
1825 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
1826 .and.(mnum.ne.5)) then
1829 d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
1834 end subroutine RESPA_vel
1835 !-----------------------------------------------------------------------------
1837 ! Applying velocity Verlet algorithm - step 1 to coordinates
1839 ! implicit real*8 (a-h,o-z)
1840 ! include 'DIMENSIONS'
1841 ! include 'COMMON.CONTROL'
1842 ! include 'COMMON.VAR'
1843 ! include 'COMMON.MD'
1844 ! include 'COMMON.CHAIN'
1845 ! include 'COMMON.DERIV'
1846 ! include 'COMMON.GEO'
1847 ! include 'COMMON.LOCAL'
1848 ! include 'COMMON.INTERACT'
1849 ! include 'COMMON.IOUNITS'
1850 ! include 'COMMON.NAMES'
1851 real(kind=8) :: adt,adt2
1852 integer :: i,j,inres,mnum
1855 write (iout,*) "VELVERLET1 START: DC"
1857 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),&
1858 (dc(j,i+nres),j=1,3)
1862 adt=d_a_old(j,0)*d_time
1864 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1865 d_t_new(j,0)=d_t_old(j,0)+adt2
1866 d_t(j,0)=d_t_old(j,0)+adt
1870 adt=d_a_old(j,i)*d_time
1872 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1873 d_t_new(j,i)=d_t_old(j,i)+adt2
1874 d_t(j,i)=d_t_old(j,i)+adt
1879 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
1880 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
1881 .and.(mnum.ne.5)) then
1884 adt=d_a_old(j,inres)*d_time
1886 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1887 d_t_new(j,inres)=d_t_old(j,inres)+adt2
1888 d_t(j,inres)=d_t_old(j,inres)+adt
1893 write (iout,*) "VELVERLET1 END: DC"
1895 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),&
1896 (dc(j,i+nres),j=1,3)
1900 end subroutine verlet1
1901 !-----------------------------------------------------------------------------
1903 ! Step 2 of the velocity Verlet algorithm: update velocities
1905 ! implicit real*8 (a-h,o-z)
1906 ! include 'DIMENSIONS'
1907 ! include 'COMMON.CONTROL'
1908 ! include 'COMMON.VAR'
1909 ! include 'COMMON.MD'
1910 ! include 'COMMON.CHAIN'
1911 ! include 'COMMON.DERIV'
1912 ! include 'COMMON.GEO'
1913 ! include 'COMMON.LOCAL'
1914 ! include 'COMMON.INTERACT'
1915 ! include 'COMMON.IOUNITS'
1916 ! include 'COMMON.NAMES'
1917 integer :: i,j,inres,mnum
1920 d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
1924 d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
1929 ! iti=iabs(itype(i,mnum))
1930 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
1931 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
1932 .and.(mnum.ne.5)) then
1935 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
1940 end subroutine verlet2
1941 !-----------------------------------------------------------------------------
1942 subroutine sddir_precalc
1943 ! Applying velocity Verlet algorithm - step 1 to coordinates
1944 ! implicit real*8 (a-h,o-z)
1945 ! include 'DIMENSIONS'
1951 ! include 'COMMON.CONTROL'
1952 ! include 'COMMON.VAR'
1953 ! include 'COMMON.MD'
1955 ! include 'COMMON.LANGEVIN'
1957 ! include 'COMMON.LANGEVIN.lang0'
1959 ! include 'COMMON.CHAIN'
1960 ! include 'COMMON.DERIV'
1961 ! include 'COMMON.GEO'
1962 ! include 'COMMON.LOCAL'
1963 ! include 'COMMON.INTERACT'
1964 ! include 'COMMON.IOUNITS'
1965 ! include 'COMMON.NAMES'
1966 ! include 'COMMON.TIME1'
1967 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1968 !el common /stochcalc/ stochforcvec
1969 real(kind=8) :: time00
1971 ! Compute friction and stochastic forces
1976 time_fric=time_fric+MPI_Wtime()-time00
1978 call stochastic_force(stochforcvec)
1979 time_stoch=time_stoch+MPI_Wtime()-time00
1982 ! Compute the acceleration due to friction forces (d_af_work) and stochastic
1983 ! forces (d_as_work)
1985 call ginv_mult(fric_work, d_af_work)
1986 call ginv_mult(stochforcvec, d_as_work)
1988 end subroutine sddir_precalc
1989 !-----------------------------------------------------------------------------
1990 subroutine sddir_verlet1
1991 ! Applying velocity Verlet algorithm - step 1 to velocities
1994 ! implicit real*8 (a-h,o-z)
1995 ! include 'DIMENSIONS'
1996 ! include 'COMMON.CONTROL'
1997 ! include 'COMMON.VAR'
1998 ! include 'COMMON.MD'
2000 ! include 'COMMON.LANGEVIN'
2002 ! include 'COMMON.LANGEVIN.lang0'
2004 ! include 'COMMON.CHAIN'
2005 ! include 'COMMON.DERIV'
2006 ! include 'COMMON.GEO'
2007 ! include 'COMMON.LOCAL'
2008 ! include 'COMMON.INTERACT'
2009 ! include 'COMMON.IOUNITS'
2010 ! include 'COMMON.NAMES'
2011 ! Revised 3/31/05 AL: correlation between random contributions to
2012 ! position and velocity increments included.
2013 real(kind=8) :: sqrt13 = 0.57735026918962576451d0 ! 1/sqrt(3)
2014 real(kind=8) :: adt,adt2
2015 integer :: i,j,ind,inres,mnum
2017 ! Add the contribution from BOTH friction and stochastic force to the
2018 ! coordinates, but ONLY the contribution from the friction forces to velocities
2021 adt=(d_a_old(j,0)+d_af_work(j))*d_time
2022 adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
2023 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
2024 d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
2025 d_t(j,0)=d_t_old(j,0)+adt
2030 adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
2031 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
2032 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
2033 d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
2034 d_t(j,i)=d_t_old(j,i)+adt
2040 ! iti=iabs(itype(i,mnum))
2041 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2042 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2043 .and.(mnum.ne.5)) then
2046 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
2047 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
2048 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
2049 d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
2050 d_t(j,inres)=d_t_old(j,inres)+adt
2056 end subroutine sddir_verlet1
2057 !-----------------------------------------------------------------------------
2058 subroutine sddir_verlet2
2059 ! Calculating the adjusted velocities for accelerations
2062 ! implicit real*8 (a-h,o-z)
2063 ! include 'DIMENSIONS'
2064 ! include 'COMMON.CONTROL'
2065 ! include 'COMMON.VAR'
2066 ! include 'COMMON.MD'
2068 ! include 'COMMON.LANGEVIN'
2070 ! include 'COMMON.LANGEVIN.lang0'
2072 ! include 'COMMON.CHAIN'
2073 ! include 'COMMON.DERIV'
2074 ! include 'COMMON.GEO'
2075 ! include 'COMMON.LOCAL'
2076 ! include 'COMMON.INTERACT'
2077 ! include 'COMMON.IOUNITS'
2078 ! include 'COMMON.NAMES'
2079 real(kind=8),dimension(6*nres) :: stochforcvec,d_as_work1 !(MAXRES6) maxres6=6*maxres
2080 real(kind=8) :: cos60 = 0.5d0, sin60 = 0.86602540378443864676d0
2081 integer :: i,j,ind,inres,mnum
2082 ! Revised 3/31/05 AL: correlation between random contributions to
2083 ! position and velocity increments included.
2084 ! The correlation coefficients are calculated at low-friction limit.
2085 ! Also, friction forces are now not calculated with new velocities.
2087 ! call friction_force
2088 call stochastic_force(stochforcvec)
2090 ! Compute the acceleration due to friction forces (d_af_work) and stochastic
2091 ! forces (d_as_work)
2093 call ginv_mult(stochforcvec, d_as_work1)
2099 d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j)) &
2100 +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
2105 d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j)) &
2106 +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
2112 ! iti=iabs(itype(i,mnum))
2113 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2114 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2115 .and.(mnum.ne.5)) then
2118 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres) &
2119 +d_af_work(ind+j))+sin60*d_as_work(ind+j) &
2120 +cos60*d_as_work1(ind+j))*d_time
2126 end subroutine sddir_verlet2
2127 !-----------------------------------------------------------------------------
2128 subroutine max_accel
2130 ! Find the maximum difference in the accelerations of the the sites
2131 ! at the beginning and the end of the time step.
2134 ! implicit real*8 (a-h,o-z)
2135 ! include 'DIMENSIONS'
2136 ! include 'COMMON.CONTROL'
2137 ! include 'COMMON.VAR'
2138 ! include 'COMMON.MD'
2139 ! include 'COMMON.CHAIN'
2140 ! include 'COMMON.DERIV'
2141 ! include 'COMMON.GEO'
2142 ! include 'COMMON.LOCAL'
2143 ! include 'COMMON.INTERACT'
2144 ! include 'COMMON.IOUNITS'
2145 real(kind=8),dimension(3) :: aux,accel,accel_old
2146 real(kind=8) :: dacc
2150 ! aux(j)=d_a(j,0)-d_a_old(j,0)
2151 accel_old(j)=d_a_old(j,0)
2158 ! 7/3/08 changed to asymmetric difference
2160 ! accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
2161 accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
2162 accel(j)=accel(j)+0.5d0*d_a(j,i)
2163 ! if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
2164 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
2165 dacc=dabs(accel(j)-accel_old(j))
2166 ! write (iout,*) i,dacc
2167 if (dacc.gt.amax) amax=dacc
2175 accel_old(j)=d_a_old(j,0)
2180 accel_old(j)=accel_old(j)+d_a_old(j,1)
2181 accel(j)=accel(j)+d_a(j,1)
2186 ! iti=iabs(itype(i,mnum))
2187 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2188 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2189 .and.(mnum.ne.5)) then
2191 ! accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
2192 accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
2193 accel(j)=accel(j)+d_a(j,i+nres)
2197 ! if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
2198 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
2199 dacc=dabs(accel(j)-accel_old(j))
2200 ! write (iout,*) "side-chain",i,dacc
2201 if (dacc.gt.amax) amax=dacc
2205 accel_old(j)=accel_old(j)+d_a_old(j,i)
2206 accel(j)=accel(j)+d_a(j,i)
2207 ! aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
2211 end subroutine max_accel
2212 !-----------------------------------------------------------------------------
2213 subroutine predict_edrift(epdrift)
2215 ! Predict the drift of the potential energy
2218 use control_data, only: lmuca
2219 ! implicit real*8 (a-h,o-z)
2220 ! include 'DIMENSIONS'
2221 ! include 'COMMON.CONTROL'
2222 ! include 'COMMON.VAR'
2223 ! include 'COMMON.MD'
2224 ! include 'COMMON.CHAIN'
2225 ! include 'COMMON.DERIV'
2226 ! include 'COMMON.GEO'
2227 ! include 'COMMON.LOCAL'
2228 ! include 'COMMON.INTERACT'
2229 ! include 'COMMON.IOUNITS'
2230 ! include 'COMMON.MUCA'
2231 real(kind=8) :: epdrift,epdriftij
2233 ! Drift of the potential energy
2239 epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
2240 if (lmuca) epdriftij=epdriftij*factor
2241 ! write (iout,*) "back",i,j,epdriftij
2242 if (epdriftij.gt.epdrift) epdrift=epdriftij
2246 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1.and.&
2247 molnum(i).ne.5) then
2250 dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
2251 if (lmuca) epdriftij=epdriftij*factor
2252 ! write (iout,*) "side",i,j,epdriftij
2253 if (epdriftij.gt.epdrift) epdrift=epdriftij
2257 epdrift=0.5d0*epdrift*d_time*d_time
2258 ! write (iout,*) "epdrift",epdrift
2260 end subroutine predict_edrift
2261 !-----------------------------------------------------------------------------
2262 subroutine verlet_bath
2264 ! Coupling to the thermostat by using the Berendsen algorithm
2267 ! implicit real*8 (a-h,o-z)
2268 ! include 'DIMENSIONS'
2269 ! include 'COMMON.CONTROL'
2270 ! include 'COMMON.VAR'
2271 ! include 'COMMON.MD'
2272 ! include 'COMMON.CHAIN'
2273 ! include 'COMMON.DERIV'
2274 ! include 'COMMON.GEO'
2275 ! include 'COMMON.LOCAL'
2276 ! include 'COMMON.INTERACT'
2277 ! include 'COMMON.IOUNITS'
2278 ! include 'COMMON.NAMES'
2279 real(kind=8) :: T_half,fact
2280 integer :: i,j,inres,mnum
2282 T_half=2.0d0/(dimen3*Rb)*EK
2283 fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
2284 ! write(iout,*) "T_half", T_half
2285 ! write(iout,*) "EK", EK
2286 ! write(iout,*) "fact", fact
2288 d_t(j,0)=fact*d_t(j,0)
2292 d_t(j,i)=fact*d_t(j,i)
2297 ! iti=iabs(itype(i,mnum))
2298 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2299 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2300 .and.(mnum.ne.5)) then
2303 d_t(j,inres)=fact*d_t(j,inres)
2308 end subroutine verlet_bath
2309 !-----------------------------------------------------------------------------
2311 ! Set up the initial conditions of a MD simulation
2314 use control, only:tcpu
2315 !el use io_basic, only:ilen
2318 use minimm, only:minim_dc,minimize,sc_move
2319 use io_config, only:readrst
2320 use io, only:statout
2321 ! implicit real*8 (a-h,o-z)
2322 ! include 'DIMENSIONS'
2325 character(len=16) :: form
2326 integer :: IERROR,ERRCODE
2328 integer :: iranmin,itrial,itmp
2329 ! include 'COMMON.SETUP'
2330 ! include 'COMMON.CONTROL'
2331 ! include 'COMMON.VAR'
2332 ! include 'COMMON.MD'
2334 ! include 'COMMON.LANGEVIN'
2336 ! include 'COMMON.LANGEVIN.lang0'
2338 ! include 'COMMON.CHAIN'
2339 ! include 'COMMON.DERIV'
2340 ! include 'COMMON.GEO'
2341 ! include 'COMMON.LOCAL'
2342 ! include 'COMMON.INTERACT'
2343 ! include 'COMMON.IOUNITS'
2344 ! include 'COMMON.NAMES'
2345 ! include 'COMMON.REMD'
2346 real(kind=8),dimension(0:n_ene) :: energia_long,energia_short
2347 real(kind=8),dimension(3) :: vcm,incr,L
2348 real(kind=8) :: xv,sigv,lowb,highb
2349 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
2350 character(len=256) :: qstr
2353 character(len=50) :: tytul
2354 logical :: file_exist
2355 !el common /gucio/ cm
2356 integer :: i,j,ipos,iq,iw,nft_sc,iretcode,nfun,itime,ierr,mnum
2357 real(kind=8) :: etot,tt0
2361 ! write(iout,*) "d_time", d_time
2362 ! Compute the standard deviations of stochastic forces for Langevin dynamics
2363 ! if the friction coefficients do not depend on surface area
2364 if (lang.gt.0 .and. .not.surfarea) then
2367 stdforcp(i)=stdfp(mnum)*dsqrt(gamp(mnum))
2371 stdforcsc(i)=stdfsc(iabs(itype(i,mnum)),mnum) &
2372 *dsqrt(gamsc(iabs(itype(i,mnum)),mnum))
2375 ! Open the pdb file for snapshotshots
2378 if (ilen(tmpdir).gt.0) &
2379 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2380 liczba(:ilen(liczba))//".pdb")
2382 file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2386 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) &
2387 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2388 liczba(:ilen(liczba))//".x")
2389 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2392 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) &
2393 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2394 liczba(:ilen(liczba))//".cx")
2395 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2401 if (ilen(tmpdir).gt.0) &
2402 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
2403 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
2405 if (ilen(tmpdir).gt.0) &
2406 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
2407 cartname=prefix(:ilen(prefix))//"_MD.cx"
2411 write (qstr,'(256(1h ))')
2414 iq = qinfrag(i,iset)*10
2415 iw = wfrag(i,iset)/100
2417 if(me.eq.king.or..not.out1file) &
2418 write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
2419 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
2424 iq = qinpair(i,iset)*10
2425 iw = wpair(i,iset)/100
2427 if(me.eq.king.or..not.out1file) &
2428 write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
2429 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
2433 ! pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
2435 ! cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
2437 ! cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
2439 ! statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
2443 if (restart1file) then
2445 inquire(file=mremd_rst_name,exist=file_exist)
2446 write (*,*) me," Before broadcast: file_exist",file_exist
2448 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,&
2451 write (*,*) me," After broadcast: file_exist",file_exist
2452 ! inquire(file=mremd_rst_name,exist=file_exist)
2453 if(me.eq.king.or..not.out1file) &
2454 write(iout,*) "Initial state read by master and distributed"
2456 if (ilen(tmpdir).gt.0) &
2457 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_' &
2458 //liczba(:ilen(liczba))//'.rst')
2459 inquire(file=rest2name,exist=file_exist)
2462 if(.not.restart1file) then
2463 if(me.eq.king.or..not.out1file) &
2464 write(iout,*) "Initial state will be read from file ",&
2465 rest2name(:ilen(rest2name))
2468 call rescale_weights(t_bath)
2470 if(me.eq.king.or..not.out1file)then
2471 if (restart1file) then
2472 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),&
2475 write(iout,*) "File ",rest2name(:ilen(rest2name)),&
2478 write(iout,*) "Initial velocities randomly generated"
2485 ! Generate initial velocities
2486 if(me.eq.king.or..not.out1file) &
2487 write(iout,*) "Initial velocities randomly generated"
2492 ! rest2name = prefix(:ilen(prefix))//'.rst'
2493 if(me.eq.king.or..not.out1file)then
2494 write (iout,*) "Initial velocities"
2496 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2497 (d_t(j,i+nres),j=1,3)
2499 ! Zeroing the total angular momentum of the system
2500 write(iout,*) "Calling the zero-angular momentum subroutine"
2503 ! Getting the potential energy and forces and velocities and accelerations
2505 ! write (iout,*) "velocity of the center of the mass:"
2506 ! write (iout,*) (vcm(j),j=1,3)
2508 d_t(j,0)=d_t(j,0)-vcm(j)
2510 ! Removing the velocity of the center of mass
2512 if(me.eq.king.or..not.out1file)then
2513 write (iout,*) "vcm right after adjustment:"
2514 write (iout,*) (vcm(j),j=1,3)
2518 if(iranconf.ne.0 .or.indpdb.gt.0.and..not.unres_pdb .or.preminim) then
2520 print *, 'Calling OVERLAP_SC'
2521 call overlap_sc(fail)
2522 print *,'after OVERLAP'
2525 print *,'call SC_MOVE'
2526 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
2527 print *,'SC_move',nft_sc,etot
2528 if(me.eq.king.or..not.out1file) &
2529 write(iout,*) 'SC_move',nft_sc,etot
2533 print *, 'Calling MINIM_DC'
2534 call minim_dc(etot,iretcode,nfun)
2536 call geom_to_var(nvar,varia)
2537 print *,'Calling MINIMIZE.'
2538 call minimize(etot,varia,iretcode,nfun)
2539 call var_to_geom(nvar,varia)
2541 if(me.eq.king.or..not.out1file) &
2542 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2544 if(iranconf.ne.0) then
2545 !c 8/22/17 AL Loop to produce a low-energy random conformation
2548 if(me.eq.king.or..not.out1file) &
2549 write (iout,*) 'Calling OVERLAP_SC'
2550 call overlap_sc(fail)
2551 endif !endif overlap
2554 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
2555 print *,'SC_move',nft_sc,etot
2556 if(me.eq.king.or..not.out1file) &
2557 write(iout,*) 'SC_move',nft_sc,etot
2561 print *, 'Calling MINIM_DC'
2562 call minim_dc(etot,iretcode,nfun)
2563 call int_from_cart1(.false.)
2565 call geom_to_var(nvar,varia)
2566 print *,'Calling MINIMIZE.'
2567 call minimize(etot,varia,iretcode,nfun)
2568 call var_to_geom(nvar,varia)
2570 if(me.eq.king.or..not.out1file) &
2571 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2573 if (isnan(etot) .or. etot.gt.4.0d6) then
2574 write (iout,*) "Energy too large",etot, &
2575 " trying another random conformation"
2578 call gen_rand_conf(itmp,*30)
2580 30 write (iout,*) 'Failed to generate random conformation', &
2582 write (*,*) 'Processor:',me, &
2583 ' Failed to generate random conformation',&
2592 write (iout,'(a,i3,a)') 'Processor:',me, &
2593 ' error in generating random conformation.'
2594 write (*,'(a,i3,a)') 'Processor:',me, &
2595 ' error in generating random conformation.'
2598 ! call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
2599 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
2609 write (iout,'(a,i3,a)') 'Processor:',me, &
2610 ' failed to generate a low-energy random conformation.'
2611 write (*,'(a,i3,a,f10.3)') 'Processor:',me, &
2612 ' failed to generate a low-energy random conformation.',etot
2616 ! call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
2617 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
2624 call chainbuild_cart
2629 kinetic_T=2.0d0/(dimen3*Rb)*EK
2630 if(me.eq.king.or..not.out1file)then
2640 write(iout,*) "before ETOTAL"
2641 call etotal(potEcomp)
2642 if (large) call enerprint(potEcomp)
2645 t_etotal=t_etotal+MPI_Wtime()-tt0
2647 t_etotal=t_etotal+tcpu()-tt0
2654 if (amax*d_time .gt. dvmax) then
2655 d_time=d_time*dvmax/amax
2656 if(me.eq.king.or..not.out1file) write (iout,*) &
2657 "Time step reduced to",d_time,&
2658 " because of too large initial acceleration."
2660 if(me.eq.king.or..not.out1file)then
2661 write(iout,*) "Potential energy and its components"
2662 call enerprint(potEcomp)
2663 ! write(iout,*) (potEcomp(i),i=0,n_ene)
2665 potE=potEcomp(0)-potEcomp(20)
2668 if (ntwe.ne.0) call statout(itime)
2669 if(me.eq.king.or..not.out1file) &
2670 write (iout,'(/a/3(a25,1pe14.5/))') "Initial:", &
2671 " Kinetic energy",EK," Potential energy",potE, &
2672 " Total energy",totE," Maximum acceleration ", &
2675 write (iout,*) "Initial coordinates"
2677 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),&
2680 write (iout,*) "Initial dC"
2682 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),&
2683 (dc(j,i+nres),j=1,3)
2685 write (iout,*) "Initial velocities"
2686 write (iout,"(13x,' backbone ',23x,' side chain')")
2688 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2689 (d_t(j,i+nres),j=1,3)
2691 write (iout,*) "Initial accelerations"
2693 ! write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2694 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),&
2695 (d_a(j,i+nres),j=1,3)
2701 d_t_old(j,i)=d_t(j,i)
2702 d_a_old(j,i)=d_a(j,i)
2704 ! write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
2713 call etotal_short(energia_short)
2714 if (large) call enerprint(potEcomp)
2717 t_eshort=t_eshort+MPI_Wtime()-tt0
2719 t_eshort=t_eshort+tcpu()-tt0
2724 if(.not.out1file .and. large) then
2725 write (iout,*) "energia_long",energia_long(0),&
2726 " energia_short",energia_short(0),&
2727 " total",energia_long(0)+energia_short(0)
2728 write (iout,*) "Initial fast-force accelerations"
2730 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
2731 (d_a(j,i+nres),j=1,3)
2734 ! 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
2737 d_a_short(j,i)=d_a(j,i)
2746 call etotal_long(energia_long)
2747 if (large) call enerprint(potEcomp)
2750 t_elong=t_elong+MPI_Wtime()-tt0
2752 t_elong=t_elong+tcpu()-tt0
2757 if(.not.out1file .and. large) then
2758 write (iout,*) "energia_long",energia_long(0)
2759 write (iout,*) "Initial slow-force accelerations"
2761 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
2762 (d_a(j,i+nres),j=1,3)
2766 t_enegrad=t_enegrad+MPI_Wtime()-tt0
2768 t_enegrad=t_enegrad+tcpu()-tt0
2772 end subroutine init_MD
2773 !-----------------------------------------------------------------------------
2774 subroutine random_vel
2776 ! implicit real*8 (a-h,o-z)
2778 use random, only:anorm_distr
2780 ! include 'DIMENSIONS'
2781 ! include 'COMMON.CONTROL'
2782 ! include 'COMMON.VAR'
2783 ! include 'COMMON.MD'
2785 ! include 'COMMON.LANGEVIN'
2787 ! include 'COMMON.LANGEVIN.lang0'
2789 ! include 'COMMON.CHAIN'
2790 ! include 'COMMON.DERIV'
2791 ! include 'COMMON.GEO'
2792 ! include 'COMMON.LOCAL'
2793 ! include 'COMMON.INTERACT'
2794 ! include 'COMMON.IOUNITS'
2795 ! include 'COMMON.NAMES'
2796 ! include 'COMMON.TIME1'
2797 real(kind=8) :: xv,sigv,lowb,highb ,Ek1
2800 real(kind=8) ,allocatable, dimension(:) :: DDU1,DDU2,DL2,DL1,xsolv,DML,rs
2801 real(kind=8) :: sumx
2803 real(kind=8) ,allocatable, dimension(:) :: rsold
2804 real (kind=8),allocatable,dimension(:,:) :: matold
2808 integer :: i,j,ii,k,ind,mark,imark,mnum
2809 ! Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
2810 ! First generate velocities in the eigenspace of the G matrix
2811 ! write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
2818 sigv=dsqrt((Rb*t_bath)/geigen(i))
2821 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
2823 write (iout,*) "i",i," ii",ii," geigen",geigen(i),&
2824 " d_t_work_new",d_t_work_new(ii)
2835 Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
2838 write (iout,*) "Ek from eigenvectors",Ek1
2839 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
2843 ! Transform velocities to UNRES coordinate space
2844 allocate (DL1(2*nres))
2845 allocate (DDU1(2*nres))
2846 allocate (DL2(2*nres))
2847 allocate (DDU2(2*nres))
2848 allocate (xsolv(2*nres))
2849 allocate (DML(2*nres))
2850 allocate (rs(2*nres))
2852 allocate (rsold(2*nres))
2853 allocate (matold(2*nres,2*nres))
2855 matold(1,1)=DMorig(1)
2856 matold(1,2)=DU1orig(1)
2857 matold(1,3)=DU2orig(1)
2858 write (*,*) DMorig(1),DU1orig(1),DU2orig(1)
2863 matold(i,j)=DMorig(i)
2864 matold(i,j-1)=DU1orig(i-1)
2866 matold(i,j-2)=DU2orig(i-2)
2874 matold(i,j+1)=DU1orig(i)
2880 matold(i,j+2)=DU2orig(i)
2884 matold(dimen,dimen)=DMorig(dimen)
2885 matold(dimen,dimen-1)=DU1orig(dimen-1)
2886 matold(dimen,dimen-2)=DU2orig(dimen-2)
2887 write(iout,*) "old gmatrix"
2888 call matout(dimen,dimen,2*nres,2*nres,matold)
2892 ! Find the ith eigenvector of the pentadiagonal inertiq matrix
2896 DML(j)=DMorig(j)-geigen(i)
2899 DML(j-1)=DMorig(j)-geigen(i)
2904 DDU1(imark-1)=DU2orig(imark-1)
2905 do j=imark+1,dimen-1
2906 DDU1(j-1)=DU1orig(j)
2914 DDU2(j)=DU2orig(j+1)
2923 write (iout,*) "DL2,DL1,DML,DDU1,DDU2"
2924 write(iout,'(10f10.5)') (DL2(k),k=3,dimen-1)
2925 write(iout,'(10f10.5)') (DL1(k),k=2,dimen-1)
2926 write(iout,'(10f10.5)') (DML(k),k=1,dimen-1)
2927 write(iout,'(10f10.5)') (DDU1(k),k=1,dimen-2)
2928 write(iout,'(10f10.5)') (DDU2(k),k=1,dimen-3)
2931 if (imark.gt.2) rs(imark-2)=-DU2orig(imark-2)
2932 if (imark.gt.1) rs(imark-1)=-DU1orig(imark-1)
2933 if (imark.lt.dimen) rs(imark)=-DU1orig(imark)
2934 if (imark.lt.dimen-1) rs(imark+1)=-DU2orig(imark)
2938 ! write (iout,*) "Vector rs"
2940 ! write (iout,*) j,rs(j)
2943 call FDIAG(dimen-1,DL2,DL1,DML,DDU1,DDU2,rs,xsolv,mark)
2950 sumx=-geigen(i)*xsolv(j)
2952 sumx=sumx+matold(j,k)*xsolv(k)
2955 sumx=sumx+matold(j,k)*xsolv(k-1)
2957 write(iout,'(i5,3f10.5)') j,sumx,rsold(j),sumx-rsold(j)
2960 sumx=-geigen(i)*xsolv(j-1)
2962 sumx=sumx+matold(j,k)*xsolv(k)
2965 sumx=sumx+matold(j,k)*xsolv(k-1)
2967 write(iout,'(i5,3f10.5)') j-1,sumx,rsold(j-1),sumx-rsold(j-1)
2971 "Solution of equations system",i," for eigenvalue",geigen(i)
2973 write(iout,'(i5,f10.5)') j,xsolv(j)
2976 do j=dimen-1,imark,-1
2981 write (iout,*) "Un-normalized eigenvector",i," for eigenvalue",geigen(i)
2983 write(iout,'(i5,f10.5)') j,xsolv(j)
2986 ! Normalize ith eigenvector
2989 sumx=sumx+xsolv(j)**2
2993 xsolv(j)=xsolv(j)/sumx
2996 write (iout,*) "Eigenvector",i," for eigenvalue",geigen(i)
2998 write(iout,'(i5,3f10.5)') j,xsolv(j)
3001 ! All done at this point for eigenvector i, exit loop
3009 write (iout,*) "Unable to find eigenvector",i
3012 ! write (iout,*) "i=",i
3014 ! write (iout,*) "k=",k
3017 ! write(iout,*) "ind",ind," ind1",3*(i-1)+k
3018 d_t_work(ind)=d_t_work(ind) &
3019 +xsolv(j)*d_t_work_new(3*(i-1)+k)
3022 enddo ! i (loop over the eigenvectors)
3025 write (iout,*) "d_t_work"
3027 write (iout,"(i5,f10.5)") i,d_t_work(i)
3032 ! if (itype(i,1).eq.10) then
3034 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
3035 iti=iabs(itype(i,mnum))
3036 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3037 if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
3038 .or.(mnum.eq.5)) then
3045 ! write (iout,*) "k",k," ii+k",ii+k," ii+j+k",ii+j+k,"EK1",Ek1
3046 Ek1=Ek1+0.5d0*mp(mnum)*((d_t_work(ii+j+k)+d_t_work_new(ii+k))/2)**2+&
3047 0.5d0*0.25d0*IP(mnum)*(d_t_work(ii+j+k)-d_t_work_new(ii+k))**2
3050 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3051 .and.(mnum.ne.5)) ii=ii+3
3052 write (iout,*) "i",i," itype",itype(i,mnum)," mass",msc(itype(i,mnum),mnum)
3053 write (iout,*) "ii",ii
3056 write (iout,*) "k",k," ii",ii,"EK1",EK1
3057 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3059 Ek1=Ek1+0.5d0*Isc(iabs(itype(i,mnum)),mnum)*(d_t_work(ii)-d_t_work(ii-3))**2
3060 Ek1=Ek1+0.5d0*msc(iabs(itype(i,mnum)),mnum)*d_t_work(ii)**2
3062 write (iout,*) "i",i," ii",ii
3064 write (iout,*) "Ek from d_t_work",Ek1
3065 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
3067 if(allocated(DDU1)) deallocate(DDU1)
3068 if(allocated(DDU2)) deallocate(DDU2)
3069 if(allocated(DL2)) deallocate(DL2)
3070 if(allocated(DL1)) deallocate(DL1)
3071 if(allocated(xsolv)) deallocate(xsolv)
3072 if(allocated(DML)) deallocate(DML)
3073 if(allocated(rs)) deallocate(rs)
3075 if(allocated(matold)) deallocate(matold)
3076 if(allocated(rsold)) deallocate(rsold)
3081 d_t(k,j)=d_t_work(ind)
3085 if (itype(j,1).ne.10 .and. itype(j,mnum).ne.ntyp1_molec(mnum)&
3086 .and.(mnum.ne.5)) then
3088 d_t(k,j+nres)=d_t_work(ind)
3094 write (iout,*) "Random velocities in the Calpha,SC space"
3096 write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3)
3099 write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3)
3106 ! if (itype(i,1).eq.10) then
3108 if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
3109 .or.(mnum.eq.5)) then
3111 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
3115 d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
3116 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
3121 write (iout,*)"Random velocities in the virtual-bond-vector space"
3123 write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3)
3126 write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3)
3129 write (iout,*) "Ek from d_t_work",Ek1
3130 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
3138 d_t_work(ind)=d_t_work(ind) &
3139 +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
3141 ! write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
3145 ! Transfer to the d_t vector
3147 d_t(j,0)=d_t_work(j)
3153 d_t(j,i)=d_t_work(ind)
3158 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3159 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3160 .and.(mnum.ne.5)) then
3163 d_t(j,i+nres)=d_t_work(ind)
3169 ! write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",&
3170 ! 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
3172 ! write(iout,*) "end init MD"
3174 end subroutine random_vel
3175 !-----------------------------------------------------------------------------
3177 subroutine sd_verlet_p_setup
3178 ! Sets up the parameters of stochastic Verlet algorithm
3179 ! implicit real*8 (a-h,o-z)
3180 ! include 'DIMENSIONS'
3181 use control, only: tcpu
3186 ! include 'COMMON.CONTROL'
3187 ! include 'COMMON.VAR'
3188 ! include 'COMMON.MD'
3190 ! include 'COMMON.LANGEVIN'
3192 ! include 'COMMON.LANGEVIN.lang0'
3194 ! include 'COMMON.CHAIN'
3195 ! include 'COMMON.DERIV'
3196 ! include 'COMMON.GEO'
3197 ! include 'COMMON.LOCAL'
3198 ! include 'COMMON.INTERACT'
3199 ! include 'COMMON.IOUNITS'
3200 ! include 'COMMON.NAMES'
3201 ! include 'COMMON.TIME1'
3202 real(kind=8),dimension(6*nres) :: emgdt !(MAXRES6) maxres6=6*maxres
3203 real(kind=8) :: pterm,vterm,rho,rhoc,vsig
3204 real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
3205 prand_vec,vrand_vec1,vrand_vec2 !(MAXRES6) maxres6=6*maxres
3206 logical :: lprn = .false.
3207 real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
3208 real(kind=8) :: ktm,gdt,egdt,gdt2,gdt3,gdt4,gdt5,gdt6,gdt7,gdt8,&
3210 integer :: i,maxres2
3217 ! AL 8/17/04 Code adapted from tinker
3219 ! Get the frictional and random terms for stochastic dynamics in the
3220 ! eigenspace of mass-scaled UNRES friction matrix
3224 gdt = fricgam(i) * d_time
3226 ! Stochastic dynamics reduces to simple MD for zero friction
3228 if (gdt .le. zero) then
3229 pfric_vec(i) = 1.0d0
3230 vfric_vec(i) = d_time
3231 afric_vec(i) = 0.5d0 * d_time * d_time
3232 prand_vec(i) = 0.0d0
3233 vrand_vec1(i) = 0.0d0
3234 vrand_vec2(i) = 0.0d0
3236 ! Analytical expressions when friction coefficient is large
3239 if (gdt .ge. gdt_radius) then
3242 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
3243 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
3244 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
3245 vterm = 1.0d0 - egdt**2
3246 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
3248 ! Use series expansions when friction coefficient is small
3259 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0 &
3260 - gdt5/120.0d0 + gdt6/720.0d0 &
3261 - gdt7/5040.0d0 + gdt8/40320.0d0 &
3262 - gdt9/362880.0d0) / fricgam(i)**2
3263 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
3264 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
3265 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0 &
3266 + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0 &
3267 + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0 &
3268 + 127.0d0*gdt9/90720.0d0
3269 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0 &
3270 - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0 &
3271 - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0 &
3272 - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
3273 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0 &
3274 - 17.0d0*gdt2/1280.0d0 &
3275 + 17.0d0*gdt3/6144.0d0 &
3276 + 40967.0d0*gdt4/34406400.0d0 &
3277 - 57203.0d0*gdt5/275251200.0d0 &
3278 - 1429487.0d0*gdt6/13212057600.0d0)
3281 ! Compute the scaling factors of random terms for the nonzero friction case
3283 ktm = 0.5d0*d_time/fricgam(i)
3284 psig = dsqrt(ktm*pterm) / fricgam(i)
3285 vsig = dsqrt(ktm*vterm)
3286 rhoc = dsqrt(1.0d0 - rho*rho)
3288 vrand_vec1(i) = vsig * rho
3289 vrand_vec2(i) = vsig * rhoc
3294 "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",&
3297 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),&
3298 afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
3302 ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
3305 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
3306 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
3307 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
3308 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
3309 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
3310 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
3313 t_sdsetup=t_sdsetup+MPI_Wtime()
3315 t_sdsetup=t_sdsetup+tcpu()-tt0
3318 end subroutine sd_verlet_p_setup
3319 !-----------------------------------------------------------------------------
3320 subroutine eigtransf1(n,ndim,ab,d,c)
3324 real(kind=8) :: ab(ndim,ndim,n),c(ndim,n),d(ndim)
3330 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
3335 end subroutine eigtransf1
3336 !-----------------------------------------------------------------------------
3337 subroutine eigtransf(n,ndim,a,b,d,c)
3341 real(kind=8) :: a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
3347 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
3352 end subroutine eigtransf
3353 !-----------------------------------------------------------------------------
3354 subroutine sd_verlet1
3356 ! Applying stochastic velocity Verlet algorithm - step 1 to velocities
3358 ! implicit real*8 (a-h,o-z)
3359 ! include 'DIMENSIONS'
3360 ! include 'COMMON.CONTROL'
3361 ! include 'COMMON.VAR'
3362 ! include 'COMMON.MD'
3364 ! include 'COMMON.LANGEVIN'
3366 ! include 'COMMON.LANGEVIN.lang0'
3368 ! include 'COMMON.CHAIN'
3369 ! include 'COMMON.DERIV'
3370 ! include 'COMMON.GEO'
3371 ! include 'COMMON.LOCAL'
3372 ! include 'COMMON.INTERACT'
3373 ! include 'COMMON.IOUNITS'
3374 ! include 'COMMON.NAMES'
3375 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
3376 !el common /stochcalc/ stochforcvec
3377 logical :: lprn = .false.
3378 real(kind=8) :: ddt1,ddt2
3379 integer :: i,j,ind,inres
3381 ! write (iout,*) "dc_old"
3383 ! write (iout,'(i5,3f10.5,5x,3f10.5)')
3384 ! & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3387 dc_work(j)=dc_old(j,0)
3388 d_t_work(j)=d_t_old(j,0)
3389 d_a_work(j)=d_a_old(j,0)
3394 dc_work(ind+j)=dc_old(j,i)
3395 d_t_work(ind+j)=d_t_old(j,i)
3396 d_a_work(ind+j)=d_a_old(j,i)
3402 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3403 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3404 .and.(mnum.ne.5)) then
3406 dc_work(ind+j)=dc_old(j,i+nres)
3407 d_t_work(ind+j)=d_t_old(j,i+nres)
3408 d_a_work(ind+j)=d_a_old(j,i+nres)
3416 "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",&
3420 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
3421 vfric_mat(i,j),afric_mat(i,j),&
3422 prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3430 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) &
3431 +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3432 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3433 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3435 d_t_work_new(i)=ddt1+0.5d0*ddt2
3436 d_t_work(i)=ddt1+ddt2
3441 d_t(j,0)=d_t_work(j)
3446 dc(j,i)=dc_work(ind+j)
3447 d_t(j,i)=d_t_work(ind+j)
3453 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3454 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3455 .and.(mnum.ne.5)) then
3458 dc(j,inres)=dc_work(ind+j)
3459 d_t(j,inres)=d_t_work(ind+j)
3465 end subroutine sd_verlet1
3466 !-----------------------------------------------------------------------------
3467 subroutine sd_verlet2
3469 ! Calculating the adjusted velocities for accelerations
3471 ! implicit real*8 (a-h,o-z)
3472 ! include 'DIMENSIONS'
3473 ! include 'COMMON.CONTROL'
3474 ! include 'COMMON.VAR'
3475 ! include 'COMMON.MD'
3477 ! include 'COMMON.LANGEVIN'
3479 ! include 'COMMON.LANGEVIN.lang0'
3481 ! include 'COMMON.CHAIN'
3482 ! include 'COMMON.DERIV'
3483 ! include 'COMMON.GEO'
3484 ! include 'COMMON.LOCAL'
3485 ! include 'COMMON.INTERACT'
3486 ! include 'COMMON.IOUNITS'
3487 ! include 'COMMON.NAMES'
3488 !el real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV !(MAXRES6) maxres6=6*maxres
3489 real(kind=8),dimension(6*nres) :: stochforcvecV !(MAXRES6) maxres6=6*maxres
3490 !el common /stochcalc/ stochforcvec
3492 real(kind=8) :: ddt1,ddt2
3493 integer :: i,j,ind,inres
3494 ! Compute the stochastic forces which contribute to velocity change
3496 call stochastic_force(stochforcvecV)
3503 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3504 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+ &
3505 vrand_mat2(i,j)*stochforcvecV(j)
3507 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
3511 d_t(j,0)=d_t_work(j)
3516 d_t(j,i)=d_t_work(ind+j)
3522 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3523 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3524 .and.(mnum.ne.5)) then
3527 d_t(j,inres)=d_t_work(ind+j)
3533 end subroutine sd_verlet2
3534 !-----------------------------------------------------------------------------
3535 subroutine sd_verlet_ciccotti_setup
3537 ! Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
3539 ! implicit real*8 (a-h,o-z)
3540 ! include 'DIMENSIONS'
3541 use control, only: tcpu
3546 ! include 'COMMON.CONTROL'
3547 ! include 'COMMON.VAR'
3548 ! include 'COMMON.MD'
3550 ! include 'COMMON.LANGEVIN'
3552 ! include 'COMMON.LANGEVIN.lang0'
3554 ! include 'COMMON.CHAIN'
3555 ! include 'COMMON.DERIV'
3556 ! include 'COMMON.GEO'
3557 ! include 'COMMON.LOCAL'
3558 ! include 'COMMON.INTERACT'
3559 ! include 'COMMON.IOUNITS'
3560 ! include 'COMMON.NAMES'
3561 ! include 'COMMON.TIME1'
3562 real(kind=8),dimension(6*nres) :: emgdt !(MAXRES6) maxres6=6*maxres
3563 real(kind=8) :: pterm,vterm,rho,rhoc,vsig
3564 real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
3565 prand_vec,vrand_vec1,vrand_vec2 !(MAXRES6) maxres6=6*maxres
3566 logical :: lprn = .false.
3567 real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
3568 real(kind=8) :: ktm,gdt,egdt,tt0
3569 integer :: i,maxres2
3576 ! AL 8/17/04 Code adapted from tinker
3578 ! Get the frictional and random terms for stochastic dynamics in the
3579 ! eigenspace of mass-scaled UNRES friction matrix
3583 write (iout,*) "i",i," fricgam",fricgam(i)
3584 gdt = fricgam(i) * d_time
3586 ! Stochastic dynamics reduces to simple MD for zero friction
3588 if (gdt .le. zero) then
3589 pfric_vec(i) = 1.0d0
3590 vfric_vec(i) = d_time
3591 afric_vec(i) = 0.5d0*d_time*d_time
3592 prand_vec(i) = afric_vec(i)
3593 vrand_vec2(i) = vfric_vec(i)
3595 ! Analytical expressions when friction coefficient is large
3600 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
3601 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
3602 prand_vec(i) = afric_vec(i)
3603 vrand_vec2(i) = vfric_vec(i)
3605 ! Compute the scaling factors of random terms for the nonzero friction case
3607 ! ktm = 0.5d0*d_time/fricgam(i)
3608 ! psig = dsqrt(ktm*pterm) / fricgam(i)
3609 ! vsig = dsqrt(ktm*vterm)
3610 ! prand_vec(i) = psig*afric_vec(i)
3611 ! vrand_vec2(i) = vsig*vfric_vec(i)
3616 "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",&
3619 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),&
3620 afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
3624 ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
3626 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
3627 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
3628 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
3629 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
3630 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
3632 t_sdsetup=t_sdsetup+MPI_Wtime()
3634 t_sdsetup=t_sdsetup+tcpu()-tt0
3637 end subroutine sd_verlet_ciccotti_setup
3638 !-----------------------------------------------------------------------------
3639 subroutine sd_verlet1_ciccotti
3641 ! Applying stochastic velocity Verlet algorithm - step 1 to velocities
3642 ! implicit real*8 (a-h,o-z)
3644 ! include 'DIMENSIONS'
3648 ! include 'COMMON.CONTROL'
3649 ! include 'COMMON.VAR'
3650 ! include 'COMMON.MD'
3652 ! include 'COMMON.LANGEVIN'
3654 ! include 'COMMON.LANGEVIN.lang0'
3656 ! include 'COMMON.CHAIN'
3657 ! include 'COMMON.DERIV'
3658 ! include 'COMMON.GEO'
3659 ! include 'COMMON.LOCAL'
3660 ! include 'COMMON.INTERACT'
3661 ! include 'COMMON.IOUNITS'
3662 ! include 'COMMON.NAMES'
3663 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
3664 !el common /stochcalc/ stochforcvec
3665 logical :: lprn = .false.
3666 real(kind=8) :: ddt1,ddt2
3667 integer :: i,j,ind,inres
3668 ! write (iout,*) "dc_old"
3670 ! write (iout,'(i5,3f10.5,5x,3f10.5)')
3671 ! & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3674 dc_work(j)=dc_old(j,0)
3675 d_t_work(j)=d_t_old(j,0)
3676 d_a_work(j)=d_a_old(j,0)
3681 dc_work(ind+j)=dc_old(j,i)
3682 d_t_work(ind+j)=d_t_old(j,i)
3683 d_a_work(ind+j)=d_a_old(j,i)
3688 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3690 dc_work(ind+j)=dc_old(j,i+nres)
3691 d_t_work(ind+j)=d_t_old(j,i+nres)
3692 d_a_work(ind+j)=d_a_old(j,i+nres)
3701 "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",&
3705 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
3706 vfric_mat(i,j),afric_mat(i,j),&
3707 prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3715 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) &
3716 +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3717 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3718 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3720 d_t_work_new(i)=ddt1+0.5d0*ddt2
3721 d_t_work(i)=ddt1+ddt2
3726 d_t(j,0)=d_t_work(j)
3731 dc(j,i)=dc_work(ind+j)
3732 d_t(j,i)=d_t_work(ind+j)
3738 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3739 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3740 .and.(mnum.ne.5)) then
3743 dc(j,inres)=dc_work(ind+j)
3744 d_t(j,inres)=d_t_work(ind+j)
3750 end subroutine sd_verlet1_ciccotti
3751 !-----------------------------------------------------------------------------
3752 subroutine sd_verlet2_ciccotti
3754 ! Calculating the adjusted velocities for accelerations
3756 ! implicit real*8 (a-h,o-z)
3757 ! include 'DIMENSIONS'
3758 ! include 'COMMON.CONTROL'
3759 ! include 'COMMON.VAR'
3760 ! include 'COMMON.MD'
3762 ! include 'COMMON.LANGEVIN'
3764 ! include 'COMMON.LANGEVIN.lang0'
3766 ! include 'COMMON.CHAIN'
3767 ! include 'COMMON.DERIV'
3768 ! include 'COMMON.GEO'
3769 ! include 'COMMON.LOCAL'
3770 ! include 'COMMON.INTERACT'
3771 ! include 'COMMON.IOUNITS'
3772 ! include 'COMMON.NAMES'
3773 !el real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV !(MAXRES6) maxres6=6*maxres
3774 real(kind=8),dimension(6*nres) :: stochforcvecV !(MAXRES6) maxres6=6*maxres
3775 !el common /stochcalc/ stochforcvec
3776 real(kind=8) :: ddt1,ddt2
3777 integer :: i,j,ind,inres
3779 ! Compute the stochastic forces which contribute to velocity change
3781 call stochastic_force(stochforcvecV)
3788 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3789 ! ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
3790 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
3792 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
3796 d_t(j,0)=d_t_work(j)
3801 d_t(j,i)=d_t_work(ind+j)
3807 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3809 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3812 d_t(j,inres)=d_t_work(ind+j)
3818 end subroutine sd_verlet2_ciccotti
3820 !-----------------------------------------------------------------------------
3822 !-----------------------------------------------------------------------------
3823 subroutine inertia_tensor
3825 ! Calculating the intertia tensor for the entire protein in order to
3826 ! remove the perpendicular components of velocity matrix which cause
3827 ! the molecule to rotate.
3830 ! implicit real*8 (a-h,o-z)
3831 ! include 'DIMENSIONS'
3832 ! include 'COMMON.CONTROL'
3833 ! include 'COMMON.VAR'
3834 ! include 'COMMON.MD'
3835 ! include 'COMMON.CHAIN'
3836 ! include 'COMMON.DERIV'
3837 ! include 'COMMON.GEO'
3838 ! include 'COMMON.LOCAL'
3839 ! include 'COMMON.INTERACT'
3840 ! include 'COMMON.IOUNITS'
3841 ! include 'COMMON.NAMES'
3843 real(kind=8),dimension(3,3) :: Im,Imcp,eigvec,Id
3844 real(kind=8),dimension(3) :: pr,eigval,L,vp,vrot
3845 real(kind=8) :: M_SC,mag,mag2,M_PEP
3846 real(kind=8),dimension(3,0:nres) :: vpp !(3,0:MAXRES)
3847 real(kind=8),dimension(3) :: vs_p,pp,incr,v
3848 real(kind=8),dimension(3,3) :: pr1,pr2
3850 !el common /gucio/ cm
3851 integer :: iti,inres,i,j,k,mnum
3862 ! calculating the center of the mass of the protein
3866 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
3867 if (itype(i,mnum).eq.ntyp1_molec(mnum)) cycle
3868 M_PEP=M_PEP+mp(mnum)
3870 cm(j)=cm(j)+(c(j,i)+0.5d0*dc(j,i))*mp(mnum)
3879 if (mnum.eq.5) cycle
3880 iti=iabs(itype(i,mnum))
3881 M_SC=M_SC+msc(iabs(iti),mnum)
3884 cm(j)=cm(j)+msc(iabs(iti),mnum)*c(j,inres)
3888 cm(j)=cm(j)/(M_SC+M_PEP)
3893 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
3895 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
3897 Im(1,1)=Im(1,1)+mp(mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
3898 Im(1,2)=Im(1,2)-mp(mnum)*pr(1)*pr(2)
3899 Im(1,3)=Im(1,3)-mp(mnum)*pr(1)*pr(3)
3900 Im(2,3)=Im(2,3)-mp(mnum)*pr(2)*pr(3)
3901 Im(2,2)=Im(2,2)+mp(mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
3902 Im(3,3)=Im(3,3)+mp(mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
3907 iti=iabs(itype(i,mnum))
3908 if (mnum.eq.5) cycle
3911 pr(j)=c(j,inres)-cm(j)
3913 Im(1,1)=Im(1,1)+msc(iabs(iti),mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
3914 Im(1,2)=Im(1,2)-msc(iabs(iti),mnum)*pr(1)*pr(2)
3915 Im(1,3)=Im(1,3)-msc(iabs(iti),mnum)*pr(1)*pr(3)
3916 Im(2,3)=Im(2,3)-msc(iabs(iti),mnum)*pr(2)*pr(3)
3917 Im(2,2)=Im(2,2)+msc(iabs(iti),mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
3918 Im(3,3)=Im(3,3)+msc(iabs(iti),mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
3923 Im(1,1)=Im(1,1)+Ip(mnum)*(1-dc_norm(1,i)*dc_norm(1,i))* &
3924 vbld(i+1)*vbld(i+1)*0.25d0
3925 Im(1,2)=Im(1,2)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(2,i))* &
3926 vbld(i+1)*vbld(i+1)*0.25d0
3927 Im(1,3)=Im(1,3)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(3,i))* &
3928 vbld(i+1)*vbld(i+1)*0.25d0
3929 Im(2,3)=Im(2,3)+Ip(mnum)*(-dc_norm(2,i)*dc_norm(3,i))* &
3930 vbld(i+1)*vbld(i+1)*0.25d0
3931 Im(2,2)=Im(2,2)+Ip(mnum)*(1-dc_norm(2,i)*dc_norm(2,i))* &
3932 vbld(i+1)*vbld(i+1)*0.25d0
3933 Im(3,3)=Im(3,3)+Ip(mnum)*(1-dc_norm(3,i)*dc_norm(3,i))* &
3934 vbld(i+1)*vbld(i+1)*0.25d0
3940 ! if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
3941 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3942 .and.(mnum.ne.5)) then
3943 iti=iabs(itype(i,mnum))
3945 Im(1,1)=Im(1,1)+Isc(iti,mnum)*(1-dc_norm(1,inres)* &
3946 dc_norm(1,inres))*vbld(inres)*vbld(inres)
3947 Im(1,2)=Im(1,2)-Isc(iti,mnum)*(dc_norm(1,inres)* &
3948 dc_norm(2,inres))*vbld(inres)*vbld(inres)
3949 Im(1,3)=Im(1,3)-Isc(iti,mnum)*(dc_norm(1,inres)* &
3950 dc_norm(3,inres))*vbld(inres)*vbld(inres)
3951 Im(2,3)=Im(2,3)-Isc(iti,mnum)*(dc_norm(2,inres)* &
3952 dc_norm(3,inres))*vbld(inres)*vbld(inres)
3953 Im(2,2)=Im(2,2)+Isc(iti,mnum)*(1-dc_norm(2,inres)* &
3954 dc_norm(2,inres))*vbld(inres)*vbld(inres)
3955 Im(3,3)=Im(3,3)+Isc(iti,mnum)*(1-dc_norm(3,inres)* &
3956 dc_norm(3,inres))*vbld(inres)*vbld(inres)
3961 ! write(iout,*) "The angular momentum before adjustment:"
3962 ! write(iout,*) (L(j),j=1,3)
3968 ! Copying the Im matrix for the djacob subroutine
3976 ! Finding the eigenvectors and eignvalues of the inertia tensor
3977 call djacob(3,3,10000,1.0d-10,Imcp,eigvec,eigval)
3978 ! write (iout,*) "Eigenvalues & Eigenvectors"
3979 ! write (iout,'(5x,3f10.5)') (eigval(i),i=1,3)
3982 ! write (iout,'(i5,3f10.5)') i,(eigvec(i,j),j=1,3)
3984 ! Constructing the diagonalized matrix
3986 if (dabs(eigval(i)).gt.1.0d-15) then
3987 Id(i,i)=1.0d0/eigval(i)
3994 Imcp(i,j)=eigvec(j,i)
4000 pr1(i,j)=pr1(i,j)+Id(i,k)*Imcp(k,j)
4007 pr2(i,j)=pr2(i,j)+eigvec(i,k)*pr1(k,j)
4011 ! Calculating the total rotational velocity of the molecule
4014 vrot(i)=vrot(i)+pr2(i,j)*L(j)
4017 ! Resetting the velocities
4019 call vecpr(vrot(1),dc(1,i),vp)
4021 d_t(j,i)=d_t(j,i)-vp(j)
4026 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4027 .and.(mnum.ne.5)) then
4028 ! if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
4029 ! if(itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
4031 call vecpr(vrot(1),dc(1,inres),vp)
4033 d_t(j,inres)=d_t(j,inres)-vp(j)
4038 ! write(iout,*) "The angular momentum after adjustment:"
4039 ! write(iout,*) (L(j),j=1,3)
4042 end subroutine inertia_tensor
4043 !-----------------------------------------------------------------------------
4044 subroutine angmom(cm,L)
4047 ! implicit real*8 (a-h,o-z)
4048 ! include 'DIMENSIONS'
4049 ! include 'COMMON.CONTROL'
4050 ! include 'COMMON.VAR'
4051 ! include 'COMMON.MD'
4052 ! include 'COMMON.CHAIN'
4053 ! include 'COMMON.DERIV'
4054 ! include 'COMMON.GEO'
4055 ! include 'COMMON.LOCAL'
4056 ! include 'COMMON.INTERACT'
4057 ! include 'COMMON.IOUNITS'
4058 ! include 'COMMON.NAMES'
4059 real(kind=8) :: mscab
4060 real(kind=8),dimension(3) :: L,cm,pr,vp,vrot,incr,v,pp
4061 integer :: iti,inres,i,j,mnum
4062 ! Calculate the angular momentum
4071 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
4073 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
4076 v(j)=incr(j)+0.5d0*d_t(j,i)
4079 incr(j)=incr(j)+d_t(j,i)
4081 call vecpr(pr(1),v(1),vp)
4083 L(j)=L(j)+mp(mnum)*vp(j)
4087 pp(j)=0.5d0*d_t(j,i)
4089 call vecpr(pr(1),pp(1),vp)
4091 L(j)=L(j)+Ip(mnum)*vp(j)
4099 iti=iabs(itype(i,mnum))
4107 pr(j)=c(j,inres)-cm(j)
4109 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4110 .and.(mnum.ne.5)) then
4112 v(j)=incr(j)+d_t(j,inres)
4119 call vecpr(pr(1),v(1),vp)
4120 ! write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),&
4121 ! " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
4123 L(j)=L(j)+mscab*vp(j)
4125 ! write (iout,*) "L",(l(j),j=1,3)
4126 if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4127 .and.(mnum.ne.5)) then
4129 v(j)=incr(j)+d_t(j,inres)
4131 call vecpr(dc(1,inres),d_t(1,inres),vp)
4133 L(j)=L(j)+Isc(iti,mnum)*vp(j)
4137 incr(j)=incr(j)+d_t(j,i)
4141 end subroutine angmom
4142 !-----------------------------------------------------------------------------
4143 subroutine vcm_vel(vcm)
4146 ! implicit real*8 (a-h,o-z)
4147 ! include 'DIMENSIONS'
4148 ! include 'COMMON.VAR'
4149 ! include 'COMMON.MD'
4150 ! include 'COMMON.CHAIN'
4151 ! include 'COMMON.DERIV'
4152 ! include 'COMMON.GEO'
4153 ! include 'COMMON.LOCAL'
4154 ! include 'COMMON.INTERACT'
4155 ! include 'COMMON.IOUNITS'
4156 real(kind=8),dimension(3) :: vcm,vv
4157 real(kind=8) :: summas,amas
4167 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
4169 summas=summas+mp(mnum)
4171 vcm(j)=vcm(j)+mp(mnum)*(vv(j)+0.5d0*d_t(j,i))
4175 amas=msc(iabs(itype(i,mnum)),mnum)
4180 if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4181 .and.(mnum.ne.5)) then
4183 vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres))
4187 vcm(j)=vcm(j)+amas*vv(j)
4191 vv(j)=vv(j)+d_t(j,i)
4194 ! write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas
4196 vcm(j)=vcm(j)/summas
4199 end subroutine vcm_vel
4200 !-----------------------------------------------------------------------------
4202 !-----------------------------------------------------------------------------
4204 ! RATTLE algorithm for velocity Verlet - step 1, UNRES
4208 ! implicit real*8 (a-h,o-z)
4209 ! include 'DIMENSIONS'
4211 ! include 'COMMON.CONTROL'
4212 ! include 'COMMON.VAR'
4213 ! include 'COMMON.MD'
4215 ! include 'COMMON.LANGEVIN'
4217 ! include 'COMMON.LANGEVIN.lang0'
4219 ! include 'COMMON.CHAIN'
4220 ! include 'COMMON.DERIV'
4221 ! include 'COMMON.GEO'
4222 ! include 'COMMON.LOCAL'
4223 ! include 'COMMON.INTERACT'
4224 ! include 'COMMON.IOUNITS'
4225 ! include 'COMMON.NAMES'
4226 ! include 'COMMON.TIME1'
4227 !el real(kind=8) :: gginv(2*nres,2*nres),&
4228 !el gdc(3,2*nres,2*nres)
4229 real(kind=8) :: dC_uncor(3,2*nres) !,&
4230 !el real(kind=8) :: Cmat(2*nres,2*nres)
4231 real(kind=8) :: x(2*nres),xcorr(3,2*nres) !maxres2=2*maxres
4232 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
4233 !el common /przechowalnia/ nbond
4234 integer :: max_rattle = 5
4235 logical :: lprn = .false., lprn1 = .false., not_done
4236 real(kind=8) :: tol_rattle = 1.0d-5
4238 integer :: ii,i,j,jj,l,ind,ind1,nres2
4241 !el /common/ przechowalnia
4243 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4244 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4245 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4247 if (lprn) write (iout,*) "RATTLE1"
4251 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4252 .and.(mnum.ne.5)) nbond=nbond+1
4254 ! Make a folded form of the Ginv-matrix
4267 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4272 if (itype(k,1).ne.10 .and. itype(k,mnum).ne.ntyp1_molec(mnum)&
4273 .and.(mnum.ne.5)) then
4277 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4285 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4296 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4300 if (itype(k,1).ne.10) then
4304 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4312 write (iout,*) "Matrix GGinv"
4313 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
4319 if (iter.gt.max_rattle) then
4320 write (iout,*) "Error - too many iterations in RATTLE."
4323 ! Calculate the matrix C = GG**(-1) dC_old o dC
4328 dC_uncor(j,ind1)=dC(j,i)
4332 if (itype(i,1).ne.10) then
4335 dC_uncor(j,ind1)=dC(j,i+nres)
4344 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
4348 if (itype(k,1).ne.10) then
4351 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
4356 ! Calculate deviations from standard virtual-bond lengths
4360 x(ind)=vbld(i+1)**2-vbl**2
4363 if (itype(i,1).ne.10) then
4365 x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4369 write (iout,*) "Coordinates and violations"
4371 write(iout,'(i5,3f10.5,5x,e15.5)') &
4372 i,(dC_uncor(j,i),j=1,3),x(i)
4374 write (iout,*) "Velocities and violations"
4378 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4379 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4383 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4384 .and.(mnum.ne.5)) then
4387 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4388 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4389 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4392 ! write (iout,*) "gdc"
4394 ! write (iout,*) "i",i
4396 ! write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4402 if (dabs(x(i)).gt.xmax) then
4406 if (xmax.lt.tol_rattle) then
4410 ! Calculate the matrix of the system of equations
4415 Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
4420 write (iout,*) "Matrix Cmat"
4421 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4423 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4424 ! Add constraint term to positions
4431 xx = xx+x(ii)*gdc(j,ind,ii)
4435 d_t_new(j,i)=d_t_new(j,i)-xx/d_time
4439 if (itype(i,1).ne.10) then
4444 xx = xx+x(ii)*gdc(j,ind,ii)
4447 dC(j,i+nres)=dC(j,i+nres)-xx
4448 d_t_new(j,i+nres)=d_t_new(j,i+nres)-xx/d_time
4452 ! Rebuild the chain using the new coordinates
4453 call chainbuild_cart
4455 write (iout,*) "New coordinates, Lagrange multipliers,",&
4456 " and differences between actual and standard bond lengths"
4460 xx=vbld(i+1)**2-vbl**2
4461 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4462 i,(dC(j,i),j=1,3),x(ind),xx
4466 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4469 xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4470 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4471 i,(dC(j,i+nres),j=1,3),x(ind),xx
4474 write (iout,*) "Velocities and violations"
4478 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4479 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4482 if (itype(i,1).ne.10) then
4484 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4485 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4486 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4493 10 write (iout,*) "Error - singularity in solving the system",&
4494 " of equations for Lagrange multipliers."
4498 "RATTLE inactive; use -DRATTLE switch at compile time."
4501 end subroutine rattle1
4502 !-----------------------------------------------------------------------------
4504 ! RATTLE algorithm for velocity Verlet - step 2, UNRES
4508 ! implicit real*8 (a-h,o-z)
4509 ! include 'DIMENSIONS'
4511 ! include 'COMMON.CONTROL'
4512 ! include 'COMMON.VAR'
4513 ! include 'COMMON.MD'
4515 ! include 'COMMON.LANGEVIN'
4517 ! include 'COMMON.LANGEVIN.lang0'
4519 ! include 'COMMON.CHAIN'
4520 ! include 'COMMON.DERIV'
4521 ! include 'COMMON.GEO'
4522 ! include 'COMMON.LOCAL'
4523 ! include 'COMMON.INTERACT'
4524 ! include 'COMMON.IOUNITS'
4525 ! include 'COMMON.NAMES'
4526 ! include 'COMMON.TIME1'
4527 !el real(kind=8) :: gginv(2*nres,2*nres),&
4528 !el gdc(3,2*nres,2*nres)
4529 real(kind=8) :: dC_uncor(3,2*nres) !,&
4530 !el Cmat(2*nres,2*nres)
4531 real(kind=8) :: x(2*nres) !maxres2=2*maxres
4532 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
4533 !el common /przechowalnia/ nbond
4534 integer :: max_rattle = 5
4535 logical :: lprn = .false., lprn1 = .false., not_done
4536 real(kind=8) :: tol_rattle = 1.0d-5
4540 !el /common/ przechowalnia
4541 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4542 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4543 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4545 if (lprn) write (iout,*) "RATTLE2"
4546 if (lprn) write (iout,*) "Velocity correction"
4547 ! Calculate the matrix G dC
4553 gdc(j,i,ind)=GGinv(i,ind)*dC(j,k)
4558 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4559 .and.(mnum.ne.5)) then
4562 gdc(j,i,ind)=GGinv(i,ind)*dC(j,k+nres)
4568 ! write (iout,*) "gdc"
4570 ! write (iout,*) "i",i
4572 ! write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4576 ! Calculate the matrix of the system of equations
4583 Cmat(ind,j)=Cmat(ind,j)+dC(k,i)*gdc(k,ind,j)
4589 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4590 .and.(mnum.ne.5)) then
4595 Cmat(ind,j)=Cmat(ind,j)+dC(k,i+nres)*gdc(k,ind,j)
4600 ! Calculate the scalar product dC o d_t_new
4604 x(ind)=scalar(d_t(1,i),dC(1,i))
4608 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4609 .and.(mnum.ne.5)) then
4611 x(ind)=scalar(d_t(1,i+nres),dC(1,i+nres))
4615 write (iout,*) "Velocities and violations"
4619 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4620 i,ind,(d_t(j,i),j=1,3),x(ind)
4624 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4625 .and.(mnum.ne.5)) then
4627 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4628 i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind)
4634 if (dabs(x(i)).gt.xmax) then
4638 if (xmax.lt.tol_rattle) then
4643 write (iout,*) "Matrix Cmat"
4644 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4646 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4647 ! Add constraint term to velocities
4654 xx = xx+x(ii)*gdc(j,ind,ii)
4656 d_t(j,i)=d_t(j,i)-xx
4661 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4662 .and.(mnum.ne.5)) then
4667 xx = xx+x(ii)*gdc(j,ind,ii)
4669 d_t(j,i+nres)=d_t(j,i+nres)-xx
4675 "New velocities, Lagrange multipliers violations"
4679 if (lprn) write (iout,'(2i5,3f10.5,5x,2e15.5)') &
4680 i,ind,(d_t(j,i),j=1,3),x(ind),scalar(d_t(1,i),dC(1,i))
4684 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4687 write (iout,'(2i5,3f10.5,5x,2e15.5)') &
4688 i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind),&
4689 scalar(d_t(1,i+nres),dC(1,i+nres))
4695 10 write (iout,*) "Error - singularity in solving the system",&
4696 " of equations for Lagrange multipliers."
4700 "RATTLE inactive; use -DRATTLE option at compile time."
4703 end subroutine rattle2
4704 !-----------------------------------------------------------------------------
4705 subroutine rattle_brown
4706 ! RATTLE/LINCS algorithm for Brownian dynamics, UNRES
4710 ! implicit real*8 (a-h,o-z)
4711 ! include 'DIMENSIONS'
4713 ! include 'COMMON.CONTROL'
4714 ! include 'COMMON.VAR'
4715 ! include 'COMMON.MD'
4717 ! include 'COMMON.LANGEVIN'
4719 ! include 'COMMON.LANGEVIN.lang0'
4721 ! include 'COMMON.CHAIN'
4722 ! include 'COMMON.DERIV'
4723 ! include 'COMMON.GEO'
4724 ! include 'COMMON.LOCAL'
4725 ! include 'COMMON.INTERACT'
4726 ! include 'COMMON.IOUNITS'
4727 ! include 'COMMON.NAMES'
4728 ! include 'COMMON.TIME1'
4729 !el real(kind=8) :: gginv(2*nres,2*nres),&
4730 !el gdc(3,2*nres,2*nres)
4731 real(kind=8) :: dC_uncor(3,2*nres) !,&
4732 !el real(kind=8) :: Cmat(2*nres,2*nres)
4733 real(kind=8) :: x(2*nres) !maxres2=2*maxres
4734 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
4735 !el common /przechowalnia/ nbond
4736 integer :: max_rattle = 5
4737 logical :: lprn = .false., lprn1 = .false., not_done
4738 real(kind=8) :: tol_rattle = 1.0d-5
4742 !el /common/ przechowalnia
4743 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4744 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4745 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4748 if (lprn) write (iout,*) "RATTLE_BROWN"
4751 if (itype(i,1).ne.10) nbond=nbond+1
4753 ! Make a folded form of the Ginv-matrix
4766 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4770 if (itype(k,1).ne.10) then
4774 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4781 if (itype(i,1).ne.10) then
4791 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4795 if (itype(k,1).ne.10) then
4799 if (j.eq.1 .and. l.eq.1)GGinv(ii,jj)=fricmat(ind,ind1)
4807 write (iout,*) "Matrix GGinv"
4808 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
4814 if (iter.gt.max_rattle) then
4815 write (iout,*) "Error - too many iterations in RATTLE."
4818 ! Calculate the matrix C = GG**(-1) dC_old o dC
4823 dC_uncor(j,ind1)=dC(j,i)
4827 if (itype(i,1).ne.10) then
4830 dC_uncor(j,ind1)=dC(j,i+nres)
4839 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
4843 if (itype(k,1).ne.10) then
4846 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
4851 ! Calculate deviations from standard virtual-bond lengths
4855 x(ind)=vbld(i+1)**2-vbl**2
4858 if (itype(i,1).ne.10) then
4860 x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4864 write (iout,*) "Coordinates and violations"
4866 write(iout,'(i5,3f10.5,5x,e15.5)') &
4867 i,(dC_uncor(j,i),j=1,3),x(i)
4869 write (iout,*) "Velocities and violations"
4873 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4874 i,ind,(d_t(j,i),j=1,3),scalar(d_t(1,i),dC_old(1,i))
4877 if (itype(i,1).ne.10) then
4879 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4880 i+nres,ind,(d_t(j,i+nres),j=1,3),&
4881 scalar(d_t(1,i+nres),dC_old(1,i+nres))
4884 write (iout,*) "gdc"
4886 write (iout,*) "i",i
4888 write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4894 if (dabs(x(i)).gt.xmax) then
4898 if (xmax.lt.tol_rattle) then
4902 ! Calculate the matrix of the system of equations
4907 Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
4912 write (iout,*) "Matrix Cmat"
4913 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4915 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4916 ! Add constraint term to positions
4923 xx = xx+x(ii)*gdc(j,ind,ii)
4926 d_t(j,i)=d_t(j,i)+xx/d_time
4931 if (itype(i,1).ne.10) then
4936 xx = xx+x(ii)*gdc(j,ind,ii)
4939 d_t(j,i+nres)=d_t(j,i+nres)+xx/d_time
4940 dC(j,i+nres)=dC(j,i+nres)+xx
4944 ! Rebuild the chain using the new coordinates
4945 call chainbuild_cart
4947 write (iout,*) "New coordinates, Lagrange multipliers,",&
4948 " and differences between actual and standard bond lengths"
4952 xx=vbld(i+1)**2-vbl**2
4953 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4954 i,(dC(j,i),j=1,3),x(ind),xx
4957 if (itype(i,1).ne.10) then
4959 xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4960 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4961 i,(dC(j,i+nres),j=1,3),x(ind),xx
4964 write (iout,*) "Velocities and violations"
4968 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4969 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4972 if (itype(i,1).ne.10) then
4974 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4975 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4976 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4983 10 write (iout,*) "Error - singularity in solving the system",&
4984 " of equations for Lagrange multipliers."
4988 "RATTLE inactive; use -DRATTLE option at compile time"
4991 end subroutine rattle_brown
4992 !-----------------------------------------------------------------------------
4994 !-----------------------------------------------------------------------------
4995 subroutine friction_force
5000 ! implicit real*8 (a-h,o-z)
5001 ! include 'DIMENSIONS'
5002 ! include 'COMMON.VAR'
5003 ! include 'COMMON.CHAIN'
5004 ! include 'COMMON.DERIV'
5005 ! include 'COMMON.GEO'
5006 ! include 'COMMON.LOCAL'
5007 ! include 'COMMON.INTERACT'
5008 ! include 'COMMON.MD'
5010 ! include 'COMMON.LANGEVIN'
5012 ! include 'COMMON.LANGEVIN.lang0'
5014 ! include 'COMMON.IOUNITS'
5015 !el real(kind=8),dimension(6*nres) :: gamvec !(MAXRES6) maxres6=6*maxres
5016 !el common /syfek/ gamvec
5017 real(kind=8) :: vv(3),vvtot(3,nres),v_work(6*nres) !,&
5018 !el ginvfric(2*nres,2*nres) !maxres2=2*maxres
5019 !el common /przechowalnia/ ginvfric
5021 logical :: lprn = .false., checkmode = .false.
5022 integer :: i,j,ind,k,nres2,nres6,mnum
5026 if(.not.allocated(gamvec)) allocate(gamvec(nres6)) !(MAXRES6)
5027 if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
5035 d_t_work(j)=d_t(j,0)
5040 d_t_work(ind+j)=d_t(j,i)
5046 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5047 .and.(mnum.ne.5)) then
5049 d_t_work(ind+j)=d_t(j,i+nres)
5055 call fricmat_mult(d_t_work,fric_work)
5057 if (.not.checkmode) return
5060 write (iout,*) "d_t_work and fric_work"
5062 write (iout,'(i3,2e15.5)') i,d_t_work(i),fric_work(i)
5066 friction(j,0)=fric_work(j)
5071 friction(j,i)=fric_work(ind+j)
5077 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5078 .and.(mnum.ne.5)) then
5079 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
5081 friction(j,i+nres)=fric_work(ind+j)
5087 write(iout,*) "Friction backbone"
5089 write(iout,'(i5,3e15.5,5x,3e15.5)') &
5090 i,(friction(j,i),j=1,3),(d_t(j,i),j=1,3)
5092 write(iout,*) "Friction side chain"
5094 write(iout,'(i5,3e15.5,5x,3e15.5)') &
5095 i,(friction(j,i+nres),j=1,3),(d_t(j,i+nres),j=1,3)
5104 vvtot(j,i)=vv(j)+0.5d0*d_t(j,i)
5105 vvtot(j,i+nres)=vv(j)+d_t(j,i+nres)
5106 vv(j)=vv(j)+d_t(j,i)
5109 write (iout,*) "vvtot backbone and sidechain"
5111 write (iout,'(i5,3e15.5,5x,3e15.5)') i,(vvtot(j,i),j=1,3),&
5112 (vvtot(j,i+nres),j=1,3)
5117 v_work(ind+j)=vvtot(j,i)
5123 v_work(ind+j)=vvtot(j,i+nres)
5127 write (iout,*) "v_work gamvec and site-based friction forces"
5129 write (iout,'(i5,3e15.5)') i,v_work(i),gamvec(i),&
5133 ! fric_work1(i)=0.0d0
5135 ! fric_work1(i)=fric_work1(i)-A(j,i)*gamvec(j)*v_work(j)
5138 ! write (iout,*) "fric_work and fric_work1"
5140 ! write (iout,'(i5,2e15.5)') i,fric_work(i),fric_work1(i)
5146 ginvfric(i,j)=ginvfric(i,j)+ginv(i,k)*fricmat(k,j)
5150 write (iout,*) "ginvfric"
5152 write (iout,'(i5,100f8.3)') i,(ginvfric(i,j),j=1,dimen)
5154 write (iout,*) "symmetry check"
5157 write (iout,*) i,j,ginvfric(i,j)-ginvfric(j,i)
5162 end subroutine friction_force
5163 !-----------------------------------------------------------------------------
5164 subroutine setup_fricmat
5168 use control_data, only:time_Bcast
5169 use control, only:tcpu
5171 ! implicit real*8 (a-h,o-z)
5175 real(kind=8) :: time00
5177 ! include 'DIMENSIONS'
5178 ! include 'COMMON.VAR'
5179 ! include 'COMMON.CHAIN'
5180 ! include 'COMMON.DERIV'
5181 ! include 'COMMON.GEO'
5182 ! include 'COMMON.LOCAL'
5183 ! include 'COMMON.INTERACT'
5184 ! include 'COMMON.MD'
5185 ! include 'COMMON.SETUP'
5186 ! include 'COMMON.TIME1'
5187 ! integer licznik /0/
5190 ! include 'COMMON.LANGEVIN'
5192 ! include 'COMMON.LANGEVIN.lang0'
5194 ! include 'COMMON.IOUNITS'
5196 integer :: i,j,ind,ind1,m
5197 logical :: lprn = .false.
5198 real(kind=8) :: dtdi !el ,gamvec(2*nres)
5199 !el real(kind=8),dimension(2*nres,2*nres) :: ginvfric,fcopy
5200 ! real(kind=8),allocatable,dimension(:,:) :: fcopy
5201 !el real(kind=8),dimension(2*nres*(2*nres+1)/2) :: Ghalf !(mmaxres2) (mmaxres2=(maxres2*(maxres2+1)/2))
5202 !el common /syfek/ gamvec
5203 real(kind=8) :: work(8*2*nres)
5204 integer :: iwork(2*nres)
5205 !el common /przechowalnia/ ginvfric,Ghalf,fcopy
5206 integer :: ii,iti,k,l,nzero,nres2,nres6,ierr,mnum
5210 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
5211 if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
5212 if (fg_rank.ne.king) goto 10
5217 if(.not.allocated(gamvec)) allocate(gamvec(nres2)) !(MAXRES2)
5218 if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
5219 if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
5220 !el allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
5221 if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !maxres2=2*maxres
5223 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
5224 ! Zeroing out fricmat
5230 ! Load the friction coefficients corresponding to peptide groups
5235 gamvec(ind1)=gamp(mnum)
5237 ! Load the friction coefficients corresponding to side chains
5241 gamsc(ntyp1_molec(j),j)=1.0d0
5248 gamvec(ii)=gamsc(iabs(iti),mnum)
5250 if (surfarea) call sdarea(gamvec)
5252 ! write (iout,*) "Matrix A and vector gamma"
5254 ! write (iout,'(i2,$)') i
5256 ! write (iout,'(f4.1,$)') A(i,j)
5258 ! write (iout,'(f8.3)') gamvec(i)
5262 write (iout,*) "Vector gamvec"
5264 write (iout,'(i5,f10.5)') i, gamvec(i)
5268 ! The friction matrix
5273 dtdi=dtdi+A(j,k)*A(j,i)*gamvec(j)
5280 write (iout,'(//a)') "Matrix fricmat"
5281 call matout2(dimen,dimen,nres2,nres2,fricmat)
5283 if (lang.eq.2 .or. lang.eq.3) then
5284 ! Mass-scale the friction matrix if non-direct integration will be performed
5290 Ginvfric(i,j)=Ginvfric(i,j)+ &
5291 Gsqrm(i,k)*Gsqrm(l,j)*fricmat(k,l)
5296 ! Diagonalize the friction matrix
5301 Ghalf(ind)=Ginvfric(i,j)
5304 call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
5307 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
5308 " mass-scaled friction matrix"
5309 call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
5311 ! Precompute matrices for tinker stochastic integrator
5318 mt1(i,j)=mt1(i,j)+fricvec(k,i)*gsqrm(k,j)
5319 mt2(i,j)=mt2(i,j)+fricvec(k,i)*gsqrp(k,j)
5325 else if (lang.eq.4) then
5326 ! Diagonalize the friction matrix
5331 Ghalf(ind)=fricmat(i,j)
5334 call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
5337 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
5339 call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
5341 ! Determine the number of zero eigenvalues of the friction matrix
5342 nzero=max0(dimen-dimen1,0)
5343 ! do while (fricgam(nzero+1).le.1.0d-5 .and. nzero.lt.dimen)
5346 write (iout,*) "Number of zero eigenvalues:",nzero
5351 fricmat(i,j)=fricmat(i,j) &
5352 +fricvec(i,k)*fricvec(j,k)/fricgam(k)
5357 write (iout,'(//a)') "Generalized inverse of fricmat"
5358 call matout(dimen,dimen,nres6,nres6,fricmat)
5363 if (nfgtasks.gt.1) then
5364 if (fg_rank.eq.0) then
5365 ! The matching BROADCAST for fg processors is called in ERGASTULUM
5371 call MPI_Bcast(10,1,MPI_INTEGER,king,FG_COMM,IERROR)
5373 time_Bcast=time_Bcast+MPI_Wtime()-time00
5375 time_Bcast=time_Bcast+tcpu()-time00
5377 ! print *,"Processor",myrank,
5378 ! & " BROADCAST iorder in SETUP_FRICMAT"
5381 write (iout,*) "setup_fricmat licznik"!,licznik !sp
5387 ! Scatter the friction matrix
5388 call MPI_Scatterv(fricmat(1,1),nginv_counts(0),&
5389 nginv_start(0),MPI_DOUBLE_PRECISION,fcopy(1,1),&
5390 myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
5393 time_scatter=time_scatter+MPI_Wtime()-time00
5394 time_scatter_fmat=time_scatter_fmat+MPI_Wtime()-time00
5396 time_scatter=time_scatter+tcpu()-time00
5397 time_scatter_fmat=time_scatter_fmat+tcpu()-time00
5401 do j=1,2*my_ng_count
5402 fricmat(j,i)=fcopy(i,j)
5405 ! write (iout,*) "My chunk of fricmat"
5406 ! call MATOUT2(my_ng_count,dimen,maxres2,maxres2,fcopy)
5410 end subroutine setup_fricmat
5411 !-----------------------------------------------------------------------------
5412 subroutine stochastic_force(stochforcvec)
5415 use random, only:anorm_distr
5416 ! implicit real*8 (a-h,o-z)
5417 ! include 'DIMENSIONS'
5418 use control, only: tcpu
5423 ! include 'COMMON.VAR'
5424 ! include 'COMMON.CHAIN'
5425 ! include 'COMMON.DERIV'
5426 ! include 'COMMON.GEO'
5427 ! include 'COMMON.LOCAL'
5428 ! include 'COMMON.INTERACT'
5429 ! include 'COMMON.MD'
5430 ! include 'COMMON.TIME1'
5432 ! include 'COMMON.LANGEVIN'
5434 ! include 'COMMON.LANGEVIN.lang0'
5436 ! include 'COMMON.IOUNITS'
5438 real(kind=8) :: x,sig,lowb,highb
5439 real(kind=8) :: ff(3),force(3,0:2*nres),zeta2,lowb2
5440 real(kind=8) :: highb2,sig2,forcvec(6*nres),stochforcvec(6*nres)
5441 real(kind=8) :: time00
5442 logical :: lprn = .false.
5443 integer :: i,j,ind,mnum
5447 stochforc(j,i)=0.0d0
5457 ! Compute the stochastic forces acting on bodies. Store in force.
5463 force(j,i)=anorm_distr(x,sig,lowb,highb)
5471 force(j,i+nres)=anorm_distr(x,sig2,lowb2,highb2)
5475 time_fsample=time_fsample+MPI_Wtime()-time00
5477 time_fsample=time_fsample+tcpu()-time00
5479 ! Compute the stochastic forces acting on virtual-bond vectors.
5485 stochforc(j,i)=ff(j)+0.5d0*force(j,i)
5488 ff(j)=ff(j)+force(j,i)
5490 ! if (itype(i+1,1).ne.ntyp1) then
5492 if (itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
5494 stochforc(j,i)=stochforc(j,i)+force(j,i+nres+1)
5495 ff(j)=ff(j)+force(j,i+nres+1)
5500 stochforc(j,0)=ff(j)+force(j,nnt+nres)
5504 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5505 .and.(mnum.ne.5)) then
5506 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
5508 stochforc(j,i+nres)=force(j,i+nres)
5514 stochforcvec(j)=stochforc(j,0)
5519 stochforcvec(ind+j)=stochforc(j,i)
5525 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5526 .and.(mnum.ne.5)) then
5527 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
5529 stochforcvec(ind+j)=stochforc(j,i+nres)
5535 write (iout,*) "stochforcvec"
5537 write(iout,'(i5,e15.5)') i,stochforcvec(i)
5539 write(iout,*) "Stochastic forces backbone"
5541 write(iout,'(i5,3e15.5)') i,(stochforc(j,i),j=1,3)
5543 write(iout,*) "Stochastic forces side chain"
5545 write(iout,'(i5,3e15.5)') &
5546 i,(stochforc(j,i+nres),j=1,3)
5554 write (iout,*) i,ind
5556 forcvec(ind+j)=force(j,i)
5561 write (iout,*) i,ind
5563 forcvec(j+ind)=force(j,i+nres)
5568 write (iout,*) "forcvec"
5572 write (iout,'(2i3,2f10.5)') i,j,force(j,i),&
5579 write (iout,'(2i3,2f10.5)') i,j,force(j,i+nres),&
5588 end subroutine stochastic_force
5589 !-----------------------------------------------------------------------------
5590 subroutine sdarea(gamvec)
5592 ! Scale the friction coefficients according to solvent accessible surface areas
5593 ! Code adapted from TINKER
5597 ! implicit real*8 (a-h,o-z)
5598 ! include 'DIMENSIONS'
5599 ! include 'COMMON.CONTROL'
5600 ! include 'COMMON.VAR'
5601 ! include 'COMMON.MD'
5603 ! include 'COMMON.LANGEVIN'
5605 ! include 'COMMON.LANGEVIN.lang0'
5607 ! include 'COMMON.CHAIN'
5608 ! include 'COMMON.DERIV'
5609 ! include 'COMMON.GEO'
5610 ! include 'COMMON.LOCAL'
5611 ! include 'COMMON.INTERACT'
5612 ! include 'COMMON.IOUNITS'
5613 ! include 'COMMON.NAMES'
5614 real(kind=8),dimension(2*nres) :: radius,gamvec !(maxres2)
5615 real(kind=8),parameter :: twosix = 1.122462048309372981d0
5616 logical :: lprn = .false.
5617 real(kind=8) :: probe,area,ratio
5618 integer :: i,j,ind,iti,mnum
5620 ! determine new friction coefficients every few SD steps
5622 ! set the atomic radii to estimates of sigma values
5624 ! print *,"Entered sdarea"
5630 ! Load peptide group radii
5633 radius(i)=pstok(mnum)
5635 ! Load side chain radii
5639 radius(i+nres)=restok(iti,mnum)
5642 ! write (iout,*) "i",i," radius",radius(i)
5645 radius(i) = radius(i) / twosix
5646 if (radius(i) .ne. 0.0d0) radius(i) = radius(i) + probe
5649 ! scale atomic friction coefficients by accessible area
5651 if (lprn) write (iout,*) &
5652 "Original gammas, surface areas, scaling factors, new gammas, ",&
5653 "std's of stochastic forces"
5656 if (radius(i).gt.0.0d0) then
5657 call surfatom (i,area,radius)
5658 ratio = dmax1(area/(4.0d0*pi*radius(i)**2),1.0d-1)
5659 if (lprn) write (iout,'(i5,3f10.5,$)') &
5660 i,gamvec(ind+1),area,ratio
5663 gamvec(ind) = ratio * gamvec(ind)
5666 stdforcp(i)=stdfp(mnum)*dsqrt(gamvec(ind))
5667 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcp(i)
5671 if (radius(i+nres).gt.0.0d0) then
5672 call surfatom (i+nres,area,radius)
5673 ratio = dmax1(area/(4.0d0*pi*radius(i+nres)**2),1.0d-1)
5674 if (lprn) write (iout,'(i5,3f10.5,$)') &
5675 i,gamvec(ind+1),area,ratio
5678 gamvec(ind) = ratio * gamvec(ind)
5681 stdforcsc(i)=stdfsc(itype(i,mnum),mnum)*dsqrt(gamvec(ind))
5682 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcsc(i)
5687 end subroutine sdarea
5688 !-----------------------------------------------------------------------------
5690 !-----------------------------------------------------------------------------
5693 ! ###################################################
5694 ! ## COPYRIGHT (C) 1996 by Jay William Ponder ##
5695 ! ## All Rights Reserved ##
5696 ! ###################################################
5698 ! ################################################################
5700 ! ## subroutine surfatom -- exposed surface area of an atom ##
5702 ! ################################################################
5705 ! "surfatom" performs an analytical computation of the surface
5706 ! area of a specified atom; a simplified version of "surface"
5708 ! literature references:
5710 ! T. J. Richmond, "Solvent Accessible Surface Area and
5711 ! Excluded Volume in Proteins", Journal of Molecular Biology,
5714 ! L. Wesson and D. Eisenberg, "Atomic Solvation Parameters
5715 ! Applied to Molecular Dynamics of Proteins in Solution",
5716 ! Protein Science, 1, 227-235 (1992)
5718 ! variables and parameters:
5720 ! ir number of atom for which area is desired
5721 ! area accessible surface area of the atom
5722 ! radius radii of each of the individual atoms
5725 subroutine surfatom(ir,area,radius)
5727 ! implicit real*8 (a-h,o-z)
5728 ! include 'DIMENSIONS'
5730 ! include 'COMMON.GEO'
5731 ! include 'COMMON.IOUNITS'
5733 integer :: nsup,nstart_sup
5734 ! double precision c,dc,dc_old,d_c_work,xloc,xrot,dc_norm
5735 ! common /chain/ c(3,maxres2+2),dc(3,0:maxres2),dc_old(3,0:maxres2),
5736 ! & xloc(3,maxres),xrot(3,maxres),dc_norm(3,0:maxres2),
5737 ! & dc_work(MAXRES6),nres,nres0
5738 integer,parameter :: maxarc=300
5742 integer :: mi,ni,narc
5743 integer :: key(maxarc)
5744 integer :: intag(maxarc)
5745 integer :: intag1(maxarc)
5746 real(kind=8) :: area,arcsum
5747 real(kind=8) :: arclen,exang
5748 real(kind=8) :: delta,delta2
5749 real(kind=8) :: eps,rmove
5750 real(kind=8) :: xr,yr,zr
5751 real(kind=8) :: rr,rrsq
5752 real(kind=8) :: rplus,rminus
5753 real(kind=8) :: axx,axy,axz
5754 real(kind=8) :: ayx,ayy
5755 real(kind=8) :: azx,azy,azz
5756 real(kind=8) :: uxj,uyj,uzj
5757 real(kind=8) :: tx,ty,tz
5758 real(kind=8) :: txb,tyb,td
5759 real(kind=8) :: tr2,tr,txr,tyr
5760 real(kind=8) :: tk1,tk2
5761 real(kind=8) :: thec,the,t,tb
5762 real(kind=8) :: txk,tyk,tzk
5763 real(kind=8) :: t1,ti,tf,tt
5764 real(kind=8) :: txj,tyj,tzj
5765 real(kind=8) :: ccsq,cc,xysq
5766 real(kind=8) :: bsqk,bk,cosine
5767 real(kind=8) :: dsqj,gi,pix2
5768 real(kind=8) :: therk,dk,gk
5769 real(kind=8) :: risqk,rik
5770 real(kind=8) :: radius(2*nres) !(maxatm) (maxatm=maxres2)
5771 real(kind=8) :: ri(maxarc),risq(maxarc)
5772 real(kind=8) :: ux(maxarc),uy(maxarc),uz(maxarc)
5773 real(kind=8) :: xc(maxarc),yc(maxarc),zc(maxarc)
5774 real(kind=8) :: xc1(maxarc),yc1(maxarc),zc1(maxarc)
5775 real(kind=8) :: dsq(maxarc),bsq(maxarc)
5776 real(kind=8) :: dsq1(maxarc),bsq1(maxarc)
5777 real(kind=8) :: arci(maxarc),arcf(maxarc)
5778 real(kind=8) :: ex(maxarc),lt(maxarc),gr(maxarc)
5779 real(kind=8) :: b(maxarc),b1(maxarc),bg(maxarc)
5780 real(kind=8) :: kent(maxarc),kout(maxarc)
5781 real(kind=8) :: ther(maxarc)
5782 logical :: moved,top
5783 logical :: omit(maxarc)
5786 ! maxatm = 2*nres !maxres2 maxres2=2*maxres
5787 ! maxlight = 8*maxatm
5790 ! maxtors = 4*maxatm
5792 ! zero out the surface area for the sphere of interest
5795 ! write (2,*) "ir",ir," radius",radius(ir)
5796 if (radius(ir) .eq. 0.0d0) return
5798 ! set the overlap significance and connectivity shift
5802 delta2 = delta * delta
5807 ! store coordinates and radius of the sphere of interest
5815 ! initialize values of some counters and summations
5824 ! test each sphere to see if it overlaps the sphere of interest
5827 if (i.eq.ir .or. radius(i).eq.0.0d0) goto 30
5828 rplus = rr + radius(i)
5830 if (abs(tx) .ge. rplus) goto 30
5832 if (abs(ty) .ge. rplus) goto 30
5834 if (abs(tz) .ge. rplus) goto 30
5836 ! check for sphere overlap by testing distance against radii
5838 xysq = tx*tx + ty*ty
5839 if (xysq .lt. delta2) then
5846 if (rplus-cc .le. delta) goto 30
5847 rminus = rr - radius(i)
5849 ! check to see if sphere of interest is completely buried
5851 if (cc-abs(rminus) .le. delta) then
5852 if (rminus .le. 0.0d0) goto 170
5856 ! check for too many overlaps with sphere of interest
5858 if (io .ge. maxarc) then
5860 20 format (/,' SURFATOM -- Increase the Value of MAXARC')
5864 ! get overlap between current sphere and sphere of interest
5873 gr(io) = (ccsq+rplus*rminus) / (2.0d0*rr*b1(io))
5879 ! case where no other spheres overlap the sphere of interest
5882 area = 4.0d0 * pi * rrsq
5886 ! case where only one sphere overlaps the sphere of interest
5889 area = pix2 * (1.0d0 + gr(1))
5890 area = mod(area,4.0d0*pi) * rrsq
5894 ! case where many spheres intersect the sphere of interest;
5895 ! sort the intersecting spheres by their degree of overlap
5897 call sort2 (io,gr,key)
5900 intag(i) = intag1(k)
5909 ! get radius of each overlap circle on surface of the sphere
5914 risq(i) = rrsq - gi*gi
5915 ri(i) = sqrt(risq(i))
5916 ther(i) = 0.5d0*pi - asin(min(1.0d0,max(-1.0d0,gr(i))))
5919 ! find boundary of inaccessible area on sphere of interest
5922 if (.not. omit(k)) then
5929 ! check to see if J circle is intersecting K circle;
5930 ! get distance between circle centers and sum of radii
5933 if (omit(j)) goto 60
5934 cc = (txk*xc(j)+tyk*yc(j)+tzk*zc(j))/(bk*b(j))
5935 cc = acos(min(1.0d0,max(-1.0d0,cc)))
5936 td = therk + ther(j)
5938 ! check to see if circles enclose separate regions
5940 if (cc .ge. td) goto 60
5942 ! check for circle J completely inside circle K
5944 if (cc+ther(j) .lt. therk) goto 40
5946 ! check for circles that are essentially parallel
5948 if (cc .gt. delta) goto 50
5953 ! check to see if sphere of interest is completely buried
5956 if (pix2-cc .le. td) goto 170
5962 ! find T value of circle intersections
5965 if (omit(k)) goto 110
5980 ! rotation matrix elements
5992 if (.not. omit(j)) then
5997 ! rotate spheres so K vector colinear with z-axis
5999 uxj = txj*axx + tyj*axy - tzj*axz
6000 uyj = tyj*ayy - txj*ayx
6001 uzj = txj*azx + tyj*azy + tzj*azz
6002 cosine = min(1.0d0,max(-1.0d0,uzj/b(j)))
6003 if (acos(cosine) .lt. therk+ther(j)) then
6004 dsqj = uxj*uxj + uyj*uyj
6009 tr2 = risqk*dsqj - tb*tb
6015 ! get T values of intersection for K circle
6018 tb = min(1.0d0,max(-1.0d0,tb))
6020 if (tyb-txr .lt. 0.0d0) tk1 = pix2 - tk1
6022 tb = min(1.0d0,max(-1.0d0,tb))
6024 if (tyb+txr .lt. 0.0d0) tk2 = pix2 - tk2
6025 thec = (rrsq*uzj-gk*bg(j)) / (rik*ri(j)*b(j))
6026 if (abs(thec) .lt. 1.0d0) then
6028 else if (thec .ge. 1.0d0) then
6030 else if (thec .le. -1.0d0) then
6034 ! see if "tk1" is entry or exit point; check t=0 point;
6035 ! "ti" is exit point, "tf" is entry point
6037 cosine = min(1.0d0,max(-1.0d0, &
6038 (uzj*gk-uxj*rik)/(b(j)*rr)))
6039 if ((acos(cosine)-ther(j))*(tk2-tk1) .le. 0.0d0) then
6047 if (narc .ge. maxarc) then
6049 70 format (/,' SURFATOM -- Increase the Value',&
6053 if (tf .le. ti) then
6074 ! special case; K circle without intersections
6076 if (narc .le. 0) goto 90
6078 ! general case; sum up arclength and set connectivity code
6080 call sort2 (narc,arci,key)
6085 if (narc .gt. 1) then
6088 if (t .lt. arci(j)) then
6089 arcsum = arcsum + arci(j) - t
6090 exang = exang + ex(ni)
6092 if (jb .ge. maxarc) then
6094 80 format (/,' SURFATOM -- Increase the Value',&
6099 kent(jb) = maxarc*i + k
6101 kout(jb) = maxarc*k + i
6110 arcsum = arcsum + pix2 - t
6112 exang = exang + ex(ni)
6115 kent(jb) = maxarc*i + k
6117 kout(jb) = maxarc*k + i
6124 arclen = arclen + gr(k)*arcsum
6127 if (arclen .eq. 0.0d0) goto 170
6128 if (jb .eq. 0) goto 150
6130 ! find number of independent boundaries and check connectivity
6134 if (kout(k) .ne. 0) then
6141 if (m .eq. kent(ii)) then
6144 if (j .eq. jb) goto 150
6156 ! attempt to fix connectivity error by moving atom slightly
6160 140 format (/,' SURFATOM -- Connectivity Error at Atom',i6)
6169 ! compute the exposed surface area for the sphere of interest
6172 area = ib*pix2 + exang + arclen
6173 area = mod(area,4.0d0*pi) * rrsq
6175 ! attempt to fix negative area by moving atom slightly
6177 if (area .lt. 0.0d0) then
6180 160 format (/,' SURFATOM -- Negative Area at Atom',i6)
6191 end subroutine surfatom
6192 !----------------------------------------------------------------
6193 !----------------------------------------------------------------
6194 subroutine alloc_MD_arrays
6195 !EL Allocation of arrays used by MD module
6197 integer :: nres2,nres6
6200 !----------------------
6204 allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
6205 allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
6206 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
6207 allocate(fricvec(nres2,nres2))
6208 allocate(pfric_mat(nres2,nres2),vfric_mat(nres2,nres2))
6209 allocate(afric_mat(nres2,nres2),prand_mat(nres2,nres2))
6210 allocate(vrand_mat1(nres2,nres2),vrand_mat2(nres2,nres2)) !(MAXRES2,MAXRES2)
6211 allocate(pfric0_mat(nres2,nres2,0:maxflag_stoch))
6212 allocate(afric0_mat(nres2,nres2,0:maxflag_stoch))
6213 allocate(vfric0_mat(nres2,nres2,0:maxflag_stoch))
6214 allocate(prand0_mat(nres2,nres2,0:maxflag_stoch))
6215 allocate(vrand0_mat1(nres2,nres2,0:maxflag_stoch))
6216 allocate(vrand0_mat2(nres2,nres2,0:maxflag_stoch)) !(MAXRES2,MAXRES2,0:maxflag_stoch)
6217 allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
6219 allocate(mt1(nres2,nres2),mt2(nres2,nres2),mt3(nres2,nres2)) !(maxres2,maxres2)
6220 !----------------------
6222 ! commom.langevin.lang0
6224 allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
6225 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
6226 allocate(fricvec(nres2,nres2)) !(MAXRES2,MAXRES2)
6227 allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
6228 allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
6231 if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2))
6232 !----------------------
6233 ! commom.hairpin in CSA module
6234 !----------------------
6235 ! common.mce in MCM_MD module
6236 !----------------------
6238 ! common /mdgrad/ in module.energy
6239 ! common /back_constr/ in module.energy
6240 ! common /qmeas/ in module.energy
6243 allocate(potEcomp(0:n_ene+4)) !(0:n_ene+4)
6245 allocate(d_t(3,0:nres2),d_a(3,0:nres2),d_t_old(3,0:nres2)) !(3,0:MAXRES2)
6246 allocate(d_a_work(nres6)) !(6*MAXRES)
6248 allocate(DM(nres2),DU1(nres2),DU2(nres2))
6249 allocate(DMorig(nres2),DU1orig(nres2),DU2orig(nres2))
6251 write (iout,*) "Before A Allocation",nfgtasks-1
6253 allocate(Gmat(nres2,nres2),A(nres2,nres2))
6254 if(.not.allocated(Ginv)) allocate(Ginv(nres2,nres2)) !in control: ergastulum
6255 allocate(Gsqrp(nres2,nres2),Gsqrm(nres2,nres2),Gvec(nres2,nres2)) !(maxres2,maxres2)
6257 allocate(Geigen(nres2)) !(maxres2)
6258 if(.not.allocated(vtot)) allocate(vtot(nres2)) !(maxres2)
6259 ! common /inertia/ in io_conf: parmread
6260 ! real(kind=8),dimension(:),allocatable :: ISC,msc !(ntyp+1)
6261 ! common /langevin/in io read_MDpar
6262 ! real(kind=8),dimension(:),allocatable :: gamsc !(ntyp1)
6263 ! real(kind=8),dimension(:),allocatable :: stdfsc !(ntyp)
6264 ! in io_conf: parmread
6265 ! real(kind=8),dimension(:),allocatable :: restok !(ntyp+1)
6266 ! common /mdpmpi/ in control: ergastulum
6267 if(.not.allocated(ng_start)) allocate(ng_start(0:nfgtasks-1))
6268 if(.not.allocated(ng_counts)) allocate(ng_counts(0:nfgtasks-1))
6269 if(.not.allocated(nginv_counts)) allocate(nginv_counts(0:nfgtasks-1)) !(0:MaxProcs-1)
6270 if(.not.allocated(nginv_start)) allocate(nginv_start(0:nfgtasks)) !(0:MaxProcs)
6271 !----------------------
6272 ! common.muca in read_muca
6273 ! common /double_muca/
6274 ! real(kind=8) :: elow,ehigh,factor,hbin,factor_min
6275 ! real(kind=8),dimension(:),allocatable :: emuca,nemuca,&
6276 ! nemuca2,hist !(4*maxres)
6277 ! common /integer_muca/
6278 ! integer :: nmuca,imtime,muca_smooth
6280 ! real(kind=8),dimension(:),allocatable :: elowi,ehighi !(maxprocs)
6281 !----------------------
6283 ! common /mdgrad/ in module.energy
6284 ! common /back_constr/ in module.energy
6285 ! common /qmeas/ in module.energy
6289 allocate(d_t_work(nres6),d_t_work_new(nres6),d_af_work(nres6))
6290 allocate(d_as_work(nres6),kinetic_force(nres6)) !(MAXRES6)
6291 allocate(d_t_new(3,0:nres2),d_a_old(3,0:nres2),d_a_short(3,0:nres2)) !,d_a !(3,0:MAXRES2)
6292 allocate(stdforcp(nres),stdforcsc(nres)) !(MAXRES)
6293 !----------------------
6295 allocate(D_ban(nres6)) !(MAXRES6) maxres6=6*maxres
6296 ! common /stochcalc/ stochforcvec
6297 allocate(stochforcvec(nres6)) !(MAXRES6) maxres6=6*maxres
6298 !----------------------
6300 end subroutine alloc_MD_arrays
6301 !-----------------------------------------------------------------------------
6302 !-----------------------------------------------------------------------------