2 !-----------------------------------------------------------------------------
15 !-----------------------------------------------------------------------------
17 ! common /mdgrad/ in module.energy
18 ! common /back_constr/ in module.energy
19 ! common /qmeas/ in module.energy
23 real(kind=8),dimension(:),allocatable :: d_t_work,&
24 d_t_work_new,d_af_work,d_as_work,kinetic_force !(MAXRES6)
25 real(kind=8),dimension(:,:),allocatable :: d_t_new,&
26 d_a_old,d_a_short!,d_a !(3,0:MAXRES2)
27 ! real(kind=8),dimension(:),allocatable :: d_a_work !(6*MAXRES)
28 ! real(kind=8),dimension(:,:),allocatable :: Gmat,Ginv,A,&
29 ! Gsqrp,Gsqrm,Gvec !(maxres2,maxres2)
30 ! real(kind=8),dimension(:),allocatable :: Geigen !(maxres2)
31 ! integer :: dimen,dimen1,dimen3
32 ! integer :: lang,count_reset_moment,count_reset_vel
33 ! logical :: reset_moment,reset_vel,rattle,RESPA
36 ! real(kind=8) :: rwat,etawat,stdfp,cPoise
37 ! real(kind=8),dimension(:),allocatable :: gamsc !(ntyp1)
38 ! real(kind=8),dimension(:),allocatable :: stdfsc !(ntyp)
39 real(kind=8),dimension(:),allocatable :: stdforcp,stdforcsc !(MAXRES)
40 !-----------------------------------------------------------------------------
44 ! ###################################################
45 ! ## COPYRIGHT (C) 1992 by Jay William Ponder ##
46 ! ## All Rights Reserved ##
47 ! ###################################################
49 ! #############################################################
51 ! ## sizes.i -- parameter values to set array dimensions ##
53 ! #############################################################
56 ! "sizes.i" sets values for critical array dimensions used
57 ! throughout the software; these parameters will fix the size
58 ! of the largest systems that can be handled; values too large
59 ! for the computer's memory and/or swap space to accomodate
60 ! will result in poor performance or outright failure
62 ! parameter: maximum allowed number of:
64 ! maxatm atoms in the molecular system
65 ! maxval atoms directly bonded to an atom
66 ! maxgrp ! user-defined groups of atoms
67 ! maxtyp force field atom type definitions
68 ! maxclass force field atom class definitions
69 ! maxkey lines in the keyword file
70 ! maxrot bonds for torsional rotation
71 ! maxvar optimization variables (vector storage)
72 ! maxopt optimization variables (matrix storage)
73 ! maxhess off-diagonal Hessian elements
74 ! maxlight sites for method of lights neighbors
75 ! maxvib vibrational frequencies
76 ! maxgeo distance geometry points
77 ! maxcell unit cells in replicated crystal
78 ! maxring 3-, 4-, or 5-membered rings
79 ! maxfix geometric restraints
80 ! maxbio biopolymer atom definitions
81 ! maxres residues in the macromolecule
82 ! maxamino amino acid residue types
83 ! maxnuc nucleic acid residue types
84 ! maxbnd covalent bonds in molecular system
85 ! maxang bond angles in molecular system
86 ! maxtors torsional angles in molecular system
87 ! maxpi atoms in conjugated pisystem
88 ! maxpib covalent bonds involving pisystem
89 ! maxpit torsional angles involving pisystem
92 !el integer maxatm,maxval,maxgrp
93 !el integer maxtyp,maxclass,maxkey
94 !el integer maxrot,maxopt
95 !el integer maxhess,maxlight,maxvib
96 !el integer maxgeo,maxcell,maxring
97 !el integer maxfix,maxbio
98 !el integer maxamino,maxnuc,maxbnd
99 !el integer maxang,maxtors,maxpi
100 !el integer maxpib,maxpit
101 ! integer :: maxatm !=2*nres !maxres2 maxres2=2*maxres
102 ! integer,parameter :: maxval=8
103 ! integer,parameter :: maxgrp=1000
104 ! integer,parameter :: maxtyp=3000
105 ! integer,parameter :: maxclass=500
106 ! integer,parameter :: maxkey=10000
107 ! integer,parameter :: maxrot=1000
108 ! integer,parameter :: maxopt=1000
109 ! integer,parameter :: maxhess=1000000
110 ! integer :: maxlight !=8*maxatm
111 ! integer,parameter :: maxvib=1000
112 ! integer,parameter :: maxgeo=1000
113 ! integer,parameter :: maxcell=10000
114 ! integer,parameter :: maxring=10000
115 ! integer,parameter :: maxfix=10000
116 ! integer,parameter :: maxbio=10000
117 ! integer,parameter :: maxamino=31
118 ! integer,parameter :: maxnuc=12
119 ! integer :: maxbnd !=2*maxatm
120 ! integer :: maxang !=3*maxatm
121 ! integer :: maxtors !=4*maxatm
122 ! integer,parameter :: maxpi=100
123 ! integer,parameter :: maxpib=2*maxpi
124 ! integer,parameter :: maxpit=4*maxpi
125 !-----------------------------------------------------------------------------
126 ! Maximum number of seed
127 ! integer,parameter :: max_seed=1
128 !-----------------------------------------------------------------------------
129 real(kind=8),dimension(:),allocatable :: stochforcvec !(MAXRES6) maxres6=6*maxres
130 ! common /stochcalc/ stochforcvec
131 !-----------------------------------------------------------------------------
132 ! common /przechowalnia/ subroutines: rattle1,rattle2,rattle_brown
133 real(kind=8),dimension(:,:),allocatable :: GGinv !(2*nres,2*nres) maxres2=2*maxres
134 real(kind=8),dimension(:,:,:),allocatable :: gdc !(3,2*nres,2*nres) maxres2=2*maxres
135 real(kind=8),dimension(:,:),allocatable :: Cmat !(2*nres,2*nres) maxres2=2*maxres
136 !-----------------------------------------------------------------------------
137 ! common /syfek/ subroutines: friction_force,setup_fricmat
138 !el real(kind=8),dimension(:),allocatable :: gamvec !(MAXRES6) or (MAXRES2)
139 !-----------------------------------------------------------------------------
140 ! common /przechowalnia/ subroutines: friction_force,setup_fricmat
141 real(kind=8),dimension(:,:),allocatable :: ginvfric !(2*nres,2*nres) !maxres2=2*maxres
142 !-----------------------------------------------------------------------------
143 ! common /przechowalnia/ subroutine: setup_fricmat
144 real(kind=8),dimension(:,:),allocatable :: fcopy !(2*nres,2*nres)
145 !-----------------------------------------------------------------------------
148 !-----------------------------------------------------------------------------
150 !-----------------------------------------------------------------------------
152 !-----------------------------------------------------------------------------
153 subroutine brown_step(itime)
154 !------------------------------------------------
155 ! Perform a single Euler integration step of Brownian dynamics
156 !------------------------------------------------
157 ! implicit real*8 (a-h,o-z)
159 use control, only: tcpu
162 ! use io_conf, only:cartprint
163 ! include 'DIMENSIONS'
167 ! include 'COMMON.CONTROL'
168 ! include 'COMMON.VAR'
169 ! include 'COMMON.MD'
171 ! include 'COMMON.LANGEVIN'
173 ! include 'COMMON.LANGEVIN.lang0'
175 ! include 'COMMON.CHAIN'
176 ! include 'COMMON.DERIV'
177 ! include 'COMMON.GEO'
178 ! include 'COMMON.LOCAL'
179 ! include 'COMMON.INTERACT'
180 ! include 'COMMON.IOUNITS'
181 ! include 'COMMON.NAMES'
182 ! include 'COMMON.TIME1'
183 real(kind=8),dimension(6*nres) :: zapas !(MAXRES6) maxres6=6*maxres
184 integer :: rstcount !ilen,
186 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
187 real(kind=8),dimension(6*nres,2*nres) :: Bmat,GBmat,Tmat !(MAXRES6,MAXRES2) (maxres2=2*maxres,maxres6=6*maxres)
188 real(kind=8),dimension(2*nres,2*nres) :: Cmat_,Cinv !(maxres2,maxres2) maxres2=2*maxres
189 real(kind=8),dimension(6*nres,6*nres) :: Pmat !(maxres6,maxres6) maxres6=6*maxres
190 ! real(kind=8),dimension(:,:),allocatable :: Bmat,GBmat,Tmat !(MAXRES6,MAXRES2) (maxres2=2*maxres,maxres6=6*maxres)
191 ! real(kind=8),dimension(:,:),allocatable :: Cmat_,Cinv !(maxres2,maxres2) maxres2=2*maxres
192 ! real(kind=8),dimension(:,:),allocatable :: Pmat !(maxres6,maxres6) maxres6=6*maxres
193 real(kind=8),dimension(6*nres) :: Td !(maxres6) maxres6=6*maxres
194 real(kind=8),dimension(2*nres) :: ppvec !(maxres2) maxres2=2*maxres
195 !el common /stochcalc/ stochforcvec
196 !el real(kind=8),dimension(3) :: cm !el
197 !el common /gucio/ cm
199 logical :: lprn = .false.,lprn1 = .false.
200 integer :: maxiter = 5
201 real(kind=8) :: difftol = 1.0d-5
202 real(kind=8) :: xx,diffmax,blen2,diffbond,tt0
203 integer :: i,j,nbond,k,ind,ind1,iter
204 integer :: nres2,nres6
208 ! if (.not.allocated(Bmat)) allocate(Bmat(nres6,nres2))
209 ! if (.not.allocated(GBmat)) allocate (GBmat(nres6,nres2))
210 ! if (.not.allocated(Tmat)) allocate (Tmat(nres6,nres2))
211 ! if (.not.allocated(Cmat_)) allocate(Cmat_(nres2,nres2))
212 ! if (.not.allocated(Cinv)) allocate (Cinv(nres2,nres2))
213 ! if (.not.allocated(Pmat)) allocate(Pmat(6*nres,6*nres))
215 if (.not.allocated(stochforcvec)) allocate(stochforcvec(nres6)) !(MAXRES6) maxres6=6*maxres
219 if (itype(i,1).ne.10) nbond=nbond+1
223 write (iout,*) "Generalized inverse of fricmat"
224 call matout(dimen,dimen,nres6,nres6,fricmat)
236 Bmat(ind+j,ind1)=dC_norm(j,i)
241 if (itype(i,1).ne.10) then
244 Bmat(ind+j,ind1)=dC_norm(j,i+nres)
250 write (iout,*) "Matrix Bmat"
251 call MATOUT(nbond,dimen,nres6,nres6,Bmat)
257 GBmat(i,j)=GBmat(i,j)+fricmat(i,k)*Bmat(k,j)
262 write (iout,*) "Matrix GBmat"
263 call MATOUT(nbond,dimen,nres6,nres2,Gbmat)
269 Cmat_(i,j)=Cmat_(i,j)+Bmat(k,i)*GBmat(k,j)
274 write (iout,*) "Matrix Cmat"
275 call MATOUT(nbond,nbond,nres2,nres2,Cmat_)
277 call matinvert(nbond,nres2,Cmat_,Cinv,osob)
279 write (iout,*) "Matrix Cinv"
280 call MATOUT(nbond,nbond,nres2,nres2,Cinv)
286 Tmat(i,j)=Tmat(i,j)+GBmat(i,k)*Cinv(k,j)
291 write (iout,*) "Matrix Tmat"
292 call MATOUT(nbond,dimen,nres6,nres2,Tmat)
302 Pmat(i,j)=Pmat(i,j)-Tmat(i,k)*Bmat(j,k)
307 write (iout,*) "Matrix Pmat"
308 call MATOUT(dimen,dimen,nres6,nres6,Pmat)
315 Td(i)=Td(i)+vbl*Tmat(i,ind)
318 if (itype(k,1).ne.10) then
320 Td(i)=Td(i)+vbldsc0(1,itype(k,1))*Tmat(i,ind)
325 write (iout,*) "Vector Td"
327 write (iout,'(i5,f10.5)') i,Td(i)
330 call stochastic_force(stochforcvec)
332 write (iout,*) "stochforcvec"
334 write (iout,*) i,stochforcvec(i)
338 zapas(j)=-gcart(j,0)+stochforcvec(j)
340 dC_work(j)=dC_old(j,0)
346 zapas(ind)=-gcart(j,i)+stochforcvec(ind)
347 dC_work(ind)=dC_old(j,i)
351 if (itype(i,1).ne.10) then
354 zapas(ind)=-gxcart(j,i)+stochforcvec(ind)
355 dC_work(ind)=dC_old(j,i+nres)
361 write (iout,*) "Initial d_t_work"
363 write (iout,*) i,d_t_work(i)
370 d_t_work(i)=d_t_work(i)+fricmat(i,j)*zapas(j)
377 zapas(i)=zapas(i)+Pmat(i,j)*(dC_work(j)+d_t_work(j)*d_time)
381 write (iout,*) "Final d_t_work and zapas"
383 write (iout,*) i,d_t_work(i),zapas(i)
397 dc_work(ind+j)=dc(j,i)
403 d_t(j,i+nres)=d_t_work(ind+j)
404 dc(j,i+nres)=zapas(ind+j)
405 dc_work(ind+j)=dc(j,i+nres)
411 write (iout,*) "Before correction for rotational lengthening"
412 write (iout,*) "New coordinates",&
413 " and differences between actual and standard bond lengths"
418 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
422 if (itype(i,1).ne.10) then
424 xx=vbld(i+nres)-vbldsc0(1,itype(i,1))
425 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
426 i,(dC(j,i+nres),j=1,3),xx
430 ! Second correction (rotational lengthening)
436 blen2 = scalar(dc(1,i),dc(1,i))
437 ppvec(ind)=2*vbl**2-blen2
438 diffbond=dabs(vbl-dsqrt(blen2))
439 if (diffbond.gt.diffmax) diffmax=diffbond
440 if (ppvec(ind).gt.0.0d0) then
441 ppvec(ind)=dsqrt(ppvec(ind))
446 write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind)
450 if (itype(i,1).ne.10) then
452 blen2 = scalar(dc(1,i+nres),dc(1,i+nres))
453 ppvec(ind)=2*vbldsc0(1,itype(i,1))**2-blen2
454 diffbond=dabs(vbldsc0(1,itype(i,1))-dsqrt(blen2))
455 if (diffbond.gt.diffmax) diffmax=diffbond
456 if (ppvec(ind).gt.0.0d0) then
457 ppvec(ind)=dsqrt(ppvec(ind))
462 write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind)
466 if (lprn) write (iout,*) "iter",iter," diffmax",diffmax
467 if (diffmax.lt.difftol) goto 10
471 Td(i)=Td(i)+ppvec(j)*Tmat(i,j)
477 zapas(i)=zapas(i)+Pmat(i,j)*dc_work(j)
488 dc_work(ind+j)=zapas(ind+j)
493 if (itype(i,1).ne.10) then
495 dc(j,i+nres)=zapas(ind+j)
496 dc_work(ind+j)=zapas(ind+j)
501 ! Building the chain from the newly calculated coordinates
504 if (large.and. mod(itime,ntwe).eq.0) then
505 write (iout,*) "Cartesian and internal coordinates: step 1"
508 write (iout,'(a)') "Potential forces"
510 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(-gcart(j,i),j=1,3),&
513 write (iout,'(a)') "Stochastic forces"
515 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(stochforc(j,i),j=1,3),&
516 (stochforc(j,i+nres),j=1,3)
518 write (iout,'(a)') "Velocities"
520 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
521 (d_t(j,i+nres),j=1,3)
526 write (iout,*) "After correction for rotational lengthening"
527 write (iout,*) "New coordinates",&
528 " and differences between actual and standard bond lengths"
533 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
537 if (itype(i,1).ne.10) then
539 xx=vbld(i+nres)-vbldsc0(1,itype(i,1))
540 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
541 i,(dC(j,i+nres),j=1,3),xx
546 ! write (iout,*) "Too many attempts at correcting the bonds"
554 ! Calculate energy and forces
556 call etotal(potEcomp)
557 potE=potEcomp(0)-potEcomp(20)
561 ! Calculate the kinetic and total energy and the kinetic temperature
564 t_enegrad=t_enegrad+MPI_Wtime()-tt0
566 t_enegrad=t_enegrad+tcpu()-tt0
569 kinetic_T=2.0d0/(dimen*Rb)*EK
571 end subroutine brown_step
572 !-----------------------------------------------------------------------------
574 !-----------------------------------------------------------------------------
575 subroutine gauss(RO,AP,MT,M,N,*)
577 ! CALCULATES (RO**(-1))*AP BY GAUSS ELIMINATION
578 ! RO IS A SQUARE MATRIX
579 ! THE CALCULATED PRODUCT IS STORED IN AP
580 ! ABNORMAL EXIT IF RO IS SINGULAR
582 integer :: MT, M, N, M1,I,J,IM,&
584 real(kind=8) :: RO(MT,M),AP(MT,N),X,RM,PR,Y
590 if(dabs(X).le.1.0D-13) return 1
602 if(DABS(RO(J,I)).LE.RM) goto 2
616 if(dabs(X).le.1.0E-13) return 1
625 8 AP(J,K)=AP(J,K)-Y*AP(I,K)
627 9 RO(J,K)=RO(J,K)-Y*RO(I,K)
631 if(dabs(X).le.1.0E-13) return 1
641 15 X=X-AP(K,J)*RO(MI,K)
646 !-----------------------------------------------------------------------------
648 !-----------------------------------------------------------------------------
649 subroutine kinetic(KE_total)
650 !----------------------------------------------------------------
651 ! This subroutine calculates the total kinetic energy of the chain
652 !-----------------------------------------------------------------
654 ! implicit real*8 (a-h,o-z)
655 ! include 'DIMENSIONS'
656 ! include 'COMMON.VAR'
657 ! include 'COMMON.CHAIN'
658 ! include 'COMMON.DERIV'
659 ! include 'COMMON.GEO'
660 ! include 'COMMON.LOCAL'
661 ! include 'COMMON.INTERACT'
662 ! include 'COMMON.MD'
663 ! include 'COMMON.IOUNITS'
664 real(kind=8) :: KE_total,mscab
666 integer :: i,j,k,iti,mnum
667 real(kind=8) :: KEt_p,KEt_sc,KEr_p,KEr_sc,incr(3),&
670 write (iout,*) "Velocities, kietic"
672 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
673 (d_t(j,i+nres),j=1,3)
678 ! write (iout,*) "ISC",(isc(itype(i,1)),i=1,nres)
679 ! The translational part for peptide virtual bonds
685 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
686 ! write (iout,*) "Kinetic trp:",i,(incr(j),j=1,3)
688 v(j)=incr(j)+0.5d0*d_t(j,i)
690 vtot(i)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
691 KEt_p=KEt_p+mp(mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
693 incr(j)=incr(j)+d_t(j,i)
696 ! write(iout,*) 'KEt_p', KEt_p
697 ! The translational part for the side chain virtual bond
698 ! Only now we can initialize incr with zeros. It must be equal
699 ! to the velocities of the first Calpha.
705 iti=iabs(itype(i,mnum))
711 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
712 if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
713 .or.(mnum.eq.5)) then
719 v(j)=incr(j)+d_t(j,nres+i)
722 ! write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3)
723 ! write (iout,*) "i",i," msc",msc(iti)," v",(v(j),j=1,3)
724 KEt_sc=KEt_sc+mscab*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
725 vtot(i+nres)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
727 incr(j)=incr(j)+d_t(j,i)
731 ! write(iout,*) 'KEt_sc', KEt_sc
732 ! The part due to stretching and rotation of the peptide groups
736 ! write (iout,*) "i",i
737 ! write (iout,*) "i",i," mag1",mag1," mag2",mag2
741 ! write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3)
742 KEr_p=KEr_p+Ip(mnum)*(incr(1)*incr(1)+incr(2)*incr(2) &
746 ! write(iout,*) 'KEr_p', KEr_p
747 ! The rotational part of the side chain virtual bond
751 iti=iabs(itype(i,mnum))
752 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
753 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
754 .and.(mnum.ne.5)) then
756 incr(j)=d_t(j,nres+i)
758 ! write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3)
759 KEr_sc=KEr_sc+Isc(iti,mnum)*(incr(1)*incr(1)+incr(2)*incr(2)+ &
763 ! The total kinetic energy
765 ! write(iout,*) 'KEr_sc', KEr_sc
766 KE_total=0.5d0*(KEt_p+KEt_sc+0.25d0*KEr_p+KEr_sc)
767 ! write (iout,*) "KE_total",KE_total
769 end subroutine kinetic
770 !-----------------------------------------------------------------------------
772 !-----------------------------------------------------------------------------
774 !------------------------------------------------
775 ! The driver for molecular dynamics subroutines
776 !------------------------------------------------
779 use control, only:tcpu,ovrtim
780 ! use io_comm, only:ilen
782 use compare, only:secondary2,hairpin
783 use io, only:cartout,statout
784 ! implicit real*8 (a-h,o-z)
785 ! include 'DIMENSIONS'
788 integer :: IERROR,ERRCODE
790 ! include 'COMMON.SETUP'
791 ! include 'COMMON.CONTROL'
792 ! include 'COMMON.VAR'
793 ! include 'COMMON.MD'
795 ! include 'COMMON.LANGEVIN'
797 ! include 'COMMON.LANGEVIN.lang0'
799 ! include 'COMMON.CHAIN'
800 ! include 'COMMON.DERIV'
801 ! include 'COMMON.GEO'
802 ! include 'COMMON.LOCAL'
803 ! include 'COMMON.INTERACT'
804 ! include 'COMMON.IOUNITS'
805 ! include 'COMMON.NAMES'
806 ! include 'COMMON.TIME1'
807 ! include 'COMMON.HAIRPIN'
808 real(kind=8),dimension(3) :: L,vcm
810 real(kind=8),dimension(6*nres) :: v_work,v_transf !(maxres6) maxres6=6*maxres
812 integer :: rstcount !ilen,
814 character(len=50) :: tytul
815 !el common /gucio/ cm
817 integer,dimension(4,nres) :: iharp !(4,nres/3)(4,maxres/3)
819 real(kind=8) :: tt0,scalfac
820 integer :: nres2,itime
825 print *,"MY tmpdir",tmpdir,ilen(tmpdir)
826 if (ilen(tmpdir).gt.0) &
827 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_" &
828 //liczba(:ilen(liczba))//'.rst')
830 if (ilen(tmpdir).gt.0) &
831 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"//'.rst')
838 write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
844 print *,"just befor setup matix",nres
845 ! Determine the inverse of the inertia matrix.
846 call setup_MD_matrices
848 print *,"AFTER SETUP MATRICES"
850 print *,"AFTER INIT MD"
853 t_MDsetup = MPI_Wtime()-tt0
855 t_MDsetup = tcpu()-tt0
858 ! Entering the MD loop
864 if (lang.eq.2 .or. lang.eq.3) then
868 call sd_verlet_p_setup
870 call sd_verlet_ciccotti_setup
874 pfric0_mat(i,j,0)=pfric_mat(i,j)
875 afric0_mat(i,j,0)=afric_mat(i,j)
876 vfric0_mat(i,j,0)=vfric_mat(i,j)
877 prand0_mat(i,j,0)=prand_mat(i,j)
878 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
879 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
884 flag_stoch(i)=.false.
888 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
890 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
894 else if (lang.eq.1 .or. lang.eq.4) then
895 print *,"before setup_fricmat"
897 print *,"after setup_fricmat"
900 t_langsetup=MPI_Wtime()-tt0
903 t_langsetup=tcpu()-tt0
906 do itime=1,n_timestep
908 if (large.and. mod(itime,ntwe).eq.0) &
909 write (iout,*) "itime",itime
911 if (lang.gt.0 .and. surfarea .and. &
912 mod(itime,reset_fricmat).eq.0) then
913 if (lang.eq.2 .or. lang.eq.3) then
917 call sd_verlet_p_setup
919 call sd_verlet_ciccotti_setup
923 pfric0_mat(i,j,0)=pfric_mat(i,j)
924 afric0_mat(i,j,0)=afric_mat(i,j)
925 vfric0_mat(i,j,0)=vfric_mat(i,j)
926 prand0_mat(i,j,0)=prand_mat(i,j)
927 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
928 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
933 flag_stoch(i)=.false.
936 else if (lang.eq.1 .or. lang.eq.4) then
939 write (iout,'(a,i10)') &
940 "Friction matrix reset based on surface area, itime",itime
942 if (reset_vel .and. tbf .and. lang.eq.0 &
943 .and. mod(itime,count_reset_vel).eq.0) then
945 write(iout,'(a,f20.2)') &
946 "Velocities reset to random values, time",totT
949 d_t_old(j,i)=d_t(j,i)
953 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
957 d_t(j,0)=d_t(j,0)-vcm(j)
960 kinetic_T=2.0d0/(dimen3*Rb)*EK
961 scalfac=dsqrt(T_bath/kinetic_T)
962 write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
965 d_t_old(j,i)=scalfac*d_t(j,i)
971 ! Time-reversible RESPA algorithm
972 ! (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
973 call RESPA_step(itime)
975 ! Variable time step algorithm.
976 call velverlet_step(itime)
980 call brown_step(itime)
982 print *,"Brown dynamics not here!"
984 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
991 if (mod(itime,ntwe).eq.0) then
994 ! call check_ecartint
1004 v_work(ind)=d_t(j,i)
1009 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1.and.mnum.ne.5) then
1012 v_work(ind)=d_t(j,i+nres)
1017 write (66,'(80f10.5)') &
1018 ((d_t(j,i),j=1,3),i=0,nres-1),((d_t(j,i+nres),j=1,3),i=1,nres)
1022 v_transf(i)=v_transf(i)+gvec(j,i)*v_work(j)
1024 v_transf(i)= v_transf(i)*dsqrt(geigen(i))
1026 write (67,'(80f10.5)') (v_transf(i),i=1,ind)
1029 if (mod(itime,ntwx).eq.0) then
1031 write (tytul,'("time",f8.2)') totT
1033 call hairpin(.true.,nharp,iharp)
1034 call secondary2(.true.)
1035 call pdbout(potE,tytul,ipdb)
1040 if (rstcount.eq.1000.or.itime.eq.n_timestep) then
1041 open(irest2,file=rest2name,status='unknown')
1042 write(irest2,*) totT,EK,potE,totE,t_bath
1044 ! AL 4/17/17: Now writing d_t(0,:) too
1046 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
1048 ! AL 4/17/17: Now writing d_c(0,:) too
1050 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
1058 t_MD=MPI_Wtime()-tt0
1062 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))') &
1064 'MD calculations setup:',t_MDsetup,&
1065 'Energy & gradient evaluation:',t_enegrad,&
1066 'Stochastic MD setup:',t_langsetup,&
1067 'Stochastic MD step setup:',t_sdsetup,&
1069 write (iout,'(/28(1h=),a25,27(1h=))') &
1070 ' End of MD calculation '
1072 write (iout,*) "time for etotal",t_etotal," elong",t_elong,&
1074 write (iout,*) "time_fric",time_fric," time_stoch",time_stoch,&
1075 " time_fricmatmult",time_fricmatmult," time_fsample ",&
1080 !-----------------------------------------------------------------------------
1081 subroutine velverlet_step(itime)
1082 !-------------------------------------------------------------------------------
1083 ! Perform a single velocity Verlet step; the time step can be rescaled if
1084 ! increments in accelerations exceed the threshold
1085 !-------------------------------------------------------------------------------
1086 ! implicit real*8 (a-h,o-z)
1087 ! include 'DIMENSIONS'
1089 use control, only:tcpu
1091 use minimm, only:minim_dc
1094 integer :: ierror,ierrcode
1095 real(kind=8) :: errcode
1097 ! include 'COMMON.SETUP'
1098 ! include 'COMMON.VAR'
1099 ! include 'COMMON.MD'
1101 ! include 'COMMON.LANGEVIN'
1103 ! include 'COMMON.LANGEVIN.lang0'
1105 ! include 'COMMON.CHAIN'
1106 ! include 'COMMON.DERIV'
1107 ! include 'COMMON.GEO'
1108 ! include 'COMMON.LOCAL'
1109 ! include 'COMMON.INTERACT'
1110 ! include 'COMMON.IOUNITS'
1111 ! include 'COMMON.NAMES'
1112 ! include 'COMMON.TIME1'
1113 ! include 'COMMON.MUCA'
1114 real(kind=8),dimension(3) :: vcm,incr
1115 real(kind=8),dimension(3) :: L
1116 integer :: count,rstcount !ilen,
1118 character(len=50) :: tytul
1119 integer :: maxcount_scale = 30
1120 !el common /gucio/ cm
1121 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1122 !el common /stochcalc/ stochforcvec
1123 integer :: icount_scale,itime_scal,i,j,ifac_time,iretcode,itime
1125 real(kind=8) :: epdrift,tt0,fac_time
1127 if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres)) !(MAXRES6) maxres6=6*maxres
1133 else if (lang.eq.2 .or. lang.eq.3) then
1135 call stochastic_force(stochforcvec)
1138 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
1140 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1147 icount_scale=icount_scale+1
1148 ! write(iout,*) "icount_scale",icount_scale,scalel
1149 if (icount_scale.gt.maxcount_scale) then
1151 "ERROR: too many attempts at scaling down the time step. ",&
1152 "amax=",amax,"epdrift=",epdrift,&
1153 "damax=",damax,"edriftmax=",edriftmax,&
1157 call MPI_Abort(MPI_COMM_WORLD,IERROR,IERRCODE)
1161 ! First step of the velocity Verlet algorithm
1166 else if (lang.eq.3) then
1168 call sd_verlet1_ciccotti
1170 else if (lang.eq.1) then
1175 ! Build the chain from the newly calculated coordinates
1176 call chainbuild_cart
1177 if (rattle) call rattle1
1179 if (large.and. mod(itime,ntwe).eq.0) then
1180 write (iout,*) "Cartesian and internal coordinates: step 1"
1185 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),&
1186 (dc(j,i+nres),j=1,3)
1188 write (iout,*) "Accelerations"
1190 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1191 (d_a(j,i+nres),j=1,3)
1193 write (iout,*) "Velocities, step 1"
1195 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1196 (d_t(j,i+nres),j=1,3)
1205 ! Calculate energy and forces
1207 call etotal(potEcomp)
1208 ! AL 4/17/17: Reduce the steps if NaNs occurred.
1209 if (potEcomp(0).gt.0.99e18 .or. isnan(potEcomp(0)).gt.0) then
1210 call enerprint(potEcomp)
1212 if (icount_scale.gt.15) then
1213 write (iout,*) "Tu jest problem",potEcomp(0),d_time
1214 ! call gen_rand_conf(1,*335)
1215 ! call minim_dc(potEcomp(0),iretcode,100)
1218 ! call etotal(potEcomp)
1219 ! write(iout,*) "needed to repara,",potEcomp
1222 ! 335 write(iout,*) "Failed genrand"
1226 if (large.and. mod(itime,ntwe).eq.0) &
1227 call enerprint(potEcomp)
1230 t_etotal=t_etotal+MPI_Wtime()-tt0
1232 t_etotal=t_etotal+tcpu()-tt0
1235 potE=potEcomp(0)-potEcomp(20)
1237 ! Get the new accelerations
1240 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1242 t_enegrad=t_enegrad+tcpu()-tt0
1244 ! Determine maximum acceleration and scale down the timestep if needed
1246 amax=amax/(itime_scal+1)**2
1247 call predict_edrift(epdrift)
1248 ! write(iout,*) "amax=",amax,damax,epdrift,edriftmax,amax/(itime_scal+1)
1250 ! write (iout,*) "before enter if",scalel,icount_scale
1251 if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then
1252 ! write(iout,*) "I enter if"
1253 ! Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
1255 ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax)) &
1257 itime_scal=itime_scal+ifac_time
1258 ! fac_time=dmin1(damax/amax,0.5d0)
1259 fac_time=0.5d0**ifac_time
1260 d_time=d_time*fac_time
1261 if (lang.eq.2 .or. lang.eq.3) then
1263 ! write (iout,*) "Calling sd_verlet_setup: 1"
1264 ! Rescale the stochastic forces and recalculate or restore
1265 ! the matrices of tinker integrator
1266 if (itime_scal.gt.maxflag_stoch) then
1267 if (large) write (iout,'(a,i5,a)') &
1268 "Calculate matrices for stochastic step;",&
1269 " itime_scal ",itime_scal
1271 call sd_verlet_p_setup
1273 call sd_verlet_ciccotti_setup
1275 write (iout,'(2a,i3,a,i3,1h.)') &
1276 "Warning: cannot store matrices for stochastic",&
1277 " integration because the index",itime_scal,&
1278 " is greater than",maxflag_stoch
1279 write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",&
1280 " integration Langevin algorithm for better efficiency."
1281 else if (flag_stoch(itime_scal)) then
1282 if (large) write (iout,'(a,i5,a,l1)') &
1283 "Restore matrices for stochastic step; itime_scal ",&
1284 itime_scal," flag ",flag_stoch(itime_scal)
1287 pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
1288 afric_mat(i,j)=afric0_mat(i,j,itime_scal)
1289 vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
1290 prand_mat(i,j)=prand0_mat(i,j,itime_scal)
1291 vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
1292 vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
1296 if (large) write (iout,'(2a,i5,a,l1)') &
1297 "Calculate & store matrices for stochastic step;",&
1298 " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
1300 call sd_verlet_p_setup
1302 call sd_verlet_ciccotti_setup
1304 flag_stoch(ifac_time)=.true.
1307 pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
1308 afric0_mat(i,j,itime_scal)=afric_mat(i,j)
1309 vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
1310 prand0_mat(i,j,itime_scal)=prand_mat(i,j)
1311 vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
1312 vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
1316 fac_time=1.0d0/dsqrt(fac_time)
1318 stochforcvec(i)=fac_time*stochforcvec(i)
1321 else if (lang.eq.1) then
1322 ! Rescale the accelerations due to stochastic forces
1323 fac_time=1.0d0/dsqrt(fac_time)
1325 d_as_work(i)=d_as_work(i)*fac_time
1328 if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)') &
1329 "itime",itime," Timestep scaled down to ",&
1330 d_time," ifac_time",ifac_time," itime_scal",itime_scal
1332 ! Second step of the velocity Verlet algorithm
1337 else if (lang.eq.3) then
1339 call sd_verlet2_ciccotti
1341 else if (lang.eq.1) then
1346 if (rattle) call rattle2
1349 if (d_time.ne.d_time0) then
1352 if (lang.eq.2 .or. lang.eq.3) then
1353 if (large) write (iout,'(a)') &
1354 "Restore original matrices for stochastic step"
1355 ! write (iout,*) "Calling sd_verlet_setup: 2"
1356 ! Restore the matrices of tinker integrator if the time step has been restored
1359 pfric_mat(i,j)=pfric0_mat(i,j,0)
1360 afric_mat(i,j)=afric0_mat(i,j,0)
1361 vfric_mat(i,j)=vfric0_mat(i,j,0)
1362 prand_mat(i,j)=prand0_mat(i,j,0)
1363 vrand_mat1(i,j)=vrand0_mat1(i,j,0)
1364 vrand_mat2(i,j)=vrand0_mat2(i,j,0)
1372 ! Calculate the kinetic and the total energy and the kinetic temperature
1376 ! call kinetic1(EK1)
1377 ! write (iout,*) "step",itime," EK",EK," EK1",EK1
1379 ! Couple the system to Berendsen bath if needed
1380 if (tbf .and. lang.eq.0) then
1383 kinetic_T=2.0d0/(dimen3*Rb)*EK
1384 ! Backup the coordinates, velocities, and accelerations
1388 d_t_old(j,i)=d_t(j,i)
1389 d_a_old(j,i)=d_a(j,i)
1393 if (mod(itime,ntwe).eq.0 .and. large) then
1394 write (iout,*) "Velocities, step 2"
1396 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1397 (d_t(j,i+nres),j=1,3)
1402 end subroutine velverlet_step
1403 !-----------------------------------------------------------------------------
1404 subroutine RESPA_step(itime)
1405 !-------------------------------------------------------------------------------
1406 ! Perform a single RESPA step.
1407 !-------------------------------------------------------------------------------
1408 ! implicit real*8 (a-h,o-z)
1409 ! include 'DIMENSIONS'
1413 use control, only:tcpu
1415 ! use io_conf, only:cartprint
1418 integer :: IERROR,ERRCODE
1420 ! include 'COMMON.SETUP'
1421 ! include 'COMMON.CONTROL'
1422 ! include 'COMMON.VAR'
1423 ! include 'COMMON.MD'
1425 ! include 'COMMON.LANGEVIN'
1427 ! include 'COMMON.LANGEVIN.lang0'
1429 ! include 'COMMON.CHAIN'
1430 ! include 'COMMON.DERIV'
1431 ! include 'COMMON.GEO'
1432 ! include 'COMMON.LOCAL'
1433 ! include 'COMMON.INTERACT'
1434 ! include 'COMMON.IOUNITS'
1435 ! include 'COMMON.NAMES'
1436 ! include 'COMMON.TIME1'
1437 real(kind=8),dimension(0:n_ene) :: energia_short,energia_long
1438 real(kind=8),dimension(3) :: L,vcm,incr
1439 real(kind=8),dimension(3,0:2*nres) :: dc_old0,d_t_old0,d_a_old0 !(3,0:maxres2) maxres2=2*maxres
1440 logical :: PRINT_AMTS_MSG = .false.
1441 integer :: count,rstcount !ilen,
1443 character(len=50) :: tytul
1444 integer :: maxcount_scale = 10
1445 !el common /gucio/ cm
1446 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1447 !el common /stochcalc/ stochforcvec
1448 integer :: itt,i,j,itsplit,itime
1450 !el common /cipiszcze/ itt
1452 real(kind=8) :: epdrift,tt0,epdriftmax
1455 if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres)) !(MAXRES6) maxres6=6*maxres
1459 if (large.and. mod(itime,ntwe).eq.0) then
1460 write (iout,*) "***************** RESPA itime",itime
1461 write (iout,*) "Cartesian and internal coordinates: step 0"
1463 call pdbout(0.0d0,"cipiszcze",iout)
1465 write (iout,*) "Accelerations from long-range forces"
1467 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1468 (d_a(j,i+nres),j=1,3)
1470 write (iout,*) "Velocities, step 0"
1472 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1473 (d_t(j,i+nres),j=1,3)
1478 ! Perform the initial RESPA step (increment velocities)
1479 ! write (iout,*) "*********************** RESPA ini"
1482 if (mod(itime,ntwe).eq.0 .and. large) then
1483 write (iout,*) "Velocities, end"
1485 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1486 (d_t(j,i+nres),j=1,3)
1490 ! Compute the short-range forces
1496 ! 7/2/2009 commented out
1498 ! call etotal_short(energia_short)
1501 ! 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
1504 d_a(j,i)=d_a_short(j,i)
1508 if (large.and. mod(itime,ntwe).eq.0) then
1509 write (iout,*) "energia_short",energia_short(0)
1510 write (iout,*) "Accelerations from short-range forces"
1512 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1513 (d_a(j,i+nres),j=1,3)
1518 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1520 t_enegrad=t_enegrad+tcpu()-tt0
1525 d_t_old(j,i)=d_t(j,i)
1526 d_a_old(j,i)=d_a(j,i)
1529 ! 6/30/08 A-MTS: attempt at increasing the split number
1532 dc_old0(j,i)=dc_old(j,i)
1533 d_t_old0(j,i)=d_t_old(j,i)
1534 d_a_old0(j,i)=d_a_old(j,i)
1537 if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
1538 if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
1545 ! write (iout,*) "itime",itime," ntime_split",ntime_split
1546 ! Split the time step
1547 d_time=d_time0/ntime_split
1548 ! Perform the short-range RESPA steps (velocity Verlet increments of
1549 ! positions and velocities using short-range forces)
1550 ! write (iout,*) "*********************** RESPA split"
1551 do itsplit=1,ntime_split
1554 else if (lang.eq.2 .or. lang.eq.3) then
1556 call stochastic_force(stochforcvec)
1559 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
1561 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1566 ! First step of the velocity Verlet algorithm
1571 else if (lang.eq.3) then
1573 call sd_verlet1_ciccotti
1575 else if (lang.eq.1) then
1580 ! Build the chain from the newly calculated coordinates
1581 call chainbuild_cart
1582 if (rattle) call rattle1
1584 if (large.and. mod(itime,ntwe).eq.0) then
1585 write (iout,*) "***** ITSPLIT",itsplit
1586 write (iout,*) "Cartesian and internal coordinates: step 1"
1587 call pdbout(0.0d0,"cipiszcze",iout)
1590 write (iout,*) "Velocities, step 1"
1592 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1593 (d_t(j,i+nres),j=1,3)
1602 ! Calculate energy and forces
1604 call etotal_short(energia_short)
1605 ! AL 4/17/17: Exit itime_split loop when energy goes infinite
1606 if (energia_short(0).gt.0.99e20 .or. isnan(energia_short(0)) ) then
1607 if (PRINT_AMTS_MSG) &
1608 write (iout,*) "Infinities/NaNs in energia_short",energia_short(0),"; increasing ntime_split to",ntime_split
1609 ntime_split=ntime_split*2
1610 if (ntime_split.gt.maxtime_split) then
1613 "Cannot rescue the run; aborting job. Retry with a smaller time step"
1615 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1618 "Cannot rescue the run; terminating. Retry with a smaller time step"
1624 if (large.and. mod(itime,ntwe).eq.0) &
1625 call enerprint(energia_short)
1628 t_eshort=t_eshort+MPI_Wtime()-tt0
1630 t_eshort=t_eshort+tcpu()-tt0
1634 ! Get the new accelerations
1636 ! 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
1639 d_a_short(j,i)=d_a(j,i)
1643 if (large.and. mod(itime,ntwe).eq.0) then
1644 write (iout,*)"energia_short",energia_short(0)
1645 write (iout,*) "Accelerations from short-range forces"
1647 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1648 (d_a(j,i+nres),j=1,3)
1653 ! Determine maximum acceleration and scale down the timestep if needed
1655 amax=amax/ntime_split**2
1656 call predict_edrift(epdrift)
1657 if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0) &
1658 write (iout,*) "amax",amax," damax",damax,&
1659 " epdrift",epdrift," epdriftmax",epdriftmax
1660 ! Exit loop and try with increased split number if the change of
1661 ! acceleration is too big
1662 if (amax.gt.damax .or. epdrift.gt.edriftmax) then
1663 if (ntime_split.lt.maxtime_split) then
1665 ntime_split=ntime_split*2
1666 ! AL 4/17/17: We should exit the itime_split loop when acceleration change is too big
1670 dc_old(j,i)=dc_old0(j,i)
1671 d_t_old(j,i)=d_t_old0(j,i)
1672 d_a_old(j,i)=d_a_old0(j,i)
1675 if (PRINT_AMTS_MSG) then
1676 write (iout,*) "acceleration/energy drift too large",amax,&
1677 epdrift," split increased to ",ntime_split," itime",itime,&
1683 "Uh-hu. Bumpy landscape. Maximum splitting number",&
1685 " already reached!!! Trying to carry on!"
1689 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1691 t_enegrad=t_enegrad+tcpu()-tt0
1693 ! Second step of the velocity Verlet algorithm
1698 else if (lang.eq.3) then
1700 call sd_verlet2_ciccotti
1702 else if (lang.eq.1) then
1707 if (rattle) call rattle2
1708 ! Backup the coordinates, velocities, and accelerations
1712 d_t_old(j,i)=d_t(j,i)
1713 d_a_old(j,i)=d_a(j,i)
1720 ! Restore the time step
1722 ! Compute long-range forces
1729 call etotal_long(energia_long)
1730 if (energia_long(0).gt.0.99e20 .or. isnan(energia_long(0))) then
1733 "Infinitied/NaNs in energia_long, Aborting MPI job."
1735 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1737 write (iout,*) "Infinitied/NaNs in energia_long, terminating."
1741 if (large.and. mod(itime,ntwe).eq.0) &
1742 call enerprint(energia_long)
1745 t_elong=t_elong+MPI_Wtime()-tt0
1747 t_elong=t_elong+tcpu()-tt0
1753 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1755 t_enegrad=t_enegrad+tcpu()-tt0
1757 ! Compute accelerations from long-range forces
1759 if (large.and. mod(itime,ntwe).eq.0) then
1760 write (iout,*) "energia_long",energia_long(0)
1761 write (iout,*) "Cartesian and internal coordinates: step 2"
1763 call pdbout(0.0d0,"cipiszcze",iout)
1765 write (iout,*) "Accelerations from long-range forces"
1767 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1768 (d_a(j,i+nres),j=1,3)
1770 write (iout,*) "Velocities, step 2"
1772 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1773 (d_t(j,i+nres),j=1,3)
1777 ! Compute the final RESPA step (increment velocities)
1778 ! write (iout,*) "*********************** RESPA fin"
1780 ! Compute the complete potential energy
1782 potEcomp(i)=energia_short(i)+energia_long(i)
1784 potE=potEcomp(0)-potEcomp(20)
1785 ! potE=energia_short(0)+energia_long(0)
1788 ! Calculate the kinetic and the total energy and the kinetic temperature
1791 ! Couple the system to Berendsen bath if needed
1792 if (tbf .and. lang.eq.0) then
1795 kinetic_T=2.0d0/(dimen3*Rb)*EK
1796 ! Backup the coordinates, velocities, and accelerations
1798 if (mod(itime,ntwe).eq.0 .and. large) then
1799 write (iout,*) "Velocities, end"
1801 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1802 (d_t(j,i+nres),j=1,3)
1807 end subroutine RESPA_step
1808 !-----------------------------------------------------------------------------
1809 subroutine RESPA_vel
1810 ! First and last RESPA step (incrementing velocities using long-range
1813 ! implicit real*8 (a-h,o-z)
1814 ! include 'DIMENSIONS'
1815 ! include 'COMMON.CONTROL'
1816 ! include 'COMMON.VAR'
1817 ! include 'COMMON.MD'
1818 ! include 'COMMON.CHAIN'
1819 ! include 'COMMON.DERIV'
1820 ! include 'COMMON.GEO'
1821 ! include 'COMMON.LOCAL'
1822 ! include 'COMMON.INTERACT'
1823 ! include 'COMMON.IOUNITS'
1824 ! include 'COMMON.NAMES'
1825 integer :: i,j,inres,mnum
1828 d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
1832 d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
1837 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
1838 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
1839 .and.(mnum.ne.5)) then
1842 d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
1847 end subroutine RESPA_vel
1848 !-----------------------------------------------------------------------------
1850 ! Applying velocity Verlet algorithm - step 1 to coordinates
1852 ! implicit real*8 (a-h,o-z)
1853 ! include 'DIMENSIONS'
1854 ! include 'COMMON.CONTROL'
1855 ! include 'COMMON.VAR'
1856 ! include 'COMMON.MD'
1857 ! include 'COMMON.CHAIN'
1858 ! include 'COMMON.DERIV'
1859 ! include 'COMMON.GEO'
1860 ! include 'COMMON.LOCAL'
1861 ! include 'COMMON.INTERACT'
1862 ! include 'COMMON.IOUNITS'
1863 ! include 'COMMON.NAMES'
1864 real(kind=8) :: adt,adt2
1865 integer :: i,j,inres,mnum
1868 write (iout,*) "VELVERLET1 START: DC"
1870 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),&
1871 (dc(j,i+nres),j=1,3)
1875 adt=d_a_old(j,0)*d_time
1877 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1878 d_t_new(j,0)=d_t_old(j,0)+adt2
1879 d_t(j,0)=d_t_old(j,0)+adt
1883 adt=d_a_old(j,i)*d_time
1885 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1886 d_t_new(j,i)=d_t_old(j,i)+adt2
1887 d_t(j,i)=d_t_old(j,i)+adt
1892 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
1893 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
1894 .and.(mnum.ne.5)) then
1897 adt=d_a_old(j,inres)*d_time
1899 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1900 d_t_new(j,inres)=d_t_old(j,inres)+adt2
1901 d_t(j,inres)=d_t_old(j,inres)+adt
1906 write (iout,*) "VELVERLET1 END: DC"
1908 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),&
1909 (dc(j,i+nres),j=1,3)
1913 end subroutine verlet1
1914 !-----------------------------------------------------------------------------
1916 ! Step 2 of the velocity Verlet algorithm: update velocities
1918 ! implicit real*8 (a-h,o-z)
1919 ! include 'DIMENSIONS'
1920 ! include 'COMMON.CONTROL'
1921 ! include 'COMMON.VAR'
1922 ! include 'COMMON.MD'
1923 ! include 'COMMON.CHAIN'
1924 ! include 'COMMON.DERIV'
1925 ! include 'COMMON.GEO'
1926 ! include 'COMMON.LOCAL'
1927 ! include 'COMMON.INTERACT'
1928 ! include 'COMMON.IOUNITS'
1929 ! include 'COMMON.NAMES'
1930 integer :: i,j,inres,mnum
1933 d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
1937 d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
1942 ! iti=iabs(itype(i,mnum))
1943 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
1944 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
1945 .and.(mnum.ne.5)) then
1948 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
1953 end subroutine verlet2
1954 !-----------------------------------------------------------------------------
1955 subroutine sddir_precalc
1956 ! Applying velocity Verlet algorithm - step 1 to coordinates
1957 ! implicit real*8 (a-h,o-z)
1958 ! include 'DIMENSIONS'
1964 ! include 'COMMON.CONTROL'
1965 ! include 'COMMON.VAR'
1966 ! include 'COMMON.MD'
1968 ! include 'COMMON.LANGEVIN'
1970 ! include 'COMMON.LANGEVIN.lang0'
1972 ! include 'COMMON.CHAIN'
1973 ! include 'COMMON.DERIV'
1974 ! include 'COMMON.GEO'
1975 ! include 'COMMON.LOCAL'
1976 ! include 'COMMON.INTERACT'
1977 ! include 'COMMON.IOUNITS'
1978 ! include 'COMMON.NAMES'
1979 ! include 'COMMON.TIME1'
1980 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1981 !el common /stochcalc/ stochforcvec
1982 real(kind=8) :: time00
1984 ! Compute friction and stochastic forces
1989 time_fric=time_fric+MPI_Wtime()-time00
1991 call stochastic_force(stochforcvec)
1992 time_stoch=time_stoch+MPI_Wtime()-time00
1995 ! Compute the acceleration due to friction forces (d_af_work) and stochastic
1996 ! forces (d_as_work)
1998 call ginv_mult(fric_work, d_af_work)
1999 call ginv_mult(stochforcvec, d_as_work)
2001 end subroutine sddir_precalc
2002 !-----------------------------------------------------------------------------
2003 subroutine sddir_verlet1
2004 ! Applying velocity Verlet algorithm - step 1 to velocities
2007 ! implicit real*8 (a-h,o-z)
2008 ! include 'DIMENSIONS'
2009 ! include 'COMMON.CONTROL'
2010 ! include 'COMMON.VAR'
2011 ! include 'COMMON.MD'
2013 ! include 'COMMON.LANGEVIN'
2015 ! include 'COMMON.LANGEVIN.lang0'
2017 ! include 'COMMON.CHAIN'
2018 ! include 'COMMON.DERIV'
2019 ! include 'COMMON.GEO'
2020 ! include 'COMMON.LOCAL'
2021 ! include 'COMMON.INTERACT'
2022 ! include 'COMMON.IOUNITS'
2023 ! include 'COMMON.NAMES'
2024 ! Revised 3/31/05 AL: correlation between random contributions to
2025 ! position and velocity increments included.
2026 real(kind=8) :: sqrt13 = 0.57735026918962576451d0 ! 1/sqrt(3)
2027 real(kind=8) :: adt,adt2
2028 integer :: i,j,ind,inres,mnum
2030 ! Add the contribution from BOTH friction and stochastic force to the
2031 ! coordinates, but ONLY the contribution from the friction forces to velocities
2034 adt=(d_a_old(j,0)+d_af_work(j))*d_time
2035 adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
2036 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
2037 d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
2038 d_t(j,0)=d_t_old(j,0)+adt
2043 adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
2044 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
2045 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
2046 d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
2047 d_t(j,i)=d_t_old(j,i)+adt
2053 ! iti=iabs(itype(i,mnum))
2054 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2055 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2056 .and.(mnum.ne.5)) then
2059 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
2060 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
2061 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
2062 d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
2063 d_t(j,inres)=d_t_old(j,inres)+adt
2069 end subroutine sddir_verlet1
2070 !-----------------------------------------------------------------------------
2071 subroutine sddir_verlet2
2072 ! Calculating the adjusted velocities for accelerations
2075 ! implicit real*8 (a-h,o-z)
2076 ! include 'DIMENSIONS'
2077 ! include 'COMMON.CONTROL'
2078 ! include 'COMMON.VAR'
2079 ! include 'COMMON.MD'
2081 ! include 'COMMON.LANGEVIN'
2083 ! include 'COMMON.LANGEVIN.lang0'
2085 ! include 'COMMON.CHAIN'
2086 ! include 'COMMON.DERIV'
2087 ! include 'COMMON.GEO'
2088 ! include 'COMMON.LOCAL'
2089 ! include 'COMMON.INTERACT'
2090 ! include 'COMMON.IOUNITS'
2091 ! include 'COMMON.NAMES'
2092 real(kind=8),dimension(6*nres) :: stochforcvec,d_as_work1 !(MAXRES6) maxres6=6*maxres
2093 real(kind=8) :: cos60 = 0.5d0, sin60 = 0.86602540378443864676d0
2094 integer :: i,j,ind,inres,mnum
2095 ! Revised 3/31/05 AL: correlation between random contributions to
2096 ! position and velocity increments included.
2097 ! The correlation coefficients are calculated at low-friction limit.
2098 ! Also, friction forces are now not calculated with new velocities.
2100 ! call friction_force
2101 call stochastic_force(stochforcvec)
2103 ! Compute the acceleration due to friction forces (d_af_work) and stochastic
2104 ! forces (d_as_work)
2106 call ginv_mult(stochforcvec, d_as_work1)
2112 d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j)) &
2113 +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
2118 d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j)) &
2119 +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
2125 ! iti=iabs(itype(i,mnum))
2126 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2127 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2128 .and.(mnum.ne.5)) then
2131 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres) &
2132 +d_af_work(ind+j))+sin60*d_as_work(ind+j) &
2133 +cos60*d_as_work1(ind+j))*d_time
2139 end subroutine sddir_verlet2
2140 !-----------------------------------------------------------------------------
2141 subroutine max_accel
2143 ! Find the maximum difference in the accelerations of the the sites
2144 ! at the beginning and the end of the time step.
2147 ! implicit real*8 (a-h,o-z)
2148 ! include 'DIMENSIONS'
2149 ! include 'COMMON.CONTROL'
2150 ! include 'COMMON.VAR'
2151 ! include 'COMMON.MD'
2152 ! include 'COMMON.CHAIN'
2153 ! include 'COMMON.DERIV'
2154 ! include 'COMMON.GEO'
2155 ! include 'COMMON.LOCAL'
2156 ! include 'COMMON.INTERACT'
2157 ! include 'COMMON.IOUNITS'
2158 real(kind=8),dimension(3) :: aux,accel,accel_old
2159 real(kind=8) :: dacc
2163 ! aux(j)=d_a(j,0)-d_a_old(j,0)
2164 accel_old(j)=d_a_old(j,0)
2171 ! 7/3/08 changed to asymmetric difference
2173 ! accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
2174 accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
2175 accel(j)=accel(j)+0.5d0*d_a(j,i)
2176 ! if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
2177 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
2178 dacc=dabs(accel(j)-accel_old(j))
2179 ! write (iout,*) i,dacc
2180 if (dacc.gt.amax) amax=dacc
2188 accel_old(j)=d_a_old(j,0)
2193 accel_old(j)=accel_old(j)+d_a_old(j,1)
2194 accel(j)=accel(j)+d_a(j,1)
2199 ! iti=iabs(itype(i,mnum))
2200 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2201 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2202 .and.(mnum.ne.5)) then
2204 ! accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
2205 accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
2206 accel(j)=accel(j)+d_a(j,i+nres)
2210 ! if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
2211 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
2212 dacc=dabs(accel(j)-accel_old(j))
2213 ! write (iout,*) "side-chain",i,dacc
2214 if (dacc.gt.amax) amax=dacc
2218 accel_old(j)=accel_old(j)+d_a_old(j,i)
2219 accel(j)=accel(j)+d_a(j,i)
2220 ! aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
2224 end subroutine max_accel
2225 !-----------------------------------------------------------------------------
2226 subroutine predict_edrift(epdrift)
2228 ! Predict the drift of the potential energy
2231 use control_data, only: lmuca
2232 ! implicit real*8 (a-h,o-z)
2233 ! include 'DIMENSIONS'
2234 ! include 'COMMON.CONTROL'
2235 ! include 'COMMON.VAR'
2236 ! include 'COMMON.MD'
2237 ! include 'COMMON.CHAIN'
2238 ! include 'COMMON.DERIV'
2239 ! include 'COMMON.GEO'
2240 ! include 'COMMON.LOCAL'
2241 ! include 'COMMON.INTERACT'
2242 ! include 'COMMON.IOUNITS'
2243 ! include 'COMMON.MUCA'
2244 real(kind=8) :: epdrift,epdriftij
2246 ! Drift of the potential energy
2252 epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
2253 if (lmuca) epdriftij=epdriftij*factor
2254 ! write (iout,*) "back",i,j,epdriftij
2255 if (epdriftij.gt.epdrift) epdrift=epdriftij
2259 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1.and.&
2260 molnum(i).ne.5) then
2263 dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
2264 if (lmuca) epdriftij=epdriftij*factor
2265 ! write (iout,*) "side",i,j,epdriftij
2266 if (epdriftij.gt.epdrift) epdrift=epdriftij
2270 epdrift=0.5d0*epdrift*d_time*d_time
2271 ! write (iout,*) "epdrift",epdrift
2273 end subroutine predict_edrift
2274 !-----------------------------------------------------------------------------
2275 subroutine verlet_bath
2277 ! Coupling to the thermostat by using the Berendsen algorithm
2280 ! implicit real*8 (a-h,o-z)
2281 ! include 'DIMENSIONS'
2282 ! include 'COMMON.CONTROL'
2283 ! include 'COMMON.VAR'
2284 ! include 'COMMON.MD'
2285 ! include 'COMMON.CHAIN'
2286 ! include 'COMMON.DERIV'
2287 ! include 'COMMON.GEO'
2288 ! include 'COMMON.LOCAL'
2289 ! include 'COMMON.INTERACT'
2290 ! include 'COMMON.IOUNITS'
2291 ! include 'COMMON.NAMES'
2292 real(kind=8) :: T_half,fact
2293 integer :: i,j,inres,mnum
2295 T_half=2.0d0/(dimen3*Rb)*EK
2296 fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
2297 ! write(iout,*) "T_half", T_half
2298 ! write(iout,*) "EK", EK
2299 ! write(iout,*) "fact", fact
2301 d_t(j,0)=fact*d_t(j,0)
2305 d_t(j,i)=fact*d_t(j,i)
2310 ! iti=iabs(itype(i,mnum))
2311 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2312 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2313 .and.(mnum.ne.5)) then
2316 d_t(j,inres)=fact*d_t(j,inres)
2321 end subroutine verlet_bath
2322 !-----------------------------------------------------------------------------
2324 ! Set up the initial conditions of a MD simulation
2327 use control, only:tcpu
2328 !el use io_basic, only:ilen
2331 use minimm, only:minim_dc,minimize,sc_move
2332 use io_config, only:readrst
2333 use io, only:statout
2334 use random, only: iran_num
2335 ! implicit real*8 (a-h,o-z)
2336 ! include 'DIMENSIONS'
2339 character(len=16) :: form
2340 integer :: IERROR,ERRCODE
2342 integer :: iranmin,itrial,itmp,n_model_try,k, &
2344 integer, dimension(:),allocatable :: list_model_try
2345 integer, dimension(0:nodes-1) :: i_start_models
2346 ! include 'COMMON.SETUP'
2347 ! include 'COMMON.CONTROL'
2348 ! include 'COMMON.VAR'
2349 ! include 'COMMON.MD'
2351 ! include 'COMMON.LANGEVIN'
2353 ! include 'COMMON.LANGEVIN.lang0'
2355 ! include 'COMMON.CHAIN'
2356 ! include 'COMMON.DERIV'
2357 ! include 'COMMON.GEO'
2358 ! include 'COMMON.LOCAL'
2359 ! include 'COMMON.INTERACT'
2360 ! include 'COMMON.IOUNITS'
2361 ! include 'COMMON.NAMES'
2362 ! include 'COMMON.REMD'
2363 real(kind=8),dimension(0:n_ene) :: energia_long,energia_short,energia
2364 real(kind=8),dimension(3) :: vcm,incr,L
2365 real(kind=8) :: xv,sigv,lowb,highb
2366 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
2367 character(len=256) :: qstr
2370 character(len=50) :: tytul
2371 logical :: file_exist
2372 !el common /gucio/ cm
2373 integer :: i,j,ipos,iq,iw,nft_sc,iretcode,nfun,ierr,mnum,itime
2374 real(kind=8) :: etot,tt0
2378 ! write(iout,*) "d_time", d_time
2379 ! Compute the standard deviations of stochastic forces for Langevin dynamics
2380 ! if the friction coefficients do not depend on surface area
2381 if (lang.gt.0 .and. .not.surfarea) then
2384 stdforcp(i)=stdfp(mnum)*dsqrt(gamp(mnum))
2388 stdforcsc(i)=stdfsc(iabs(itype(i,mnum)),mnum) &
2389 *dsqrt(gamsc(iabs(itype(i,mnum)),mnum))
2393 ! Open the pdb file for snapshotshots
2396 if (ilen(tmpdir).gt.0) &
2397 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2398 liczba(:ilen(liczba))//".pdb")
2400 file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2404 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) &
2405 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2406 liczba(:ilen(liczba))//".x")
2407 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2410 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) &
2411 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2412 liczba(:ilen(liczba))//".cx")
2413 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2419 if (ilen(tmpdir).gt.0) &
2420 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
2421 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
2423 if (ilen(tmpdir).gt.0) &
2424 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
2425 cartname=prefix(:ilen(prefix))//"_MD.cx"
2429 write (qstr,'(256(1h ))')
2432 iq = qinfrag(i,iset)*10
2433 iw = wfrag(i,iset)/100
2435 if(me.eq.king.or..not.out1file) &
2436 write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
2437 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
2442 iq = qinpair(i,iset)*10
2443 iw = wpair(i,iset)/100
2445 if(me.eq.king.or..not.out1file) &
2446 write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
2447 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
2451 ! pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
2453 ! cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
2455 ! cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
2457 ! statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
2461 if (restart1file) then
2463 inquire(file=mremd_rst_name,exist=file_exist)
2464 write (*,*) me," Before broadcast: file_exist",file_exist
2466 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,&
2469 write (*,*) me," After broadcast: file_exist",file_exist
2470 ! inquire(file=mremd_rst_name,exist=file_exist)
2471 if(me.eq.king.or..not.out1file) &
2472 write(iout,*) "Initial state read by master and distributed"
2474 if (ilen(tmpdir).gt.0) &
2475 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_' &
2476 //liczba(:ilen(liczba))//'.rst')
2477 inquire(file=rest2name,exist=file_exist)
2480 if(.not.restart1file) then
2481 if(me.eq.king.or..not.out1file) &
2482 write(iout,*) "Initial state will be read from file ",&
2483 rest2name(:ilen(rest2name))
2486 call rescale_weights(t_bath)
2488 if(me.eq.king.or..not.out1file)then
2489 if (restart1file) then
2490 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),&
2493 write(iout,*) "File ",rest2name(:ilen(rest2name)),&
2496 write(iout,*) "Initial velocities randomly generated"
2503 ! Generate initial velocities
2504 if(me.eq.king.or..not.out1file) &
2505 write(iout,*) "Initial velocities randomly generated"
2510 ! rest2name = prefix(:ilen(prefix))//'.rst'
2511 if(me.eq.king.or..not.out1file)then
2512 write (iout,*) "Initial velocities"
2514 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2515 (d_t(j,i+nres),j=1,3)
2517 ! Zeroing the total angular momentum of the system
2518 write(iout,*) "Calling the zero-angular momentum subroutine"
2521 ! Getting the potential energy and forces and velocities and accelerations
2523 ! write (iout,*) "velocity of the center of the mass:"
2524 ! write (iout,*) (vcm(j),j=1,3)
2526 d_t(j,0)=d_t(j,0)-vcm(j)
2528 ! Removing the velocity of the center of mass
2530 if(me.eq.king.or..not.out1file)then
2531 write (iout,*) "vcm right after adjustment:"
2532 write (iout,*) (vcm(j),j=1,3)
2537 if(iranconf.ne.0 .or.indpdb.gt.0.and..not.unres_pdb .or.preminim) then
2539 print *, 'Calling OVERLAP_SC'
2540 call overlap_sc(fail)
2541 print *,'after OVERLAP'
2544 print *,'call SC_MOVE'
2545 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
2546 print *,'SC_move',nft_sc,etot
2547 if(me.eq.king.or..not.out1file) &
2548 write(iout,*) 'SC_move',nft_sc,etot
2552 print *, 'Calling MINIM_DC'
2553 call minim_dc(etot,iretcode,nfun)
2555 call geom_to_var(nvar,varia)
2556 print *,'Calling MINIMIZE.'
2557 call minimize(etot,varia,iretcode,nfun)
2558 call var_to_geom(nvar,varia)
2560 if(me.eq.king.or..not.out1file) &
2561 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2563 if(iranconf.ne.0) then
2564 !c 8/22/17 AL Loop to produce a low-energy random conformation
2567 if(me.eq.king.or..not.out1file) &
2568 write (iout,*) 'Calling OVERLAP_SC'
2569 call overlap_sc(fail)
2570 endif !endif overlap
2573 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
2574 print *,'SC_move',nft_sc,etot
2575 if(me.eq.king.or..not.out1file) &
2576 write(iout,*) 'SC_move',nft_sc,etot
2580 print *, 'Calling MINIM_DC'
2581 call minim_dc(etot,iretcode,nfun)
2582 call int_from_cart1(.false.)
2584 call geom_to_var(nvar,varia)
2585 print *,'Calling MINIMIZE.'
2586 call minimize(etot,varia,iretcode,nfun)
2587 call var_to_geom(nvar,varia)
2589 if(me.eq.king.or..not.out1file) &
2590 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2592 if (isnan(etot) .or. etot.gt.4.0d6) then
2593 write (iout,*) "Energy too large",etot, &
2594 " trying another random conformation"
2597 call gen_rand_conf(itmp,*30)
2599 30 write (iout,*) 'Failed to generate random conformation', &
2601 write (*,*) 'Processor:',me, &
2602 ' Failed to generate random conformation',&
2611 write (iout,'(a,i3,a)') 'Processor:',me, &
2612 ' error in generating random conformation.'
2613 write (*,'(a,i3,a)') 'Processor:',me, &
2614 ' error in generating random conformation.'
2617 ! call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
2618 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
2628 write (iout,'(a,i3,a)') 'Processor:',me, &
2629 ' failed to generate a low-energy random conformation.'
2630 write (*,'(a,i3,a,f10.3)') 'Processor:',me, &
2631 ' failed to generate a low-energy random conformation.',etot
2635 ! call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
2636 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
2641 else if (preminim) then
2642 if (start_from_model) then
2646 do while (fail .and. n_model_try.lt.nmodel_start)
2647 write (iout,*) "n_model_try",n_model_try
2649 i_model=iran_num(1,nmodel_start)
2651 if (i_model.eq.list_model_try(k)) exit
2653 if (k.gt.n_model_try) exit
2655 n_model_try=n_model_try+1
2656 list_model_try(n_model_try)=i_model
2657 if (me.eq.king .or. .not. out1file) &
2658 write (iout,*) 'Trying to start from model ',&
2659 pdbfiles_chomo(i_model)(:ilen(pdbfiles_chomo(i_model)))
2662 c(j,i)=chomo(j,i,i_model)
2665 call int_from_cart(.true.,.false.)
2666 call sc_loc_geom(.false.)
2670 dc(j,i)=c(j,i+1)-c(j,i)
2671 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
2676 dc(j,i+nres)=c(j,i+nres)-c(j,i)
2677 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
2680 if (me.eq.king.or..not.out1file) then
2681 write (iout,*) "Energies before removing overlaps"
2682 call etotal(energia(0))
2683 call enerprint(energia(0))
2685 ! Remove SC overlaps if requested
2687 write (iout,*) 'Calling OVERLAP_SC'
2688 call overlap_sc(fail)
2691 "Failed to remove overlap from model",i_model
2695 if (me.eq.king.or..not.out1file) then
2696 write (iout,*) "Energies after removing overlaps"
2697 call etotal(energia(0))
2698 call enerprint(energia(0))
2701 ! Search for better SC rotamers if requested
2703 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
2704 print *,'SC_move',nft_sc,etot
2705 if (me.eq.king.or..not.out1file)&
2706 write(iout,*) 'SC_move',nft_sc,etot
2708 call etotal(energia(0))
2711 call MPI_Gather(i_model,1,MPI_INTEGER,i_start_models(0),&
2712 1,MPI_INTEGER,king,CG_COMM,IERROR)
2713 if (n_model_try.gt.nmodel_start .and.&
2714 (me.eq.king .or. out1file)) then
2716 "All models have irreparable overlaps. Trying randoms starts."
2718 i_model=nmodel_start+1
2722 ! Remove SC overlaps if requested
2724 write (iout,*) 'Calling OVERLAP_SC'
2725 call overlap_sc(fail)
2728 "Failed to remove overlap"
2731 if (me.eq.king.or..not.out1file) then
2732 write (iout,*) "Energies after removing overlaps"
2733 call etotal(energia(0))
2734 call enerprint(energia(0))
2737 ! 8/22/17 AL Minimize initial structure
2739 if (me.eq.king.or..not.out1file) write(iout,*)&
2740 'Minimizing initial PDB structure: Calling MINIM_DC'
2741 call minim_dc(etot,iretcode,nfun)
2743 call geom_to_var(nvar,varia)
2744 if(me.eq.king.or..not.out1file) write (iout,*)&
2745 'Minimizing initial PDB structure: Calling MINIMIZE.'
2746 call minimize(etot,varia,iretcode,nfun)
2747 call var_to_geom(nvar,varia)
2749 if (me.eq.king.or..not.out1file)&
2750 write(iout,*) 'LBFGS return code is ',status,' eval ',nfun
2751 if(me.eq.king.or..not.out1file)&
2752 write(iout,*) 'LBFGS return code is ',status,' eval ',nfun
2754 if (me.eq.king.or..not.out1file)&
2755 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2756 if(me.eq.king.or..not.out1file)&
2757 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2761 if (nmodel_start.gt.0 .and. me.eq.king) then
2762 write (iout,'(a)') "Task Starting model"
2764 if (i_start_models(i).gt.nmodel_start) then
2765 write (iout,'(i4,2x,a)') i,"RANDOM STRUCTURE"
2767 write(iout,'(i4,2x,a)')i,pdbfiles_chomo(i_start_models(i)) &
2768 (:ilen(pdbfiles_chomo(i_start_models(i))))
2773 call chainbuild_cart
2778 kinetic_T=2.0d0/(dimen3*Rb)*EK
2779 if(me.eq.king.or..not.out1file)then
2789 write(iout,*) "before ETOTAL"
2790 call etotal(potEcomp)
2791 if (large) call enerprint(potEcomp)
2794 t_etotal=t_etotal+MPI_Wtime()-tt0
2796 t_etotal=t_etotal+tcpu()-tt0
2803 if (amax*d_time .gt. dvmax) then
2804 d_time=d_time*dvmax/amax
2805 if(me.eq.king.or..not.out1file) write (iout,*) &
2806 "Time step reduced to",d_time,&
2807 " because of too large initial acceleration."
2809 if(me.eq.king.or..not.out1file)then
2810 write(iout,*) "Potential energy and its components"
2811 call enerprint(potEcomp)
2812 ! write(iout,*) (potEcomp(i),i=0,n_ene)
2814 potE=potEcomp(0)-potEcomp(20)
2818 if (ntwe.ne.0) call statout(itime)
2819 if(me.eq.king.or..not.out1file) &
2820 write (iout,'(/a/3(a25,1pe14.5/))') "Initial:", &
2821 " Kinetic energy",EK," Potential energy",potE, &
2822 " Total energy",totE," Maximum acceleration ", &
2825 write (iout,*) "Initial coordinates"
2827 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),&
2830 write (iout,*) "Initial dC"
2832 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),&
2833 (dc(j,i+nres),j=1,3)
2835 write (iout,*) "Initial velocities"
2836 write (iout,"(13x,' backbone ',23x,' side chain')")
2838 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2839 (d_t(j,i+nres),j=1,3)
2841 write (iout,*) "Initial accelerations"
2843 ! write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2844 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),&
2845 (d_a(j,i+nres),j=1,3)
2851 d_t_old(j,i)=d_t(j,i)
2852 d_a_old(j,i)=d_a(j,i)
2854 ! write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
2863 call etotal_short(energia_short)
2864 if (large) call enerprint(potEcomp)
2867 t_eshort=t_eshort+MPI_Wtime()-tt0
2869 t_eshort=t_eshort+tcpu()-tt0
2874 if(.not.out1file .and. large) then
2875 write (iout,*) "energia_long",energia_long(0),&
2876 " energia_short",energia_short(0),&
2877 " total",energia_long(0)+energia_short(0)
2878 write (iout,*) "Initial fast-force accelerations"
2880 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
2881 (d_a(j,i+nres),j=1,3)
2884 ! 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
2887 d_a_short(j,i)=d_a(j,i)
2896 call etotal_long(energia_long)
2897 if (large) call enerprint(potEcomp)
2900 t_elong=t_elong+MPI_Wtime()-tt0
2902 t_elong=t_elong+tcpu()-tt0
2907 if(.not.out1file .and. large) then
2908 write (iout,*) "energia_long",energia_long(0)
2909 write (iout,*) "Initial slow-force accelerations"
2911 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
2912 (d_a(j,i+nres),j=1,3)
2916 t_enegrad=t_enegrad+MPI_Wtime()-tt0
2918 t_enegrad=t_enegrad+tcpu()-tt0
2922 end subroutine init_MD
2923 !-----------------------------------------------------------------------------
2924 subroutine random_vel
2926 ! implicit real*8 (a-h,o-z)
2928 use random, only:anorm_distr
2930 ! include 'DIMENSIONS'
2931 ! include 'COMMON.CONTROL'
2932 ! include 'COMMON.VAR'
2933 ! include 'COMMON.MD'
2935 ! include 'COMMON.LANGEVIN'
2937 ! include 'COMMON.LANGEVIN.lang0'
2939 ! include 'COMMON.CHAIN'
2940 ! include 'COMMON.DERIV'
2941 ! include 'COMMON.GEO'
2942 ! include 'COMMON.LOCAL'
2943 ! include 'COMMON.INTERACT'
2944 ! include 'COMMON.IOUNITS'
2945 ! include 'COMMON.NAMES'
2946 ! include 'COMMON.TIME1'
2947 real(kind=8) :: xv,sigv,lowb,highb ,Ek1
2950 real(kind=8) ,allocatable, dimension(:) :: DDU1,DDU2,DL2,DL1,xsolv,DML,rs
2951 real(kind=8) :: sumx
2953 real(kind=8) ,allocatable, dimension(:) :: rsold
2954 real (kind=8),allocatable,dimension(:,:) :: matold
2958 integer :: i,j,ii,k,ind,mark,imark,mnum
2959 ! Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
2960 ! First generate velocities in the eigenspace of the G matrix
2961 ! write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
2968 sigv=dsqrt((Rb*t_bath)/geigen(i))
2971 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
2973 write (iout,*) "i",i," ii",ii," geigen",geigen(i),&
2974 " d_t_work_new",d_t_work_new(ii)
2985 Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
2988 write (iout,*) "Ek from eigenvectors",Ek1
2989 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
2993 ! Transform velocities to UNRES coordinate space
2994 allocate (DL1(2*nres))
2995 allocate (DDU1(2*nres))
2996 allocate (DL2(2*nres))
2997 allocate (DDU2(2*nres))
2998 allocate (xsolv(2*nres))
2999 allocate (DML(2*nres))
3000 allocate (rs(2*nres))
3002 allocate (rsold(2*nres))
3003 allocate (matold(2*nres,2*nres))
3005 matold(1,1)=DMorig(1)
3006 matold(1,2)=DU1orig(1)
3007 matold(1,3)=DU2orig(1)
3008 write (*,*) DMorig(1),DU1orig(1),DU2orig(1)
3013 matold(i,j)=DMorig(i)
3014 matold(i,j-1)=DU1orig(i-1)
3016 matold(i,j-2)=DU2orig(i-2)
3024 matold(i,j+1)=DU1orig(i)
3030 matold(i,j+2)=DU2orig(i)
3034 matold(dimen,dimen)=DMorig(dimen)
3035 matold(dimen,dimen-1)=DU1orig(dimen-1)
3036 matold(dimen,dimen-2)=DU2orig(dimen-2)
3037 write(iout,*) "old gmatrix"
3038 call matout(dimen,dimen,2*nres,2*nres,matold)
3042 ! Find the ith eigenvector of the pentadiagonal inertiq matrix
3046 DML(j)=DMorig(j)-geigen(i)
3049 DML(j-1)=DMorig(j)-geigen(i)
3054 DDU1(imark-1)=DU2orig(imark-1)
3055 do j=imark+1,dimen-1
3056 DDU1(j-1)=DU1orig(j)
3064 DDU2(j)=DU2orig(j+1)
3073 write (iout,*) "DL2,DL1,DML,DDU1,DDU2"
3074 write(iout,'(10f10.5)') (DL2(k),k=3,dimen-1)
3075 write(iout,'(10f10.5)') (DL1(k),k=2,dimen-1)
3076 write(iout,'(10f10.5)') (DML(k),k=1,dimen-1)
3077 write(iout,'(10f10.5)') (DDU1(k),k=1,dimen-2)
3078 write(iout,'(10f10.5)') (DDU2(k),k=1,dimen-3)
3081 if (imark.gt.2) rs(imark-2)=-DU2orig(imark-2)
3082 if (imark.gt.1) rs(imark-1)=-DU1orig(imark-1)
3083 if (imark.lt.dimen) rs(imark)=-DU1orig(imark)
3084 if (imark.lt.dimen-1) rs(imark+1)=-DU2orig(imark)
3088 ! write (iout,*) "Vector rs"
3090 ! write (iout,*) j,rs(j)
3093 call FDIAG(dimen-1,DL2,DL1,DML,DDU1,DDU2,rs,xsolv,mark)
3100 sumx=-geigen(i)*xsolv(j)
3102 sumx=sumx+matold(j,k)*xsolv(k)
3105 sumx=sumx+matold(j,k)*xsolv(k-1)
3107 write(iout,'(i5,3f10.5)') j,sumx,rsold(j),sumx-rsold(j)
3110 sumx=-geigen(i)*xsolv(j-1)
3112 sumx=sumx+matold(j,k)*xsolv(k)
3115 sumx=sumx+matold(j,k)*xsolv(k-1)
3117 write(iout,'(i5,3f10.5)') j-1,sumx,rsold(j-1),sumx-rsold(j-1)
3121 "Solution of equations system",i," for eigenvalue",geigen(i)
3123 write(iout,'(i5,f10.5)') j,xsolv(j)
3126 do j=dimen-1,imark,-1
3131 write (iout,*) "Un-normalized eigenvector",i," for eigenvalue",geigen(i)
3133 write(iout,'(i5,f10.5)') j,xsolv(j)
3136 ! Normalize ith eigenvector
3139 sumx=sumx+xsolv(j)**2
3143 xsolv(j)=xsolv(j)/sumx
3146 write (iout,*) "Eigenvector",i," for eigenvalue",geigen(i)
3148 write(iout,'(i5,3f10.5)') j,xsolv(j)
3151 ! All done at this point for eigenvector i, exit loop
3159 write (iout,*) "Unable to find eigenvector",i
3162 ! write (iout,*) "i=",i
3164 ! write (iout,*) "k=",k
3167 ! write(iout,*) "ind",ind," ind1",3*(i-1)+k
3168 d_t_work(ind)=d_t_work(ind) &
3169 +xsolv(j)*d_t_work_new(3*(i-1)+k)
3172 enddo ! i (loop over the eigenvectors)
3175 write (iout,*) "d_t_work"
3177 write (iout,"(i5,f10.5)") i,d_t_work(i)
3182 ! if (itype(i,1).eq.10) then
3184 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
3185 iti=iabs(itype(i,mnum))
3186 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3187 if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
3188 .or.(mnum.eq.5)) then
3195 ! write (iout,*) "k",k," ii+k",ii+k," ii+j+k",ii+j+k,"EK1",Ek1
3196 Ek1=Ek1+0.5d0*mp(mnum)*((d_t_work(ii+j+k)+d_t_work_new(ii+k))/2)**2+&
3197 0.5d0*0.25d0*IP(mnum)*(d_t_work(ii+j+k)-d_t_work_new(ii+k))**2
3200 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3201 .and.(mnum.ne.5)) ii=ii+3
3202 write (iout,*) "i",i," itype",itype(i,mnum)," mass",msc(itype(i,mnum),mnum)
3203 write (iout,*) "ii",ii
3206 write (iout,*) "k",k," ii",ii,"EK1",EK1
3207 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3209 Ek1=Ek1+0.5d0*Isc(iabs(itype(i,mnum)),mnum)*(d_t_work(ii)-d_t_work(ii-3))**2
3210 Ek1=Ek1+0.5d0*msc(iabs(itype(i,mnum)),mnum)*d_t_work(ii)**2
3212 write (iout,*) "i",i," ii",ii
3214 write (iout,*) "Ek from d_t_work",Ek1
3215 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
3217 if(allocated(DDU1)) deallocate(DDU1)
3218 if(allocated(DDU2)) deallocate(DDU2)
3219 if(allocated(DL2)) deallocate(DL2)
3220 if(allocated(DL1)) deallocate(DL1)
3221 if(allocated(xsolv)) deallocate(xsolv)
3222 if(allocated(DML)) deallocate(DML)
3223 if(allocated(rs)) deallocate(rs)
3225 if(allocated(matold)) deallocate(matold)
3226 if(allocated(rsold)) deallocate(rsold)
3231 d_t(k,j)=d_t_work(ind)
3235 if (itype(j,1).ne.10 .and. itype(j,mnum).ne.ntyp1_molec(mnum)&
3236 .and.(mnum.ne.5)) then
3238 d_t(k,j+nres)=d_t_work(ind)
3244 write (iout,*) "Random velocities in the Calpha,SC space"
3246 write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3)
3249 write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3)
3256 ! if (itype(i,1).eq.10) then
3258 if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
3259 .or.(mnum.eq.5)) then
3261 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
3265 d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
3266 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
3271 write (iout,*)"Random velocities in the virtual-bond-vector space"
3273 write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3)
3276 write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3)
3279 write (iout,*) "Ek from d_t_work",Ek1
3280 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
3288 d_t_work(ind)=d_t_work(ind) &
3289 +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
3291 ! write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
3295 ! Transfer to the d_t vector
3297 d_t(j,0)=d_t_work(j)
3303 d_t(j,i)=d_t_work(ind)
3308 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3309 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3310 .and.(mnum.ne.5)) then
3313 d_t(j,i+nres)=d_t_work(ind)
3319 ! write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",&
3320 ! 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
3322 ! write(iout,*) "end init MD"
3324 end subroutine random_vel
3325 !-----------------------------------------------------------------------------
3327 subroutine sd_verlet_p_setup
3328 ! Sets up the parameters of stochastic Verlet algorithm
3329 ! implicit real*8 (a-h,o-z)
3330 ! include 'DIMENSIONS'
3331 use control, only: tcpu
3336 ! include 'COMMON.CONTROL'
3337 ! include 'COMMON.VAR'
3338 ! include 'COMMON.MD'
3340 ! include 'COMMON.LANGEVIN'
3342 ! include 'COMMON.LANGEVIN.lang0'
3344 ! include 'COMMON.CHAIN'
3345 ! include 'COMMON.DERIV'
3346 ! include 'COMMON.GEO'
3347 ! include 'COMMON.LOCAL'
3348 ! include 'COMMON.INTERACT'
3349 ! include 'COMMON.IOUNITS'
3350 ! include 'COMMON.NAMES'
3351 ! include 'COMMON.TIME1'
3352 real(kind=8),dimension(6*nres) :: emgdt !(MAXRES6) maxres6=6*maxres
3353 real(kind=8) :: pterm,vterm,rho,rhoc,vsig
3354 real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
3355 prand_vec,vrand_vec1,vrand_vec2 !(MAXRES6) maxres6=6*maxres
3356 logical :: lprn = .false.
3357 real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
3358 real(kind=8) :: ktm,gdt,egdt,gdt2,gdt3,gdt4,gdt5,gdt6,gdt7,gdt8,&
3360 integer :: i,maxres2
3367 ! AL 8/17/04 Code adapted from tinker
3369 ! Get the frictional and random terms for stochastic dynamics in the
3370 ! eigenspace of mass-scaled UNRES friction matrix
3374 gdt = fricgam(i) * d_time
3376 ! Stochastic dynamics reduces to simple MD for zero friction
3378 if (gdt .le. zero) then
3379 pfric_vec(i) = 1.0d0
3380 vfric_vec(i) = d_time
3381 afric_vec(i) = 0.5d0 * d_time * d_time
3382 prand_vec(i) = 0.0d0
3383 vrand_vec1(i) = 0.0d0
3384 vrand_vec2(i) = 0.0d0
3386 ! Analytical expressions when friction coefficient is large
3389 if (gdt .ge. gdt_radius) then
3392 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
3393 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
3394 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
3395 vterm = 1.0d0 - egdt**2
3396 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
3398 ! Use series expansions when friction coefficient is small
3409 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0 &
3410 - gdt5/120.0d0 + gdt6/720.0d0 &
3411 - gdt7/5040.0d0 + gdt8/40320.0d0 &
3412 - gdt9/362880.0d0) / fricgam(i)**2
3413 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
3414 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
3415 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0 &
3416 + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0 &
3417 + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0 &
3418 + 127.0d0*gdt9/90720.0d0
3419 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0 &
3420 - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0 &
3421 - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0 &
3422 - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
3423 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0 &
3424 - 17.0d0*gdt2/1280.0d0 &
3425 + 17.0d0*gdt3/6144.0d0 &
3426 + 40967.0d0*gdt4/34406400.0d0 &
3427 - 57203.0d0*gdt5/275251200.0d0 &
3428 - 1429487.0d0*gdt6/13212057600.0d0)
3431 ! Compute the scaling factors of random terms for the nonzero friction case
3433 ktm = 0.5d0*d_time/fricgam(i)
3434 psig = dsqrt(ktm*pterm) / fricgam(i)
3435 vsig = dsqrt(ktm*vterm)
3436 rhoc = dsqrt(1.0d0 - rho*rho)
3438 vrand_vec1(i) = vsig * rho
3439 vrand_vec2(i) = vsig * rhoc
3444 "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",&
3447 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),&
3448 afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
3452 ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
3455 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
3456 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
3457 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
3458 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
3459 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
3460 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
3463 t_sdsetup=t_sdsetup+MPI_Wtime()
3465 t_sdsetup=t_sdsetup+tcpu()-tt0
3468 end subroutine sd_verlet_p_setup
3469 !-----------------------------------------------------------------------------
3470 subroutine eigtransf1(n,ndim,ab,d,c)
3474 real(kind=8) :: ab(ndim,ndim,n),c(ndim,n),d(ndim)
3480 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
3485 end subroutine eigtransf1
3486 !-----------------------------------------------------------------------------
3487 subroutine eigtransf(n,ndim,a,b,d,c)
3491 real(kind=8) :: a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
3497 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
3502 end subroutine eigtransf
3503 !-----------------------------------------------------------------------------
3504 subroutine sd_verlet1
3506 ! Applying stochastic velocity Verlet algorithm - step 1 to velocities
3508 ! implicit real*8 (a-h,o-z)
3509 ! include 'DIMENSIONS'
3510 ! include 'COMMON.CONTROL'
3511 ! include 'COMMON.VAR'
3512 ! include 'COMMON.MD'
3514 ! include 'COMMON.LANGEVIN'
3516 ! include 'COMMON.LANGEVIN.lang0'
3518 ! include 'COMMON.CHAIN'
3519 ! include 'COMMON.DERIV'
3520 ! include 'COMMON.GEO'
3521 ! include 'COMMON.LOCAL'
3522 ! include 'COMMON.INTERACT'
3523 ! include 'COMMON.IOUNITS'
3524 ! include 'COMMON.NAMES'
3525 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
3526 !el common /stochcalc/ stochforcvec
3527 logical :: lprn = .false.
3528 real(kind=8) :: ddt1,ddt2
3529 integer :: i,j,ind,inres
3531 ! write (iout,*) "dc_old"
3533 ! write (iout,'(i5,3f10.5,5x,3f10.5)')
3534 ! & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3537 dc_work(j)=dc_old(j,0)
3538 d_t_work(j)=d_t_old(j,0)
3539 d_a_work(j)=d_a_old(j,0)
3544 dc_work(ind+j)=dc_old(j,i)
3545 d_t_work(ind+j)=d_t_old(j,i)
3546 d_a_work(ind+j)=d_a_old(j,i)
3552 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3553 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3554 .and.(mnum.ne.5)) then
3556 dc_work(ind+j)=dc_old(j,i+nres)
3557 d_t_work(ind+j)=d_t_old(j,i+nres)
3558 d_a_work(ind+j)=d_a_old(j,i+nres)
3566 "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",&
3570 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
3571 vfric_mat(i,j),afric_mat(i,j),&
3572 prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3580 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) &
3581 +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3582 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3583 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3585 d_t_work_new(i)=ddt1+0.5d0*ddt2
3586 d_t_work(i)=ddt1+ddt2
3591 d_t(j,0)=d_t_work(j)
3596 dc(j,i)=dc_work(ind+j)
3597 d_t(j,i)=d_t_work(ind+j)
3603 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3604 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3605 .and.(mnum.ne.5)) then
3608 dc(j,inres)=dc_work(ind+j)
3609 d_t(j,inres)=d_t_work(ind+j)
3615 end subroutine sd_verlet1
3616 !-----------------------------------------------------------------------------
3617 subroutine sd_verlet2
3619 ! Calculating the adjusted velocities for accelerations
3621 ! implicit real*8 (a-h,o-z)
3622 ! include 'DIMENSIONS'
3623 ! include 'COMMON.CONTROL'
3624 ! include 'COMMON.VAR'
3625 ! include 'COMMON.MD'
3627 ! include 'COMMON.LANGEVIN'
3629 ! include 'COMMON.LANGEVIN.lang0'
3631 ! include 'COMMON.CHAIN'
3632 ! include 'COMMON.DERIV'
3633 ! include 'COMMON.GEO'
3634 ! include 'COMMON.LOCAL'
3635 ! include 'COMMON.INTERACT'
3636 ! include 'COMMON.IOUNITS'
3637 ! include 'COMMON.NAMES'
3638 !el real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV !(MAXRES6) maxres6=6*maxres
3639 real(kind=8),dimension(6*nres) :: stochforcvecV !(MAXRES6) maxres6=6*maxres
3640 !el common /stochcalc/ stochforcvec
3642 real(kind=8) :: ddt1,ddt2
3643 integer :: i,j,ind,inres
3644 ! Compute the stochastic forces which contribute to velocity change
3646 call stochastic_force(stochforcvecV)
3653 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3654 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+ &
3655 vrand_mat2(i,j)*stochforcvecV(j)
3657 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
3661 d_t(j,0)=d_t_work(j)
3666 d_t(j,i)=d_t_work(ind+j)
3672 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3673 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3674 .and.(mnum.ne.5)) then
3677 d_t(j,inres)=d_t_work(ind+j)
3683 end subroutine sd_verlet2
3684 !-----------------------------------------------------------------------------
3685 subroutine sd_verlet_ciccotti_setup
3687 ! Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
3689 ! implicit real*8 (a-h,o-z)
3690 ! include 'DIMENSIONS'
3691 use control, only: tcpu
3696 ! include 'COMMON.CONTROL'
3697 ! include 'COMMON.VAR'
3698 ! include 'COMMON.MD'
3700 ! include 'COMMON.LANGEVIN'
3702 ! include 'COMMON.LANGEVIN.lang0'
3704 ! include 'COMMON.CHAIN'
3705 ! include 'COMMON.DERIV'
3706 ! include 'COMMON.GEO'
3707 ! include 'COMMON.LOCAL'
3708 ! include 'COMMON.INTERACT'
3709 ! include 'COMMON.IOUNITS'
3710 ! include 'COMMON.NAMES'
3711 ! include 'COMMON.TIME1'
3712 real(kind=8),dimension(6*nres) :: emgdt !(MAXRES6) maxres6=6*maxres
3713 real(kind=8) :: pterm,vterm,rho,rhoc,vsig
3714 real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
3715 prand_vec,vrand_vec1,vrand_vec2 !(MAXRES6) maxres6=6*maxres
3716 logical :: lprn = .false.
3717 real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
3718 real(kind=8) :: ktm,gdt,egdt,tt0
3719 integer :: i,maxres2
3726 ! AL 8/17/04 Code adapted from tinker
3728 ! Get the frictional and random terms for stochastic dynamics in the
3729 ! eigenspace of mass-scaled UNRES friction matrix
3733 write (iout,*) "i",i," fricgam",fricgam(i)
3734 gdt = fricgam(i) * d_time
3736 ! Stochastic dynamics reduces to simple MD for zero friction
3738 if (gdt .le. zero) then
3739 pfric_vec(i) = 1.0d0
3740 vfric_vec(i) = d_time
3741 afric_vec(i) = 0.5d0*d_time*d_time
3742 prand_vec(i) = afric_vec(i)
3743 vrand_vec2(i) = vfric_vec(i)
3745 ! Analytical expressions when friction coefficient is large
3750 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
3751 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
3752 prand_vec(i) = afric_vec(i)
3753 vrand_vec2(i) = vfric_vec(i)
3755 ! Compute the scaling factors of random terms for the nonzero friction case
3757 ! ktm = 0.5d0*d_time/fricgam(i)
3758 ! psig = dsqrt(ktm*pterm) / fricgam(i)
3759 ! vsig = dsqrt(ktm*vterm)
3760 ! prand_vec(i) = psig*afric_vec(i)
3761 ! vrand_vec2(i) = vsig*vfric_vec(i)
3766 "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",&
3769 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),&
3770 afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
3774 ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
3776 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
3777 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
3778 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
3779 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
3780 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
3782 t_sdsetup=t_sdsetup+MPI_Wtime()
3784 t_sdsetup=t_sdsetup+tcpu()-tt0
3787 end subroutine sd_verlet_ciccotti_setup
3788 !-----------------------------------------------------------------------------
3789 subroutine sd_verlet1_ciccotti
3791 ! Applying stochastic velocity Verlet algorithm - step 1 to velocities
3792 ! implicit real*8 (a-h,o-z)
3794 ! include 'DIMENSIONS'
3798 ! include 'COMMON.CONTROL'
3799 ! include 'COMMON.VAR'
3800 ! include 'COMMON.MD'
3802 ! include 'COMMON.LANGEVIN'
3804 ! include 'COMMON.LANGEVIN.lang0'
3806 ! include 'COMMON.CHAIN'
3807 ! include 'COMMON.DERIV'
3808 ! include 'COMMON.GEO'
3809 ! include 'COMMON.LOCAL'
3810 ! include 'COMMON.INTERACT'
3811 ! include 'COMMON.IOUNITS'
3812 ! include 'COMMON.NAMES'
3813 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
3814 !el common /stochcalc/ stochforcvec
3815 logical :: lprn = .false.
3816 real(kind=8) :: ddt1,ddt2
3817 integer :: i,j,ind,inres
3818 ! write (iout,*) "dc_old"
3820 ! write (iout,'(i5,3f10.5,5x,3f10.5)')
3821 ! & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3824 dc_work(j)=dc_old(j,0)
3825 d_t_work(j)=d_t_old(j,0)
3826 d_a_work(j)=d_a_old(j,0)
3831 dc_work(ind+j)=dc_old(j,i)
3832 d_t_work(ind+j)=d_t_old(j,i)
3833 d_a_work(ind+j)=d_a_old(j,i)
3838 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3840 dc_work(ind+j)=dc_old(j,i+nres)
3841 d_t_work(ind+j)=d_t_old(j,i+nres)
3842 d_a_work(ind+j)=d_a_old(j,i+nres)
3851 "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",&
3855 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
3856 vfric_mat(i,j),afric_mat(i,j),&
3857 prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3865 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) &
3866 +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3867 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3868 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3870 d_t_work_new(i)=ddt1+0.5d0*ddt2
3871 d_t_work(i)=ddt1+ddt2
3876 d_t(j,0)=d_t_work(j)
3881 dc(j,i)=dc_work(ind+j)
3882 d_t(j,i)=d_t_work(ind+j)
3888 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3889 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3890 .and.(mnum.ne.5)) then
3893 dc(j,inres)=dc_work(ind+j)
3894 d_t(j,inres)=d_t_work(ind+j)
3900 end subroutine sd_verlet1_ciccotti
3901 !-----------------------------------------------------------------------------
3902 subroutine sd_verlet2_ciccotti
3904 ! Calculating the adjusted velocities for accelerations
3906 ! implicit real*8 (a-h,o-z)
3907 ! include 'DIMENSIONS'
3908 ! include 'COMMON.CONTROL'
3909 ! include 'COMMON.VAR'
3910 ! include 'COMMON.MD'
3912 ! include 'COMMON.LANGEVIN'
3914 ! include 'COMMON.LANGEVIN.lang0'
3916 ! include 'COMMON.CHAIN'
3917 ! include 'COMMON.DERIV'
3918 ! include 'COMMON.GEO'
3919 ! include 'COMMON.LOCAL'
3920 ! include 'COMMON.INTERACT'
3921 ! include 'COMMON.IOUNITS'
3922 ! include 'COMMON.NAMES'
3923 !el real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV !(MAXRES6) maxres6=6*maxres
3924 real(kind=8),dimension(6*nres) :: stochforcvecV !(MAXRES6) maxres6=6*maxres
3925 !el common /stochcalc/ stochforcvec
3926 real(kind=8) :: ddt1,ddt2
3927 integer :: i,j,ind,inres
3929 ! Compute the stochastic forces which contribute to velocity change
3931 call stochastic_force(stochforcvecV)
3938 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3939 ! ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
3940 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
3942 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
3946 d_t(j,0)=d_t_work(j)
3951 d_t(j,i)=d_t_work(ind+j)
3957 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3959 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3962 d_t(j,inres)=d_t_work(ind+j)
3968 end subroutine sd_verlet2_ciccotti
3970 !-----------------------------------------------------------------------------
3972 !-----------------------------------------------------------------------------
3973 subroutine inertia_tensor
3975 ! Calculating the intertia tensor for the entire protein in order to
3976 ! remove the perpendicular components of velocity matrix which cause
3977 ! the molecule to rotate.
3980 ! implicit real*8 (a-h,o-z)
3981 ! include 'DIMENSIONS'
3982 ! include 'COMMON.CONTROL'
3983 ! include 'COMMON.VAR'
3984 ! include 'COMMON.MD'
3985 ! include 'COMMON.CHAIN'
3986 ! include 'COMMON.DERIV'
3987 ! include 'COMMON.GEO'
3988 ! include 'COMMON.LOCAL'
3989 ! include 'COMMON.INTERACT'
3990 ! include 'COMMON.IOUNITS'
3991 ! include 'COMMON.NAMES'
3993 real(kind=8),dimension(3,3) :: Im,Imcp,eigvec,Id
3994 real(kind=8),dimension(3) :: pr,eigval,L,vp,vrot
3995 real(kind=8) :: M_SC,mag,mag2,M_PEP
3996 real(kind=8),dimension(3,0:nres) :: vpp !(3,0:MAXRES)
3997 real(kind=8),dimension(3) :: vs_p,pp,incr,v
3998 real(kind=8),dimension(3,3) :: pr1,pr2
4000 !el common /gucio/ cm
4001 integer :: iti,inres,i,j,k,mnum
4012 ! calculating the center of the mass of the protein
4016 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
4017 if (itype(i,mnum).eq.ntyp1_molec(mnum)) cycle
4018 M_PEP=M_PEP+mp(mnum)
4020 cm(j)=cm(j)+(c(j,i)+0.5d0*dc(j,i))*mp(mnum)
4029 if (mnum.eq.5) cycle
4030 iti=iabs(itype(i,mnum))
4031 M_SC=M_SC+msc(iabs(iti),mnum)
4034 cm(j)=cm(j)+msc(iabs(iti),mnum)*c(j,inres)
4038 cm(j)=cm(j)/(M_SC+M_PEP)
4043 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
4045 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
4047 Im(1,1)=Im(1,1)+mp(mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
4048 Im(1,2)=Im(1,2)-mp(mnum)*pr(1)*pr(2)
4049 Im(1,3)=Im(1,3)-mp(mnum)*pr(1)*pr(3)
4050 Im(2,3)=Im(2,3)-mp(mnum)*pr(2)*pr(3)
4051 Im(2,2)=Im(2,2)+mp(mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
4052 Im(3,3)=Im(3,3)+mp(mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
4057 iti=iabs(itype(i,mnum))
4058 if (mnum.eq.5) cycle
4061 pr(j)=c(j,inres)-cm(j)
4063 Im(1,1)=Im(1,1)+msc(iabs(iti),mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
4064 Im(1,2)=Im(1,2)-msc(iabs(iti),mnum)*pr(1)*pr(2)
4065 Im(1,3)=Im(1,3)-msc(iabs(iti),mnum)*pr(1)*pr(3)
4066 Im(2,3)=Im(2,3)-msc(iabs(iti),mnum)*pr(2)*pr(3)
4067 Im(2,2)=Im(2,2)+msc(iabs(iti),mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
4068 Im(3,3)=Im(3,3)+msc(iabs(iti),mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
4073 Im(1,1)=Im(1,1)+Ip(mnum)*(1-dc_norm(1,i)*dc_norm(1,i))* &
4074 vbld(i+1)*vbld(i+1)*0.25d0
4075 Im(1,2)=Im(1,2)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(2,i))* &
4076 vbld(i+1)*vbld(i+1)*0.25d0
4077 Im(1,3)=Im(1,3)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(3,i))* &
4078 vbld(i+1)*vbld(i+1)*0.25d0
4079 Im(2,3)=Im(2,3)+Ip(mnum)*(-dc_norm(2,i)*dc_norm(3,i))* &
4080 vbld(i+1)*vbld(i+1)*0.25d0
4081 Im(2,2)=Im(2,2)+Ip(mnum)*(1-dc_norm(2,i)*dc_norm(2,i))* &
4082 vbld(i+1)*vbld(i+1)*0.25d0
4083 Im(3,3)=Im(3,3)+Ip(mnum)*(1-dc_norm(3,i)*dc_norm(3,i))* &
4084 vbld(i+1)*vbld(i+1)*0.25d0
4090 ! if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
4091 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4092 .and.(mnum.ne.5)) then
4093 iti=iabs(itype(i,mnum))
4095 Im(1,1)=Im(1,1)+Isc(iti,mnum)*(1-dc_norm(1,inres)* &
4096 dc_norm(1,inres))*vbld(inres)*vbld(inres)
4097 Im(1,2)=Im(1,2)-Isc(iti,mnum)*(dc_norm(1,inres)* &
4098 dc_norm(2,inres))*vbld(inres)*vbld(inres)
4099 Im(1,3)=Im(1,3)-Isc(iti,mnum)*(dc_norm(1,inres)* &
4100 dc_norm(3,inres))*vbld(inres)*vbld(inres)
4101 Im(2,3)=Im(2,3)-Isc(iti,mnum)*(dc_norm(2,inres)* &
4102 dc_norm(3,inres))*vbld(inres)*vbld(inres)
4103 Im(2,2)=Im(2,2)+Isc(iti,mnum)*(1-dc_norm(2,inres)* &
4104 dc_norm(2,inres))*vbld(inres)*vbld(inres)
4105 Im(3,3)=Im(3,3)+Isc(iti,mnum)*(1-dc_norm(3,inres)* &
4106 dc_norm(3,inres))*vbld(inres)*vbld(inres)
4111 ! write(iout,*) "The angular momentum before adjustment:"
4112 ! write(iout,*) (L(j),j=1,3)
4118 ! Copying the Im matrix for the djacob subroutine
4126 ! Finding the eigenvectors and eignvalues of the inertia tensor
4127 call djacob(3,3,10000,1.0d-10,Imcp,eigvec,eigval)
4128 ! write (iout,*) "Eigenvalues & Eigenvectors"
4129 ! write (iout,'(5x,3f10.5)') (eigval(i),i=1,3)
4132 ! write (iout,'(i5,3f10.5)') i,(eigvec(i,j),j=1,3)
4134 ! Constructing the diagonalized matrix
4136 if (dabs(eigval(i)).gt.1.0d-15) then
4137 Id(i,i)=1.0d0/eigval(i)
4144 Imcp(i,j)=eigvec(j,i)
4150 pr1(i,j)=pr1(i,j)+Id(i,k)*Imcp(k,j)
4157 pr2(i,j)=pr2(i,j)+eigvec(i,k)*pr1(k,j)
4161 ! Calculating the total rotational velocity of the molecule
4164 vrot(i)=vrot(i)+pr2(i,j)*L(j)
4167 ! Resetting the velocities
4169 call vecpr(vrot(1),dc(1,i),vp)
4171 d_t(j,i)=d_t(j,i)-vp(j)
4176 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4177 .and.(mnum.ne.5)) then
4178 ! if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
4179 ! if(itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
4181 call vecpr(vrot(1),dc(1,inres),vp)
4183 d_t(j,inres)=d_t(j,inres)-vp(j)
4188 ! write(iout,*) "The angular momentum after adjustment:"
4189 ! write(iout,*) (L(j),j=1,3)
4192 end subroutine inertia_tensor
4193 !-----------------------------------------------------------------------------
4194 subroutine angmom(cm,L)
4197 ! implicit real*8 (a-h,o-z)
4198 ! include 'DIMENSIONS'
4199 ! include 'COMMON.CONTROL'
4200 ! include 'COMMON.VAR'
4201 ! include 'COMMON.MD'
4202 ! include 'COMMON.CHAIN'
4203 ! include 'COMMON.DERIV'
4204 ! include 'COMMON.GEO'
4205 ! include 'COMMON.LOCAL'
4206 ! include 'COMMON.INTERACT'
4207 ! include 'COMMON.IOUNITS'
4208 ! include 'COMMON.NAMES'
4209 real(kind=8) :: mscab
4210 real(kind=8),dimension(3) :: L,cm,pr,vp,vrot,incr,v,pp
4211 integer :: iti,inres,i,j,mnum
4212 ! Calculate the angular momentum
4221 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
4223 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
4226 v(j)=incr(j)+0.5d0*d_t(j,i)
4229 incr(j)=incr(j)+d_t(j,i)
4231 call vecpr(pr(1),v(1),vp)
4233 L(j)=L(j)+mp(mnum)*vp(j)
4237 pp(j)=0.5d0*d_t(j,i)
4239 call vecpr(pr(1),pp(1),vp)
4241 L(j)=L(j)+Ip(mnum)*vp(j)
4249 iti=iabs(itype(i,mnum))
4257 pr(j)=c(j,inres)-cm(j)
4259 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4260 .and.(mnum.ne.5)) then
4262 v(j)=incr(j)+d_t(j,inres)
4269 call vecpr(pr(1),v(1),vp)
4270 ! write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),&
4271 ! " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
4273 L(j)=L(j)+mscab*vp(j)
4275 ! write (iout,*) "L",(l(j),j=1,3)
4276 if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4277 .and.(mnum.ne.5)) then
4279 v(j)=incr(j)+d_t(j,inres)
4281 call vecpr(dc(1,inres),d_t(1,inres),vp)
4283 L(j)=L(j)+Isc(iti,mnum)*vp(j)
4287 incr(j)=incr(j)+d_t(j,i)
4291 end subroutine angmom
4292 !-----------------------------------------------------------------------------
4293 subroutine vcm_vel(vcm)
4296 ! implicit real*8 (a-h,o-z)
4297 ! include 'DIMENSIONS'
4298 ! include 'COMMON.VAR'
4299 ! include 'COMMON.MD'
4300 ! include 'COMMON.CHAIN'
4301 ! include 'COMMON.DERIV'
4302 ! include 'COMMON.GEO'
4303 ! include 'COMMON.LOCAL'
4304 ! include 'COMMON.INTERACT'
4305 ! include 'COMMON.IOUNITS'
4306 real(kind=8),dimension(3) :: vcm,vv
4307 real(kind=8) :: summas,amas
4317 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
4319 summas=summas+mp(mnum)
4321 vcm(j)=vcm(j)+mp(mnum)*(vv(j)+0.5d0*d_t(j,i))
4325 amas=msc(iabs(itype(i,mnum)),mnum)
4330 if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4331 .and.(mnum.ne.5)) then
4333 vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres))
4337 vcm(j)=vcm(j)+amas*vv(j)
4341 vv(j)=vv(j)+d_t(j,i)
4344 ! write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas
4346 vcm(j)=vcm(j)/summas
4349 end subroutine vcm_vel
4350 !-----------------------------------------------------------------------------
4352 !-----------------------------------------------------------------------------
4354 ! RATTLE algorithm for velocity Verlet - step 1, UNRES
4358 ! implicit real*8 (a-h,o-z)
4359 ! include 'DIMENSIONS'
4361 ! include 'COMMON.CONTROL'
4362 ! include 'COMMON.VAR'
4363 ! include 'COMMON.MD'
4365 ! include 'COMMON.LANGEVIN'
4367 ! include 'COMMON.LANGEVIN.lang0'
4369 ! include 'COMMON.CHAIN'
4370 ! include 'COMMON.DERIV'
4371 ! include 'COMMON.GEO'
4372 ! include 'COMMON.LOCAL'
4373 ! include 'COMMON.INTERACT'
4374 ! include 'COMMON.IOUNITS'
4375 ! include 'COMMON.NAMES'
4376 ! include 'COMMON.TIME1'
4377 !el real(kind=8) :: gginv(2*nres,2*nres),&
4378 !el gdc(3,2*nres,2*nres)
4379 real(kind=8) :: dC_uncor(3,2*nres) !,&
4380 !el real(kind=8) :: Cmat(2*nres,2*nres)
4381 real(kind=8) :: x(2*nres),xcorr(3,2*nres) !maxres2=2*maxres
4382 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
4383 !el common /przechowalnia/ nbond
4384 integer :: max_rattle = 5
4385 logical :: lprn = .false., lprn1 = .false., not_done
4386 real(kind=8) :: tol_rattle = 1.0d-5
4388 integer :: ii,i,j,jj,l,ind,ind1,nres2
4391 !el /common/ przechowalnia
4393 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4394 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4395 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4397 if (lprn) write (iout,*) "RATTLE1"
4401 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4402 .and.(mnum.ne.5)) nbond=nbond+1
4404 ! Make a folded form of the Ginv-matrix
4417 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4422 if (itype(k,1).ne.10 .and. itype(k,mnum).ne.ntyp1_molec(mnum)&
4423 .and.(mnum.ne.5)) then
4427 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4435 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4446 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4450 if (itype(k,1).ne.10) then
4454 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4462 write (iout,*) "Matrix GGinv"
4463 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
4469 if (iter.gt.max_rattle) then
4470 write (iout,*) "Error - too many iterations in RATTLE."
4473 ! Calculate the matrix C = GG**(-1) dC_old o dC
4478 dC_uncor(j,ind1)=dC(j,i)
4482 if (itype(i,1).ne.10) then
4485 dC_uncor(j,ind1)=dC(j,i+nres)
4494 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
4498 if (itype(k,1).ne.10) then
4501 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
4506 ! Calculate deviations from standard virtual-bond lengths
4510 x(ind)=vbld(i+1)**2-vbl**2
4513 if (itype(i,1).ne.10) then
4515 x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4519 write (iout,*) "Coordinates and violations"
4521 write(iout,'(i5,3f10.5,5x,e15.5)') &
4522 i,(dC_uncor(j,i),j=1,3),x(i)
4524 write (iout,*) "Velocities and violations"
4528 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4529 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4533 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4534 .and.(mnum.ne.5)) then
4537 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4538 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4539 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4542 ! write (iout,*) "gdc"
4544 ! write (iout,*) "i",i
4546 ! write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4552 if (dabs(x(i)).gt.xmax) then
4556 if (xmax.lt.tol_rattle) then
4560 ! Calculate the matrix of the system of equations
4565 Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
4570 write (iout,*) "Matrix Cmat"
4571 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4573 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4574 ! Add constraint term to positions
4581 xx = xx+x(ii)*gdc(j,ind,ii)
4585 d_t_new(j,i)=d_t_new(j,i)-xx/d_time
4589 if (itype(i,1).ne.10) then
4594 xx = xx+x(ii)*gdc(j,ind,ii)
4597 dC(j,i+nres)=dC(j,i+nres)-xx
4598 d_t_new(j,i+nres)=d_t_new(j,i+nres)-xx/d_time
4602 ! Rebuild the chain using the new coordinates
4603 call chainbuild_cart
4605 write (iout,*) "New coordinates, Lagrange multipliers,",&
4606 " and differences between actual and standard bond lengths"
4610 xx=vbld(i+1)**2-vbl**2
4611 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4612 i,(dC(j,i),j=1,3),x(ind),xx
4616 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4619 xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4620 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4621 i,(dC(j,i+nres),j=1,3),x(ind),xx
4624 write (iout,*) "Velocities and violations"
4628 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4629 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4632 if (itype(i,1).ne.10) then
4634 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4635 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4636 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4643 10 write (iout,*) "Error - singularity in solving the system",&
4644 " of equations for Lagrange multipliers."
4648 "RATTLE inactive; use -DRATTLE switch at compile time."
4651 end subroutine rattle1
4652 !-----------------------------------------------------------------------------
4654 ! RATTLE algorithm for velocity Verlet - step 2, UNRES
4658 ! implicit real*8 (a-h,o-z)
4659 ! include 'DIMENSIONS'
4661 ! include 'COMMON.CONTROL'
4662 ! include 'COMMON.VAR'
4663 ! include 'COMMON.MD'
4665 ! include 'COMMON.LANGEVIN'
4667 ! include 'COMMON.LANGEVIN.lang0'
4669 ! include 'COMMON.CHAIN'
4670 ! include 'COMMON.DERIV'
4671 ! include 'COMMON.GEO'
4672 ! include 'COMMON.LOCAL'
4673 ! include 'COMMON.INTERACT'
4674 ! include 'COMMON.IOUNITS'
4675 ! include 'COMMON.NAMES'
4676 ! include 'COMMON.TIME1'
4677 !el real(kind=8) :: gginv(2*nres,2*nres),&
4678 !el gdc(3,2*nres,2*nres)
4679 real(kind=8) :: dC_uncor(3,2*nres) !,&
4680 !el Cmat(2*nres,2*nres)
4681 real(kind=8) :: x(2*nres) !maxres2=2*maxres
4682 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
4683 !el common /przechowalnia/ nbond
4684 integer :: max_rattle = 5
4685 logical :: lprn = .false., lprn1 = .false., not_done
4686 real(kind=8) :: tol_rattle = 1.0d-5
4690 !el /common/ przechowalnia
4691 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4692 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4693 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4695 if (lprn) write (iout,*) "RATTLE2"
4696 if (lprn) write (iout,*) "Velocity correction"
4697 ! Calculate the matrix G dC
4703 gdc(j,i,ind)=GGinv(i,ind)*dC(j,k)
4708 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4709 .and.(mnum.ne.5)) then
4712 gdc(j,i,ind)=GGinv(i,ind)*dC(j,k+nres)
4718 ! write (iout,*) "gdc"
4720 ! write (iout,*) "i",i
4722 ! write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4726 ! Calculate the matrix of the system of equations
4733 Cmat(ind,j)=Cmat(ind,j)+dC(k,i)*gdc(k,ind,j)
4739 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4740 .and.(mnum.ne.5)) then
4745 Cmat(ind,j)=Cmat(ind,j)+dC(k,i+nres)*gdc(k,ind,j)
4750 ! Calculate the scalar product dC o d_t_new
4754 x(ind)=scalar(d_t(1,i),dC(1,i))
4758 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4759 .and.(mnum.ne.5)) then
4761 x(ind)=scalar(d_t(1,i+nres),dC(1,i+nres))
4765 write (iout,*) "Velocities and violations"
4769 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4770 i,ind,(d_t(j,i),j=1,3),x(ind)
4774 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4775 .and.(mnum.ne.5)) then
4777 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4778 i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind)
4784 if (dabs(x(i)).gt.xmax) then
4788 if (xmax.lt.tol_rattle) then
4793 write (iout,*) "Matrix Cmat"
4794 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4796 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4797 ! Add constraint term to velocities
4804 xx = xx+x(ii)*gdc(j,ind,ii)
4806 d_t(j,i)=d_t(j,i)-xx
4811 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4812 .and.(mnum.ne.5)) then
4817 xx = xx+x(ii)*gdc(j,ind,ii)
4819 d_t(j,i+nres)=d_t(j,i+nres)-xx
4825 "New velocities, Lagrange multipliers violations"
4829 if (lprn) write (iout,'(2i5,3f10.5,5x,2e15.5)') &
4830 i,ind,(d_t(j,i),j=1,3),x(ind),scalar(d_t(1,i),dC(1,i))
4834 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4837 write (iout,'(2i5,3f10.5,5x,2e15.5)') &
4838 i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind),&
4839 scalar(d_t(1,i+nres),dC(1,i+nres))
4845 10 write (iout,*) "Error - singularity in solving the system",&
4846 " of equations for Lagrange multipliers."
4850 "RATTLE inactive; use -DRATTLE option at compile time."
4853 end subroutine rattle2
4854 !-----------------------------------------------------------------------------
4855 subroutine rattle_brown
4856 ! RATTLE/LINCS algorithm for Brownian dynamics, UNRES
4860 ! implicit real*8 (a-h,o-z)
4861 ! include 'DIMENSIONS'
4863 ! include 'COMMON.CONTROL'
4864 ! include 'COMMON.VAR'
4865 ! include 'COMMON.MD'
4867 ! include 'COMMON.LANGEVIN'
4869 ! include 'COMMON.LANGEVIN.lang0'
4871 ! include 'COMMON.CHAIN'
4872 ! include 'COMMON.DERIV'
4873 ! include 'COMMON.GEO'
4874 ! include 'COMMON.LOCAL'
4875 ! include 'COMMON.INTERACT'
4876 ! include 'COMMON.IOUNITS'
4877 ! include 'COMMON.NAMES'
4878 ! include 'COMMON.TIME1'
4879 !el real(kind=8) :: gginv(2*nres,2*nres),&
4880 !el gdc(3,2*nres,2*nres)
4881 real(kind=8) :: dC_uncor(3,2*nres) !,&
4882 !el real(kind=8) :: Cmat(2*nres,2*nres)
4883 real(kind=8) :: x(2*nres) !maxres2=2*maxres
4884 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
4885 !el common /przechowalnia/ nbond
4886 integer :: max_rattle = 5
4887 logical :: lprn = .false., lprn1 = .false., not_done
4888 real(kind=8) :: tol_rattle = 1.0d-5
4892 !el /common/ przechowalnia
4893 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4894 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4895 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4898 if (lprn) write (iout,*) "RATTLE_BROWN"
4901 if (itype(i,1).ne.10) nbond=nbond+1
4903 ! Make a folded form of the Ginv-matrix
4916 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4920 if (itype(k,1).ne.10) then
4924 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4931 if (itype(i,1).ne.10) then
4941 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4945 if (itype(k,1).ne.10) then
4949 if (j.eq.1 .and. l.eq.1)GGinv(ii,jj)=fricmat(ind,ind1)
4957 write (iout,*) "Matrix GGinv"
4958 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
4964 if (iter.gt.max_rattle) then
4965 write (iout,*) "Error - too many iterations in RATTLE."
4968 ! Calculate the matrix C = GG**(-1) dC_old o dC
4973 dC_uncor(j,ind1)=dC(j,i)
4977 if (itype(i,1).ne.10) then
4980 dC_uncor(j,ind1)=dC(j,i+nres)
4989 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
4993 if (itype(k,1).ne.10) then
4996 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
5001 ! Calculate deviations from standard virtual-bond lengths
5005 x(ind)=vbld(i+1)**2-vbl**2
5008 if (itype(i,1).ne.10) then
5010 x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
5014 write (iout,*) "Coordinates and violations"
5016 write(iout,'(i5,3f10.5,5x,e15.5)') &
5017 i,(dC_uncor(j,i),j=1,3),x(i)
5019 write (iout,*) "Velocities and violations"
5023 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5024 i,ind,(d_t(j,i),j=1,3),scalar(d_t(1,i),dC_old(1,i))
5027 if (itype(i,1).ne.10) then
5029 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5030 i+nres,ind,(d_t(j,i+nres),j=1,3),&
5031 scalar(d_t(1,i+nres),dC_old(1,i+nres))
5034 write (iout,*) "gdc"
5036 write (iout,*) "i",i
5038 write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
5044 if (dabs(x(i)).gt.xmax) then
5048 if (xmax.lt.tol_rattle) then
5052 ! Calculate the matrix of the system of equations
5057 Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
5062 write (iout,*) "Matrix Cmat"
5063 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
5065 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
5066 ! Add constraint term to positions
5073 xx = xx+x(ii)*gdc(j,ind,ii)
5076 d_t(j,i)=d_t(j,i)+xx/d_time
5081 if (itype(i,1).ne.10) then
5086 xx = xx+x(ii)*gdc(j,ind,ii)
5089 d_t(j,i+nres)=d_t(j,i+nres)+xx/d_time
5090 dC(j,i+nres)=dC(j,i+nres)+xx
5094 ! Rebuild the chain using the new coordinates
5095 call chainbuild_cart
5097 write (iout,*) "New coordinates, Lagrange multipliers,",&
5098 " and differences between actual and standard bond lengths"
5102 xx=vbld(i+1)**2-vbl**2
5103 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
5104 i,(dC(j,i),j=1,3),x(ind),xx
5107 if (itype(i,1).ne.10) then
5109 xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
5110 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
5111 i,(dC(j,i+nres),j=1,3),x(ind),xx
5114 write (iout,*) "Velocities and violations"
5118 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5119 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
5122 if (itype(i,1).ne.10) then
5124 write (iout,'(2i5,3f10.5,5x,e15.5)') &
5125 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
5126 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
5133 10 write (iout,*) "Error - singularity in solving the system",&
5134 " of equations for Lagrange multipliers."
5138 "RATTLE inactive; use -DRATTLE option at compile time"
5141 end subroutine rattle_brown
5142 !-----------------------------------------------------------------------------
5144 !-----------------------------------------------------------------------------
5145 subroutine friction_force
5150 ! implicit real*8 (a-h,o-z)
5151 ! include 'DIMENSIONS'
5152 ! include 'COMMON.VAR'
5153 ! include 'COMMON.CHAIN'
5154 ! include 'COMMON.DERIV'
5155 ! include 'COMMON.GEO'
5156 ! include 'COMMON.LOCAL'
5157 ! include 'COMMON.INTERACT'
5158 ! include 'COMMON.MD'
5160 ! include 'COMMON.LANGEVIN'
5162 ! include 'COMMON.LANGEVIN.lang0'
5164 ! include 'COMMON.IOUNITS'
5165 !el real(kind=8),dimension(6*nres) :: gamvec !(MAXRES6) maxres6=6*maxres
5166 !el common /syfek/ gamvec
5167 real(kind=8) :: vv(3),vvtot(3,nres),v_work(6*nres) !,&
5168 !el ginvfric(2*nres,2*nres) !maxres2=2*maxres
5169 !el common /przechowalnia/ ginvfric
5171 logical :: lprn = .false., checkmode = .false.
5172 integer :: i,j,ind,k,nres2,nres6,mnum
5176 if(.not.allocated(gamvec)) allocate(gamvec(nres6)) !(MAXRES6)
5177 if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
5185 d_t_work(j)=d_t(j,0)
5190 d_t_work(ind+j)=d_t(j,i)
5196 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5197 .and.(mnum.ne.5)) then
5199 d_t_work(ind+j)=d_t(j,i+nres)
5205 call fricmat_mult(d_t_work,fric_work)
5207 if (.not.checkmode) return
5210 write (iout,*) "d_t_work and fric_work"
5212 write (iout,'(i3,2e15.5)') i,d_t_work(i),fric_work(i)
5216 friction(j,0)=fric_work(j)
5221 friction(j,i)=fric_work(ind+j)
5227 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5228 .and.(mnum.ne.5)) then
5229 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
5231 friction(j,i+nres)=fric_work(ind+j)
5237 write(iout,*) "Friction backbone"
5239 write(iout,'(i5,3e15.5,5x,3e15.5)') &
5240 i,(friction(j,i),j=1,3),(d_t(j,i),j=1,3)
5242 write(iout,*) "Friction side chain"
5244 write(iout,'(i5,3e15.5,5x,3e15.5)') &
5245 i,(friction(j,i+nres),j=1,3),(d_t(j,i+nres),j=1,3)
5254 vvtot(j,i)=vv(j)+0.5d0*d_t(j,i)
5255 vvtot(j,i+nres)=vv(j)+d_t(j,i+nres)
5256 vv(j)=vv(j)+d_t(j,i)
5259 write (iout,*) "vvtot backbone and sidechain"
5261 write (iout,'(i5,3e15.5,5x,3e15.5)') i,(vvtot(j,i),j=1,3),&
5262 (vvtot(j,i+nres),j=1,3)
5267 v_work(ind+j)=vvtot(j,i)
5273 v_work(ind+j)=vvtot(j,i+nres)
5277 write (iout,*) "v_work gamvec and site-based friction forces"
5279 write (iout,'(i5,3e15.5)') i,v_work(i),gamvec(i),&
5283 ! fric_work1(i)=0.0d0
5285 ! fric_work1(i)=fric_work1(i)-A(j,i)*gamvec(j)*v_work(j)
5288 ! write (iout,*) "fric_work and fric_work1"
5290 ! write (iout,'(i5,2e15.5)') i,fric_work(i),fric_work1(i)
5296 ginvfric(i,j)=ginvfric(i,j)+ginv(i,k)*fricmat(k,j)
5300 write (iout,*) "ginvfric"
5302 write (iout,'(i5,100f8.3)') i,(ginvfric(i,j),j=1,dimen)
5304 write (iout,*) "symmetry check"
5307 write (iout,*) i,j,ginvfric(i,j)-ginvfric(j,i)
5312 end subroutine friction_force
5313 !-----------------------------------------------------------------------------
5314 subroutine setup_fricmat
5318 use control_data, only:time_Bcast
5319 use control, only:tcpu
5321 ! implicit real*8 (a-h,o-z)
5325 real(kind=8) :: time00
5327 ! include 'DIMENSIONS'
5328 ! include 'COMMON.VAR'
5329 ! include 'COMMON.CHAIN'
5330 ! include 'COMMON.DERIV'
5331 ! include 'COMMON.GEO'
5332 ! include 'COMMON.LOCAL'
5333 ! include 'COMMON.INTERACT'
5334 ! include 'COMMON.MD'
5335 ! include 'COMMON.SETUP'
5336 ! include 'COMMON.TIME1'
5337 ! integer licznik /0/
5340 ! include 'COMMON.LANGEVIN'
5342 ! include 'COMMON.LANGEVIN.lang0'
5344 ! include 'COMMON.IOUNITS'
5346 integer :: i,j,ind,ind1,m
5347 logical :: lprn = .false.
5348 real(kind=8) :: dtdi !el ,gamvec(2*nres)
5349 !el real(kind=8),dimension(2*nres,2*nres) :: ginvfric,fcopy
5350 ! real(kind=8),allocatable,dimension(:,:) :: fcopy
5351 !el real(kind=8),dimension(2*nres*(2*nres+1)/2) :: Ghalf !(mmaxres2) (mmaxres2=(maxres2*(maxres2+1)/2))
5352 !el common /syfek/ gamvec
5353 real(kind=8) :: work(8*2*nres)
5354 integer :: iwork(2*nres)
5355 !el common /przechowalnia/ ginvfric,Ghalf,fcopy
5356 integer :: ii,iti,k,l,nzero,nres2,nres6,ierr,mnum
5360 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
5361 if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
5362 if (fg_rank.ne.king) goto 10
5367 if(.not.allocated(gamvec)) allocate(gamvec(nres2)) !(MAXRES2)
5368 if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
5369 if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
5370 !el allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
5371 if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !maxres2=2*maxres
5373 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
5374 ! Zeroing out fricmat
5380 ! Load the friction coefficients corresponding to peptide groups
5385 gamvec(ind1)=gamp(mnum)
5387 ! Load the friction coefficients corresponding to side chains
5391 gamsc(ntyp1_molec(j),j)=1.0d0
5398 gamvec(ii)=gamsc(iabs(iti),mnum)
5400 if (surfarea) call sdarea(gamvec)
5402 ! write (iout,*) "Matrix A and vector gamma"
5404 ! write (iout,'(i2,$)') i
5406 ! write (iout,'(f4.1,$)') A(i,j)
5408 ! write (iout,'(f8.3)') gamvec(i)
5412 write (iout,*) "Vector gamvec"
5414 write (iout,'(i5,f10.5)') i, gamvec(i)
5418 ! The friction matrix
5423 dtdi=dtdi+A(j,k)*A(j,i)*gamvec(j)
5430 write (iout,'(//a)') "Matrix fricmat"
5431 call matout2(dimen,dimen,nres2,nres2,fricmat)
5433 if (lang.eq.2 .or. lang.eq.3) then
5434 ! Mass-scale the friction matrix if non-direct integration will be performed
5440 Ginvfric(i,j)=Ginvfric(i,j)+ &
5441 Gsqrm(i,k)*Gsqrm(l,j)*fricmat(k,l)
5446 ! Diagonalize the friction matrix
5451 Ghalf(ind)=Ginvfric(i,j)
5454 call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
5457 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
5458 " mass-scaled friction matrix"
5459 call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
5461 ! Precompute matrices for tinker stochastic integrator
5468 mt1(i,j)=mt1(i,j)+fricvec(k,i)*gsqrm(k,j)
5469 mt2(i,j)=mt2(i,j)+fricvec(k,i)*gsqrp(k,j)
5475 else if (lang.eq.4) then
5476 ! Diagonalize the friction matrix
5481 Ghalf(ind)=fricmat(i,j)
5484 call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
5487 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
5489 call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
5491 ! Determine the number of zero eigenvalues of the friction matrix
5492 nzero=max0(dimen-dimen1,0)
5493 ! do while (fricgam(nzero+1).le.1.0d-5 .and. nzero.lt.dimen)
5496 write (iout,*) "Number of zero eigenvalues:",nzero
5501 fricmat(i,j)=fricmat(i,j) &
5502 +fricvec(i,k)*fricvec(j,k)/fricgam(k)
5507 write (iout,'(//a)') "Generalized inverse of fricmat"
5508 call matout(dimen,dimen,nres6,nres6,fricmat)
5513 if (nfgtasks.gt.1) then
5514 if (fg_rank.eq.0) then
5515 ! The matching BROADCAST for fg processors is called in ERGASTULUM
5521 call MPI_Bcast(10,1,MPI_INTEGER,king,FG_COMM,IERROR)
5523 time_Bcast=time_Bcast+MPI_Wtime()-time00
5525 time_Bcast=time_Bcast+tcpu()-time00
5527 ! print *,"Processor",myrank,
5528 ! & " BROADCAST iorder in SETUP_FRICMAT"
5531 write (iout,*) "setup_fricmat licznik"!,licznik !sp
5537 ! Scatter the friction matrix
5538 call MPI_Scatterv(fricmat(1,1),nginv_counts(0),&
5539 nginv_start(0),MPI_DOUBLE_PRECISION,fcopy(1,1),&
5540 myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
5543 time_scatter=time_scatter+MPI_Wtime()-time00
5544 time_scatter_fmat=time_scatter_fmat+MPI_Wtime()-time00
5546 time_scatter=time_scatter+tcpu()-time00
5547 time_scatter_fmat=time_scatter_fmat+tcpu()-time00
5551 do j=1,2*my_ng_count
5552 fricmat(j,i)=fcopy(i,j)
5555 ! write (iout,*) "My chunk of fricmat"
5556 ! call MATOUT2(my_ng_count,dimen,maxres2,maxres2,fcopy)
5560 end subroutine setup_fricmat
5561 !-----------------------------------------------------------------------------
5562 subroutine stochastic_force(stochforcvec)
5565 use random, only:anorm_distr
5566 ! implicit real*8 (a-h,o-z)
5567 ! include 'DIMENSIONS'
5568 use control, only: tcpu
5573 ! include 'COMMON.VAR'
5574 ! include 'COMMON.CHAIN'
5575 ! include 'COMMON.DERIV'
5576 ! include 'COMMON.GEO'
5577 ! include 'COMMON.LOCAL'
5578 ! include 'COMMON.INTERACT'
5579 ! include 'COMMON.MD'
5580 ! include 'COMMON.TIME1'
5582 ! include 'COMMON.LANGEVIN'
5584 ! include 'COMMON.LANGEVIN.lang0'
5586 ! include 'COMMON.IOUNITS'
5588 real(kind=8) :: x,sig,lowb,highb
5589 real(kind=8) :: ff(3),force(3,0:2*nres),zeta2,lowb2
5590 real(kind=8) :: highb2,sig2,forcvec(6*nres),stochforcvec(6*nres)
5591 real(kind=8) :: time00
5592 logical :: lprn = .false.
5593 integer :: i,j,ind,mnum
5597 stochforc(j,i)=0.0d0
5607 ! Compute the stochastic forces acting on bodies. Store in force.
5613 force(j,i)=anorm_distr(x,sig,lowb,highb)
5621 force(j,i+nres)=anorm_distr(x,sig2,lowb2,highb2)
5625 time_fsample=time_fsample+MPI_Wtime()-time00
5627 time_fsample=time_fsample+tcpu()-time00
5629 ! Compute the stochastic forces acting on virtual-bond vectors.
5635 stochforc(j,i)=ff(j)+0.5d0*force(j,i)
5638 ff(j)=ff(j)+force(j,i)
5640 ! if (itype(i+1,1).ne.ntyp1) then
5642 if (itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
5644 stochforc(j,i)=stochforc(j,i)+force(j,i+nres+1)
5645 ff(j)=ff(j)+force(j,i+nres+1)
5650 stochforc(j,0)=ff(j)+force(j,nnt+nres)
5654 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5655 .and.(mnum.ne.5)) then
5656 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
5658 stochforc(j,i+nres)=force(j,i+nres)
5664 stochforcvec(j)=stochforc(j,0)
5669 stochforcvec(ind+j)=stochforc(j,i)
5675 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5676 .and.(mnum.ne.5)) then
5677 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
5679 stochforcvec(ind+j)=stochforc(j,i+nres)
5685 write (iout,*) "stochforcvec"
5687 write(iout,'(i5,e15.5)') i,stochforcvec(i)
5689 write(iout,*) "Stochastic forces backbone"
5691 write(iout,'(i5,3e15.5)') i,(stochforc(j,i),j=1,3)
5693 write(iout,*) "Stochastic forces side chain"
5695 write(iout,'(i5,3e15.5)') &
5696 i,(stochforc(j,i+nres),j=1,3)
5704 write (iout,*) i,ind
5706 forcvec(ind+j)=force(j,i)
5711 write (iout,*) i,ind
5713 forcvec(j+ind)=force(j,i+nres)
5718 write (iout,*) "forcvec"
5722 write (iout,'(2i3,2f10.5)') i,j,force(j,i),&
5729 write (iout,'(2i3,2f10.5)') i,j,force(j,i+nres),&
5738 end subroutine stochastic_force
5739 !-----------------------------------------------------------------------------
5740 subroutine sdarea(gamvec)
5742 ! Scale the friction coefficients according to solvent accessible surface areas
5743 ! Code adapted from TINKER
5747 ! implicit real*8 (a-h,o-z)
5748 ! include 'DIMENSIONS'
5749 ! include 'COMMON.CONTROL'
5750 ! include 'COMMON.VAR'
5751 ! include 'COMMON.MD'
5753 ! include 'COMMON.LANGEVIN'
5755 ! include 'COMMON.LANGEVIN.lang0'
5757 ! include 'COMMON.CHAIN'
5758 ! include 'COMMON.DERIV'
5759 ! include 'COMMON.GEO'
5760 ! include 'COMMON.LOCAL'
5761 ! include 'COMMON.INTERACT'
5762 ! include 'COMMON.IOUNITS'
5763 ! include 'COMMON.NAMES'
5764 real(kind=8),dimension(2*nres) :: radius,gamvec !(maxres2)
5765 real(kind=8),parameter :: twosix = 1.122462048309372981d0
5766 logical :: lprn = .false.
5767 real(kind=8) :: probe,area,ratio
5768 integer :: i,j,ind,iti,mnum
5770 ! determine new friction coefficients every few SD steps
5772 ! set the atomic radii to estimates of sigma values
5774 ! print *,"Entered sdarea"
5780 ! Load peptide group radii
5783 radius(i)=pstok(mnum)
5785 ! Load side chain radii
5789 radius(i+nres)=restok(iti,mnum)
5792 ! write (iout,*) "i",i," radius",radius(i)
5795 radius(i) = radius(i) / twosix
5796 if (radius(i) .ne. 0.0d0) radius(i) = radius(i) + probe
5799 ! scale atomic friction coefficients by accessible area
5801 if (lprn) write (iout,*) &
5802 "Original gammas, surface areas, scaling factors, new gammas, ",&
5803 "std's of stochastic forces"
5806 if (radius(i).gt.0.0d0) then
5807 call surfatom (i,area,radius)
5808 ratio = dmax1(area/(4.0d0*pi*radius(i)**2),1.0d-1)
5809 if (lprn) write (iout,'(i5,3f10.5,$)') &
5810 i,gamvec(ind+1),area,ratio
5813 gamvec(ind) = ratio * gamvec(ind)
5816 stdforcp(i)=stdfp(mnum)*dsqrt(gamvec(ind))
5817 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcp(i)
5821 if (radius(i+nres).gt.0.0d0) then
5822 call surfatom (i+nres,area,radius)
5823 ratio = dmax1(area/(4.0d0*pi*radius(i+nres)**2),1.0d-1)
5824 if (lprn) write (iout,'(i5,3f10.5,$)') &
5825 i,gamvec(ind+1),area,ratio
5828 gamvec(ind) = ratio * gamvec(ind)
5831 stdforcsc(i)=stdfsc(itype(i,mnum),mnum)*dsqrt(gamvec(ind))
5832 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcsc(i)
5837 end subroutine sdarea
5838 !-----------------------------------------------------------------------------
5840 !-----------------------------------------------------------------------------
5843 ! ###################################################
5844 ! ## COPYRIGHT (C) 1996 by Jay William Ponder ##
5845 ! ## All Rights Reserved ##
5846 ! ###################################################
5848 ! ################################################################
5850 ! ## subroutine surfatom -- exposed surface area of an atom ##
5852 ! ################################################################
5855 ! "surfatom" performs an analytical computation of the surface
5856 ! area of a specified atom; a simplified version of "surface"
5858 ! literature references:
5860 ! T. J. Richmond, "Solvent Accessible Surface Area and
5861 ! Excluded Volume in Proteins", Journal of Molecular Biology,
5864 ! L. Wesson and D. Eisenberg, "Atomic Solvation Parameters
5865 ! Applied to Molecular Dynamics of Proteins in Solution",
5866 ! Protein Science, 1, 227-235 (1992)
5868 ! variables and parameters:
5870 ! ir number of atom for which area is desired
5871 ! area accessible surface area of the atom
5872 ! radius radii of each of the individual atoms
5875 subroutine surfatom(ir,area,radius)
5877 ! implicit real*8 (a-h,o-z)
5878 ! include 'DIMENSIONS'
5880 ! include 'COMMON.GEO'
5881 ! include 'COMMON.IOUNITS'
5883 integer :: nsup,nstart_sup
5884 ! double precision c,dc,dc_old,d_c_work,xloc,xrot,dc_norm
5885 ! common /chain/ c(3,maxres2+2),dc(3,0:maxres2),dc_old(3,0:maxres2),
5886 ! & xloc(3,maxres),xrot(3,maxres),dc_norm(3,0:maxres2),
5887 ! & dc_work(MAXRES6),nres,nres0
5888 integer,parameter :: maxarc=300
5892 integer :: mi,ni,narc
5893 integer :: key(maxarc)
5894 integer :: intag(maxarc)
5895 integer :: intag1(maxarc)
5896 real(kind=8) :: area,arcsum
5897 real(kind=8) :: arclen,exang
5898 real(kind=8) :: delta,delta2
5899 real(kind=8) :: eps,rmove
5900 real(kind=8) :: xr,yr,zr
5901 real(kind=8) :: rr,rrsq
5902 real(kind=8) :: rplus,rminus
5903 real(kind=8) :: axx,axy,axz
5904 real(kind=8) :: ayx,ayy
5905 real(kind=8) :: azx,azy,azz
5906 real(kind=8) :: uxj,uyj,uzj
5907 real(kind=8) :: tx,ty,tz
5908 real(kind=8) :: txb,tyb,td
5909 real(kind=8) :: tr2,tr,txr,tyr
5910 real(kind=8) :: tk1,tk2
5911 real(kind=8) :: thec,the,t,tb
5912 real(kind=8) :: txk,tyk,tzk
5913 real(kind=8) :: t1,ti,tf,tt
5914 real(kind=8) :: txj,tyj,tzj
5915 real(kind=8) :: ccsq,cc,xysq
5916 real(kind=8) :: bsqk,bk,cosine
5917 real(kind=8) :: dsqj,gi,pix2
5918 real(kind=8) :: therk,dk,gk
5919 real(kind=8) :: risqk,rik
5920 real(kind=8) :: radius(2*nres) !(maxatm) (maxatm=maxres2)
5921 real(kind=8) :: ri(maxarc),risq(maxarc)
5922 real(kind=8) :: ux(maxarc),uy(maxarc),uz(maxarc)
5923 real(kind=8) :: xc(maxarc),yc(maxarc),zc(maxarc)
5924 real(kind=8) :: xc1(maxarc),yc1(maxarc),zc1(maxarc)
5925 real(kind=8) :: dsq(maxarc),bsq(maxarc)
5926 real(kind=8) :: dsq1(maxarc),bsq1(maxarc)
5927 real(kind=8) :: arci(maxarc),arcf(maxarc)
5928 real(kind=8) :: ex(maxarc),lt(maxarc),gr(maxarc)
5929 real(kind=8) :: b(maxarc),b1(maxarc),bg(maxarc)
5930 real(kind=8) :: kent(maxarc),kout(maxarc)
5931 real(kind=8) :: ther(maxarc)
5932 logical :: moved,top
5933 logical :: omit(maxarc)
5936 ! maxatm = 2*nres !maxres2 maxres2=2*maxres
5937 ! maxlight = 8*maxatm
5940 ! maxtors = 4*maxatm
5942 ! zero out the surface area for the sphere of interest
5945 ! write (2,*) "ir",ir," radius",radius(ir)
5946 if (radius(ir) .eq. 0.0d0) return
5948 ! set the overlap significance and connectivity shift
5952 delta2 = delta * delta
5957 ! store coordinates and radius of the sphere of interest
5965 ! initialize values of some counters and summations
5974 ! test each sphere to see if it overlaps the sphere of interest
5977 if (i.eq.ir .or. radius(i).eq.0.0d0) goto 30
5978 rplus = rr + radius(i)
5980 if (abs(tx) .ge. rplus) goto 30
5982 if (abs(ty) .ge. rplus) goto 30
5984 if (abs(tz) .ge. rplus) goto 30
5986 ! check for sphere overlap by testing distance against radii
5988 xysq = tx*tx + ty*ty
5989 if (xysq .lt. delta2) then
5996 if (rplus-cc .le. delta) goto 30
5997 rminus = rr - radius(i)
5999 ! check to see if sphere of interest is completely buried
6001 if (cc-abs(rminus) .le. delta) then
6002 if (rminus .le. 0.0d0) goto 170
6006 ! check for too many overlaps with sphere of interest
6008 if (io .ge. maxarc) then
6010 20 format (/,' SURFATOM -- Increase the Value of MAXARC')
6014 ! get overlap between current sphere and sphere of interest
6023 gr(io) = (ccsq+rplus*rminus) / (2.0d0*rr*b1(io))
6029 ! case where no other spheres overlap the sphere of interest
6032 area = 4.0d0 * pi * rrsq
6036 ! case where only one sphere overlaps the sphere of interest
6039 area = pix2 * (1.0d0 + gr(1))
6040 area = mod(area,4.0d0*pi) * rrsq
6044 ! case where many spheres intersect the sphere of interest;
6045 ! sort the intersecting spheres by their degree of overlap
6047 call sort2 (io,gr,key)
6050 intag(i) = intag1(k)
6059 ! get radius of each overlap circle on surface of the sphere
6064 risq(i) = rrsq - gi*gi
6065 ri(i) = sqrt(risq(i))
6066 ther(i) = 0.5d0*pi - asin(min(1.0d0,max(-1.0d0,gr(i))))
6069 ! find boundary of inaccessible area on sphere of interest
6072 if (.not. omit(k)) then
6079 ! check to see if J circle is intersecting K circle;
6080 ! get distance between circle centers and sum of radii
6083 if (omit(j)) goto 60
6084 cc = (txk*xc(j)+tyk*yc(j)+tzk*zc(j))/(bk*b(j))
6085 cc = acos(min(1.0d0,max(-1.0d0,cc)))
6086 td = therk + ther(j)
6088 ! check to see if circles enclose separate regions
6090 if (cc .ge. td) goto 60
6092 ! check for circle J completely inside circle K
6094 if (cc+ther(j) .lt. therk) goto 40
6096 ! check for circles that are essentially parallel
6098 if (cc .gt. delta) goto 50
6103 ! check to see if sphere of interest is completely buried
6106 if (pix2-cc .le. td) goto 170
6112 ! find T value of circle intersections
6115 if (omit(k)) goto 110
6130 ! rotation matrix elements
6142 if (.not. omit(j)) then
6147 ! rotate spheres so K vector colinear with z-axis
6149 uxj = txj*axx + tyj*axy - tzj*axz
6150 uyj = tyj*ayy - txj*ayx
6151 uzj = txj*azx + tyj*azy + tzj*azz
6152 cosine = min(1.0d0,max(-1.0d0,uzj/b(j)))
6153 if (acos(cosine) .lt. therk+ther(j)) then
6154 dsqj = uxj*uxj + uyj*uyj
6159 tr2 = risqk*dsqj - tb*tb
6165 ! get T values of intersection for K circle
6168 tb = min(1.0d0,max(-1.0d0,tb))
6170 if (tyb-txr .lt. 0.0d0) tk1 = pix2 - tk1
6172 tb = min(1.0d0,max(-1.0d0,tb))
6174 if (tyb+txr .lt. 0.0d0) tk2 = pix2 - tk2
6175 thec = (rrsq*uzj-gk*bg(j)) / (rik*ri(j)*b(j))
6176 if (abs(thec) .lt. 1.0d0) then
6178 else if (thec .ge. 1.0d0) then
6180 else if (thec .le. -1.0d0) then
6184 ! see if "tk1" is entry or exit point; check t=0 point;
6185 ! "ti" is exit point, "tf" is entry point
6187 cosine = min(1.0d0,max(-1.0d0, &
6188 (uzj*gk-uxj*rik)/(b(j)*rr)))
6189 if ((acos(cosine)-ther(j))*(tk2-tk1) .le. 0.0d0) then
6197 if (narc .ge. maxarc) then
6199 70 format (/,' SURFATOM -- Increase the Value',&
6203 if (tf .le. ti) then
6224 ! special case; K circle without intersections
6226 if (narc .le. 0) goto 90
6228 ! general case; sum up arclength and set connectivity code
6230 call sort2 (narc,arci,key)
6235 if (narc .gt. 1) then
6238 if (t .lt. arci(j)) then
6239 arcsum = arcsum + arci(j) - t
6240 exang = exang + ex(ni)
6242 if (jb .ge. maxarc) then
6244 80 format (/,' SURFATOM -- Increase the Value',&
6249 kent(jb) = maxarc*i + k
6251 kout(jb) = maxarc*k + i
6260 arcsum = arcsum + pix2 - t
6262 exang = exang + ex(ni)
6265 kent(jb) = maxarc*i + k
6267 kout(jb) = maxarc*k + i
6274 arclen = arclen + gr(k)*arcsum
6277 if (arclen .eq. 0.0d0) goto 170
6278 if (jb .eq. 0) goto 150
6280 ! find number of independent boundaries and check connectivity
6284 if (kout(k) .ne. 0) then
6291 if (m .eq. kent(ii)) then
6294 if (j .eq. jb) goto 150
6306 ! attempt to fix connectivity error by moving atom slightly
6310 140 format (/,' SURFATOM -- Connectivity Error at Atom',i6)
6319 ! compute the exposed surface area for the sphere of interest
6322 area = ib*pix2 + exang + arclen
6323 area = mod(area,4.0d0*pi) * rrsq
6325 ! attempt to fix negative area by moving atom slightly
6327 if (area .lt. 0.0d0) then
6330 160 format (/,' SURFATOM -- Negative Area at Atom',i6)
6341 end subroutine surfatom
6342 !----------------------------------------------------------------
6343 !----------------------------------------------------------------
6344 subroutine alloc_MD_arrays
6345 !EL Allocation of arrays used by MD module
6347 integer :: nres2,nres6
6350 !----------------------
6354 allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
6355 allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
6356 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
6357 allocate(fricvec(nres2,nres2))
6358 allocate(pfric_mat(nres2,nres2),vfric_mat(nres2,nres2))
6359 allocate(afric_mat(nres2,nres2),prand_mat(nres2,nres2))
6360 allocate(vrand_mat1(nres2,nres2),vrand_mat2(nres2,nres2)) !(MAXRES2,MAXRES2)
6361 allocate(pfric0_mat(nres2,nres2,0:maxflag_stoch))
6362 allocate(afric0_mat(nres2,nres2,0:maxflag_stoch))
6363 allocate(vfric0_mat(nres2,nres2,0:maxflag_stoch))
6364 allocate(prand0_mat(nres2,nres2,0:maxflag_stoch))
6365 allocate(vrand0_mat1(nres2,nres2,0:maxflag_stoch))
6366 allocate(vrand0_mat2(nres2,nres2,0:maxflag_stoch)) !(MAXRES2,MAXRES2,0:maxflag_stoch)
6367 allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
6369 allocate(mt1(nres2,nres2),mt2(nres2,nres2),mt3(nres2,nres2)) !(maxres2,maxres2)
6370 !----------------------
6372 ! commom.langevin.lang0
6374 allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
6375 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
6376 allocate(fricvec(nres2,nres2)) !(MAXRES2,MAXRES2)
6377 allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
6378 allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
6381 if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2))
6382 !----------------------
6383 ! commom.hairpin in CSA module
6384 !----------------------
6385 ! common.mce in MCM_MD module
6386 !----------------------
6388 ! common /mdgrad/ in module.energy
6389 ! common /back_constr/ in module.energy
6390 ! common /qmeas/ in module.energy
6393 allocate(potEcomp(0:n_ene+4)) !(0:n_ene+4)
6395 allocate(d_t(3,0:nres2),d_a(3,0:nres2),d_t_old(3,0:nres2)) !(3,0:MAXRES2)
6396 allocate(d_a_work(nres6)) !(6*MAXRES)
6398 allocate(DM(nres2),DU1(nres2),DU2(nres2))
6399 allocate(DMorig(nres2),DU1orig(nres2),DU2orig(nres2))
6401 write (iout,*) "Before A Allocation",nfgtasks-1
6403 allocate(Gmat(nres2,nres2),A(nres2,nres2))
6404 if(.not.allocated(Ginv)) allocate(Ginv(nres2,nres2)) !in control: ergastulum
6405 allocate(Gsqrp(nres2,nres2),Gsqrm(nres2,nres2),Gvec(nres2,nres2)) !(maxres2,maxres2)
6407 allocate(Geigen(nres2)) !(maxres2)
6408 if(.not.allocated(vtot)) allocate(vtot(nres2)) !(maxres2)
6409 ! common /inertia/ in io_conf: parmread
6410 ! real(kind=8),dimension(:),allocatable :: ISC,msc !(ntyp+1)
6411 ! common /langevin/in io read_MDpar
6412 ! real(kind=8),dimension(:),allocatable :: gamsc !(ntyp1)
6413 ! real(kind=8),dimension(:),allocatable :: stdfsc !(ntyp)
6414 ! in io_conf: parmread
6415 ! real(kind=8),dimension(:),allocatable :: restok !(ntyp+1)
6416 ! common /mdpmpi/ in control: ergastulum
6417 if(.not.allocated(ng_start)) allocate(ng_start(0:nfgtasks-1))
6418 if(.not.allocated(ng_counts)) allocate(ng_counts(0:nfgtasks-1))
6419 if(.not.allocated(nginv_counts)) allocate(nginv_counts(0:nfgtasks-1)) !(0:MaxProcs-1)
6420 if(.not.allocated(nginv_start)) allocate(nginv_start(0:nfgtasks)) !(0:MaxProcs)
6421 !----------------------
6422 ! common.muca in read_muca
6423 ! common /double_muca/
6424 ! real(kind=8) :: elow,ehigh,factor,hbin,factor_min
6425 ! real(kind=8),dimension(:),allocatable :: emuca,nemuca,&
6426 ! nemuca2,hist !(4*maxres)
6427 ! common /integer_muca/
6428 ! integer :: nmuca,imtime,muca_smooth
6430 ! real(kind=8),dimension(:),allocatable :: elowi,ehighi !(maxprocs)
6431 !----------------------
6433 ! common /mdgrad/ in module.energy
6434 ! common /back_constr/ in module.energy
6435 ! common /qmeas/ in module.energy
6439 allocate(d_t_work(nres6),d_t_work_new(nres6),d_af_work(nres6))
6440 allocate(d_as_work(nres6),kinetic_force(nres6)) !(MAXRES6)
6441 allocate(d_t_new(3,0:nres2),d_a_old(3,0:nres2),d_a_short(3,0:nres2)) !,d_a !(3,0:MAXRES2)
6442 allocate(stdforcp(nres),stdforcsc(nres)) !(MAXRES)
6443 !----------------------
6445 allocate(D_ban(nres6)) !(MAXRES6) maxres6=6*maxres
6446 ! common /stochcalc/ stochforcvec
6447 allocate(stochforcvec(nres6)) !(MAXRES6) maxres6=6*maxres
6448 !----------------------
6450 end subroutine alloc_MD_arrays
6451 !-----------------------------------------------------------------------------
6452 !-----------------------------------------------------------------------------