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 !el 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(6*nres) :: Td !(maxres6) maxres6=6*maxres
191 real(kind=8),dimension(2*nres) :: ppvec !(maxres2) maxres2=2*maxres
192 !el common /stochcalc/ stochforcvec
193 !el real(kind=8),dimension(3) :: cm !el
194 !el common /gucio/ cm
196 logical :: lprn = .false.,lprn1 = .false.
197 integer :: maxiter = 5
198 real(kind=8) :: difftol = 1.0d-5
199 real(kind=8) :: xx,diffmax,blen2,diffbond,tt0
200 integer :: i,j,nbond,k,ind,ind1,iter
201 integer :: nres2,nres6
206 if (.not.allocated(stochforcvec)) allocate(stochforcvec(nres6)) !(MAXRES6) maxres6=6*maxres
210 if (itype(i).ne.10) nbond=nbond+1
214 write (iout,*) "Generalized inverse of fricmat"
215 call matout(dimen,dimen,nres6,nres6,fricmat)
227 Bmat(ind+j,ind1)=dC_norm(j,i)
232 if (itype(i).ne.10) then
235 Bmat(ind+j,ind1)=dC_norm(j,i+nres)
241 write (iout,*) "Matrix Bmat"
242 call MATOUT(nbond,dimen,nres6,nres6,Bmat)
248 GBmat(i,j)=GBmat(i,j)+fricmat(i,k)*Bmat(k,j)
253 write (iout,*) "Matrix GBmat"
254 call MATOUT(nbond,dimen,nres6,nres2,Gbmat)
260 Cmat_(i,j)=Cmat_(i,j)+Bmat(k,i)*GBmat(k,j)
265 write (iout,*) "Matrix Cmat"
266 call MATOUT(nbond,nbond,nres2,nres2,Cmat_)
268 call matinvert(nbond,nres2,Cmat_,Cinv,osob)
270 write (iout,*) "Matrix Cinv"
271 call MATOUT(nbond,nbond,nres2,nres2,Cinv)
277 Tmat(i,j)=Tmat(i,j)+GBmat(i,k)*Cinv(k,j)
282 write (iout,*) "Matrix Tmat"
283 call MATOUT(nbond,dimen,nres6,nres2,Tmat)
293 Pmat(i,j)=Pmat(i,j)-Tmat(i,k)*Bmat(j,k)
298 write (iout,*) "Matrix Pmat"
299 call MATOUT(dimen,dimen,nres6,nres6,Pmat)
306 Td(i)=Td(i)+vbl*Tmat(i,ind)
309 if (itype(k).ne.10) then
311 Td(i)=Td(i)+vbldsc0(1,itype(k))*Tmat(i,ind)
316 write (iout,*) "Vector Td"
318 write (iout,'(i5,f10.5)') i,Td(i)
321 call stochastic_force(stochforcvec)
323 write (iout,*) "stochforcvec"
325 write (iout,*) i,stochforcvec(i)
329 zapas(j)=-gcart(j,0)+stochforcvec(j)
331 dC_work(j)=dC_old(j,0)
337 zapas(ind)=-gcart(j,i)+stochforcvec(ind)
338 dC_work(ind)=dC_old(j,i)
342 if (itype(i).ne.10) then
345 zapas(ind)=-gxcart(j,i)+stochforcvec(ind)
346 dC_work(ind)=dC_old(j,i+nres)
352 write (iout,*) "Initial d_t_work"
354 write (iout,*) i,d_t_work(i)
361 d_t_work(i)=d_t_work(i)+fricmat(i,j)*zapas(j)
368 zapas(i)=zapas(i)+Pmat(i,j)*(dC_work(j)+d_t_work(j)*d_time)
372 write (iout,*) "Final d_t_work and zapas"
374 write (iout,*) i,d_t_work(i),zapas(i)
388 dc_work(ind+j)=dc(j,i)
394 d_t(j,i+nres)=d_t_work(ind+j)
395 dc(j,i+nres)=zapas(ind+j)
396 dc_work(ind+j)=dc(j,i+nres)
402 write (iout,*) "Before correction for rotational lengthening"
403 write (iout,*) "New coordinates",&
404 " and differences between actual and standard bond lengths"
409 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
413 if (itype(i).ne.10) then
415 xx=vbld(i+nres)-vbldsc0(1,itype(i))
416 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
417 i,(dC(j,i+nres),j=1,3),xx
421 ! Second correction (rotational lengthening)
427 blen2 = scalar(dc(1,i),dc(1,i))
428 ppvec(ind)=2*vbl**2-blen2
429 diffbond=dabs(vbl-dsqrt(blen2))
430 if (diffbond.gt.diffmax) diffmax=diffbond
431 if (ppvec(ind).gt.0.0d0) then
432 ppvec(ind)=dsqrt(ppvec(ind))
437 write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind)
441 if (itype(i).ne.10) then
443 blen2 = scalar(dc(1,i+nres),dc(1,i+nres))
444 ppvec(ind)=2*vbldsc0(1,itype(i))**2-blen2
445 diffbond=dabs(vbldsc0(1,itype(i))-dsqrt(blen2))
446 if (diffbond.gt.diffmax) diffmax=diffbond
447 if (ppvec(ind).gt.0.0d0) then
448 ppvec(ind)=dsqrt(ppvec(ind))
453 write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind)
457 if (lprn) write (iout,*) "iter",iter," diffmax",diffmax
458 if (diffmax.lt.difftol) goto 10
462 Td(i)=Td(i)+ppvec(j)*Tmat(i,j)
468 zapas(i)=zapas(i)+Pmat(i,j)*dc_work(j)
479 dc_work(ind+j)=zapas(ind+j)
484 if (itype(i).ne.10) then
486 dc(j,i+nres)=zapas(ind+j)
487 dc_work(ind+j)=zapas(ind+j)
492 ! Building the chain from the newly calculated coordinates
495 if (large.and. mod(itime,ntwe).eq.0) then
496 write (iout,*) "Cartesian and internal coordinates: step 1"
499 write (iout,'(a)') "Potential forces"
501 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(-gcart(j,i),j=1,3),&
504 write (iout,'(a)') "Stochastic forces"
506 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(stochforc(j,i),j=1,3),&
507 (stochforc(j,i+nres),j=1,3)
509 write (iout,'(a)') "Velocities"
511 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
512 (d_t(j,i+nres),j=1,3)
517 write (iout,*) "After correction for rotational lengthening"
518 write (iout,*) "New coordinates",&
519 " and differences between actual and standard bond lengths"
524 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
528 if (itype(i).ne.10) then
530 xx=vbld(i+nres)-vbldsc0(1,itype(i))
531 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
532 i,(dC(j,i+nres),j=1,3),xx
537 ! write (iout,*) "Too many attempts at correcting the bonds"
545 ! Calculate energy and forces
547 call etotal(potEcomp)
548 potE=potEcomp(0)-potEcomp(20)
551 ! Calculate the kinetic and total energy and the kinetic temperature
554 t_enegrad=t_enegrad+MPI_Wtime()-tt0
556 t_enegrad=t_enegrad+tcpu()-tt0
559 kinetic_T=2.0d0/(dimen*Rb)*EK
561 end subroutine brown_step
562 !-----------------------------------------------------------------------------
564 !-----------------------------------------------------------------------------
565 subroutine gauss(RO,AP,MT,M,N,*)
567 ! CALCULATES (RO**(-1))*AP BY GAUSS ELIMINATION
568 ! RO IS A SQUARE MATRIX
569 ! THE CALCULATED PRODUCT IS STORED IN AP
570 ! ABNORMAL EXIT IF RO IS SINGULAR
572 integer :: MT, M, N, M1,I,J,IM,&
574 real(kind=8) :: RO(MT,M),AP(MT,N),X,RM,PR,Y
580 if(dabs(X).le.1.0D-13) return 1
592 if(DABS(RO(J,I)).LE.RM) goto 2
606 if(dabs(X).le.1.0E-13) return 1
615 8 AP(J,K)=AP(J,K)-Y*AP(I,K)
617 9 RO(J,K)=RO(J,K)-Y*RO(I,K)
621 if(dabs(X).le.1.0E-13) return 1
631 15 X=X-AP(K,J)*RO(MI,K)
636 !-----------------------------------------------------------------------------
638 !-----------------------------------------------------------------------------
639 subroutine kinetic(KE_total)
640 !----------------------------------------------------------------
641 ! This subroutine calculates the total kinetic energy of the chain
642 !-----------------------------------------------------------------
644 ! implicit real*8 (a-h,o-z)
645 ! include 'DIMENSIONS'
646 ! include 'COMMON.VAR'
647 ! include 'COMMON.CHAIN'
648 ! include 'COMMON.DERIV'
649 ! include 'COMMON.GEO'
650 ! include 'COMMON.LOCAL'
651 ! include 'COMMON.INTERACT'
652 ! include 'COMMON.MD'
653 ! include 'COMMON.IOUNITS'
654 real(kind=8) :: KE_total
657 real(kind=8) :: KEt_p,KEt_sc,KEr_p,KEr_sc,incr(3),&
660 write (iout,*) "Velocities, kietic"
662 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
663 (d_t(j,i+nres),j=1,3)
668 ! write (iout,*) "ISC",(isc(itype(i)),i=1,nres)
669 ! The translational part for peptide virtual bonds
674 ! write (iout,*) "Kinetic trp:",i,(incr(j),j=1,3)
676 v(j)=incr(j)+0.5d0*d_t(j,i)
678 vtot(i)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
679 KEt_p=KEt_p+(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
681 incr(j)=incr(j)+d_t(j,i)
684 ! write(iout,*) 'KEt_p', KEt_p
685 ! The translational part for the side chain virtual bond
686 ! Only now we can initialize incr with zeros. It must be equal
687 ! to the velocities of the first Calpha.
693 if (itype(i).eq.10) then
699 v(j)=incr(j)+d_t(j,nres+i)
702 ! write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3)
703 ! write (iout,*) "i",i," msc",msc(iti)," v",(v(j),j=1,3)
704 KEt_sc=KEt_sc+msc(iti)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
705 vtot(i+nres)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
707 incr(j)=incr(j)+d_t(j,i)
711 ! write(iout,*) 'KEt_sc', KEt_sc
712 ! The part due to stretching and rotation of the peptide groups
715 ! write (iout,*) "i",i
716 ! write (iout,*) "i",i," mag1",mag1," mag2",mag2
720 ! write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3)
721 KEr_p=KEr_p+(incr(1)*incr(1)+incr(2)*incr(2) &
725 ! write(iout,*) 'KEr_p', KEr_p
726 ! The rotational part of the side chain virtual bond
730 if (itype(i).ne.10) then
732 incr(j)=d_t(j,nres+i)
734 ! write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3)
735 KEr_sc=KEr_sc+Isc(iti)*(incr(1)*incr(1)+incr(2)*incr(2)+ &
739 ! The total kinetic energy
741 ! write(iout,*) 'KEr_sc', KEr_sc
742 KE_total=0.5d0*(mp*KEt_p+KEt_sc+0.25d0*Ip*KEr_p+KEr_sc)
743 ! write (iout,*) "KE_total",KE_total
745 end subroutine kinetic
746 !-----------------------------------------------------------------------------
748 !-----------------------------------------------------------------------------
750 !------------------------------------------------
751 ! The driver for molecular dynamics subroutines
752 !------------------------------------------------
755 use control, only:tcpu,ovrtim
756 ! use io_comm, only:ilen
758 use compare, only:secondary2,hairpin
759 use io, only:cartout,statout
760 ! implicit real*8 (a-h,o-z)
761 ! include 'DIMENSIONS'
764 integer :: IERROR,ERRCODE
766 ! include 'COMMON.SETUP'
767 ! include 'COMMON.CONTROL'
768 ! include 'COMMON.VAR'
769 ! include 'COMMON.MD'
771 ! include 'COMMON.LANGEVIN'
773 ! include 'COMMON.LANGEVIN.lang0'
775 ! include 'COMMON.CHAIN'
776 ! include 'COMMON.DERIV'
777 ! include 'COMMON.GEO'
778 ! include 'COMMON.LOCAL'
779 ! include 'COMMON.INTERACT'
780 ! include 'COMMON.IOUNITS'
781 ! include 'COMMON.NAMES'
782 ! include 'COMMON.TIME1'
783 ! include 'COMMON.HAIRPIN'
784 real(kind=8),dimension(3) :: L,vcm
786 real(kind=8),dimension(6*nres) :: v_work,v_transf !(maxres6) maxres6=6*maxres
788 integer :: rstcount !ilen,
790 character(len=50) :: tytul
791 !el common /gucio/ cm
792 integer :: itime,i,j,nharp
793 integer,dimension(4,nres/3) :: iharp !(4,nres/3)(4,maxres/3)
795 real(kind=8) :: tt0,scalfac
800 if (ilen(tmpdir).gt.0) &
801 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_" &
802 //liczba(:ilen(liczba))//'.rst')
804 if (ilen(tmpdir).gt.0) &
805 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"//'.rst')
812 write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
818 ! Determine the inverse of the inertia matrix.
819 call setup_MD_matrices
823 t_MDsetup = MPI_Wtime()-tt0
825 t_MDsetup = tcpu()-tt0
828 ! Entering the MD loop
834 if (lang.eq.2 .or. lang.eq.3) then
838 call sd_verlet_p_setup
840 call sd_verlet_ciccotti_setup
844 pfric0_mat(i,j,0)=pfric_mat(i,j)
845 afric0_mat(i,j,0)=afric_mat(i,j)
846 vfric0_mat(i,j,0)=vfric_mat(i,j)
847 prand0_mat(i,j,0)=prand_mat(i,j)
848 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
849 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
854 flag_stoch(i)=.false.
858 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
860 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
864 else if (lang.eq.1 .or. lang.eq.4) then
868 t_langsetup=MPI_Wtime()-tt0
871 t_langsetup=tcpu()-tt0
874 do itime=1,n_timestep
876 if (large.and. mod(itime,ntwe).eq.0) &
877 write (iout,*) "itime",itime
879 if (lang.gt.0 .and. surfarea .and. &
880 mod(itime,reset_fricmat).eq.0) then
881 if (lang.eq.2 .or. lang.eq.3) then
885 call sd_verlet_p_setup
887 call sd_verlet_ciccotti_setup
891 pfric0_mat(i,j,0)=pfric_mat(i,j)
892 afric0_mat(i,j,0)=afric_mat(i,j)
893 vfric0_mat(i,j,0)=vfric_mat(i,j)
894 prand0_mat(i,j,0)=prand_mat(i,j)
895 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
896 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
901 flag_stoch(i)=.false.
904 else if (lang.eq.1 .or. lang.eq.4) then
907 write (iout,'(a,i10)') &
908 "Friction matrix reset based on surface area, itime",itime
910 if (reset_vel .and. tbf .and. lang.eq.0 &
911 .and. mod(itime,count_reset_vel).eq.0) then
913 write(iout,'(a,f20.2)') &
914 "Velocities reset to random values, time",totT
917 d_t_old(j,i)=d_t(j,i)
921 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
925 d_t(j,0)=d_t(j,0)-vcm(j)
928 kinetic_T=2.0d0/(dimen3*Rb)*EK
929 scalfac=dsqrt(T_bath/kinetic_T)
930 write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
933 d_t_old(j,i)=scalfac*d_t(j,i)
939 ! Time-reversible RESPA algorithm
940 ! (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
941 call RESPA_step(itime)
943 ! Variable time step algorithm.
944 call velverlet_step(itime)
948 call brown_step(itime)
950 print *,"Brown dynamics not here!"
952 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
958 if (mod(itime,ntwe).eq.0) call statout(itime)
971 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
974 v_work(ind)=d_t(j,i+nres)
979 write (66,'(80f10.5)') &
980 ((d_t(j,i),j=1,3),i=0,nres-1),((d_t(j,i+nres),j=1,3),i=1,nres)
984 v_transf(i)=v_transf(i)+gvec(j,i)*v_work(j)
986 v_transf(i)= v_transf(i)*dsqrt(geigen(i))
988 write (67,'(80f10.5)') (v_transf(i),i=1,ind)
991 if (mod(itime,ntwx).eq.0) then
992 write (tytul,'("time",f8.2)') totT
994 call hairpin(.true.,nharp,iharp)
995 call secondary2(.true.)
996 call pdbout(potE,tytul,ipdb)
1001 if (rstcount.eq.1000.or.itime.eq.n_timestep) then
1002 open(irest2,file=rest2name,status='unknown')
1003 write(irest2,*) totT,EK,potE,totE,t_bath
1004 ! AL 4/17/17: Now writing d_t(0,:) too
1006 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
1008 ! AL 4/17/17: Now writing d_c(0,:) too
1010 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
1018 t_MD=MPI_Wtime()-tt0
1022 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))') &
1024 'MD calculations setup:',t_MDsetup,&
1025 'Energy & gradient evaluation:',t_enegrad,&
1026 'Stochastic MD setup:',t_langsetup,&
1027 'Stochastic MD step setup:',t_sdsetup,&
1029 write (iout,'(/28(1h=),a25,27(1h=))') &
1030 ' End of MD calculation '
1032 write (iout,*) "time for etotal",t_etotal," elong",t_elong,&
1034 write (iout,*) "time_fric",time_fric," time_stoch",time_stoch,&
1035 " time_fricmatmult",time_fricmatmult," time_fsample ",&
1040 !-----------------------------------------------------------------------------
1041 subroutine velverlet_step(itime)
1042 !-------------------------------------------------------------------------------
1043 ! Perform a single velocity Verlet step; the time step can be rescaled if
1044 ! increments in accelerations exceed the threshold
1045 !-------------------------------------------------------------------------------
1046 ! implicit real*8 (a-h,o-z)
1047 ! include 'DIMENSIONS'
1049 use control, only:tcpu
1053 integer :: ierror,ierrcode
1054 real(kind=8) :: errcode
1056 ! include 'COMMON.SETUP'
1057 ! include 'COMMON.VAR'
1058 ! include 'COMMON.MD'
1060 ! include 'COMMON.LANGEVIN'
1062 ! include 'COMMON.LANGEVIN.lang0'
1064 ! include 'COMMON.CHAIN'
1065 ! include 'COMMON.DERIV'
1066 ! include 'COMMON.GEO'
1067 ! include 'COMMON.LOCAL'
1068 ! include 'COMMON.INTERACT'
1069 ! include 'COMMON.IOUNITS'
1070 ! include 'COMMON.NAMES'
1071 ! include 'COMMON.TIME1'
1072 ! include 'COMMON.MUCA'
1073 real(kind=8),dimension(3) :: vcm,incr
1074 real(kind=8),dimension(3) :: L
1075 integer :: count,rstcount !ilen,
1077 character(len=50) :: tytul
1078 integer :: maxcount_scale = 20
1079 !el common /gucio/ cm
1080 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1081 !el common /stochcalc/ stochforcvec
1082 integer :: itime,icount_scale,itime_scal,i,j,ifac_time
1084 real(kind=8) :: epdrift,tt0,fac_time
1086 if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres)) !(MAXRES6) maxres6=6*maxres
1092 else if (lang.eq.2 .or. lang.eq.3) then
1094 call stochastic_force(stochforcvec)
1097 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
1099 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1106 icount_scale=icount_scale+1
1107 if (icount_scale.gt.maxcount_scale) then
1109 "ERROR: too many attempts at scaling down the time step. ",&
1110 "amax=",amax,"epdrift=",epdrift,&
1111 "damax=",damax,"edriftmax=",edriftmax,&
1115 call MPI_Abort(MPI_COMM_WORLD,IERROR,IERRCODE)
1119 ! First step of the velocity Verlet algorithm
1124 else if (lang.eq.3) then
1126 call sd_verlet1_ciccotti
1128 else if (lang.eq.1) then
1133 ! Build the chain from the newly calculated coordinates
1134 call chainbuild_cart
1135 if (rattle) call rattle1
1137 if (large.and. mod(itime,ntwe).eq.0) then
1138 write (iout,*) "Cartesian and internal coordinates: step 1"
1143 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),&
1144 (dc(j,i+nres),j=1,3)
1146 write (iout,*) "Accelerations"
1148 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1149 (d_a(j,i+nres),j=1,3)
1151 write (iout,*) "Velocities, step 1"
1153 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1154 (d_t(j,i+nres),j=1,3)
1163 ! Calculate energy and forces
1165 call etotal(potEcomp)
1166 if (large.and. mod(itime,ntwe).eq.0) &
1167 call enerprint(potEcomp)
1170 t_etotal=t_etotal+MPI_Wtime()-tt0
1172 t_etotal=t_etotal+tcpu()-tt0
1175 potE=potEcomp(0)-potEcomp(20)
1177 ! Get the new accelerations
1180 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1182 t_enegrad=t_enegrad+tcpu()-tt0
1184 ! Determine maximum acceleration and scale down the timestep if needed
1186 amax=amax/(itime_scal+1)**2
1187 call predict_edrift(epdrift)
1188 if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then
1189 ! Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
1191 ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax)) &
1193 itime_scal=itime_scal+ifac_time
1194 ! fac_time=dmin1(damax/amax,0.5d0)
1195 fac_time=0.5d0**ifac_time
1196 d_time=d_time*fac_time
1197 if (lang.eq.2 .or. lang.eq.3) then
1199 ! write (iout,*) "Calling sd_verlet_setup: 1"
1200 ! Rescale the stochastic forces and recalculate or restore
1201 ! the matrices of tinker integrator
1202 if (itime_scal.gt.maxflag_stoch) then
1203 if (large) write (iout,'(a,i5,a)') &
1204 "Calculate matrices for stochastic step;",&
1205 " itime_scal ",itime_scal
1207 call sd_verlet_p_setup
1209 call sd_verlet_ciccotti_setup
1211 write (iout,'(2a,i3,a,i3,1h.)') &
1212 "Warning: cannot store matrices for stochastic",&
1213 " integration because the index",itime_scal,&
1214 " is greater than",maxflag_stoch
1215 write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",&
1216 " integration Langevin algorithm for better efficiency."
1217 else if (flag_stoch(itime_scal)) then
1218 if (large) write (iout,'(a,i5,a,l1)') &
1219 "Restore matrices for stochastic step; itime_scal ",&
1220 itime_scal," flag ",flag_stoch(itime_scal)
1223 pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
1224 afric_mat(i,j)=afric0_mat(i,j,itime_scal)
1225 vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
1226 prand_mat(i,j)=prand0_mat(i,j,itime_scal)
1227 vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
1228 vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
1232 if (large) write (iout,'(2a,i5,a,l1)') &
1233 "Calculate & store matrices for stochastic step;",&
1234 " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
1236 call sd_verlet_p_setup
1238 call sd_verlet_ciccotti_setup
1240 flag_stoch(ifac_time)=.true.
1243 pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
1244 afric0_mat(i,j,itime_scal)=afric_mat(i,j)
1245 vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
1246 prand0_mat(i,j,itime_scal)=prand_mat(i,j)
1247 vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
1248 vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
1252 fac_time=1.0d0/dsqrt(fac_time)
1254 stochforcvec(i)=fac_time*stochforcvec(i)
1257 else if (lang.eq.1) then
1258 ! Rescale the accelerations due to stochastic forces
1259 fac_time=1.0d0/dsqrt(fac_time)
1261 d_as_work(i)=d_as_work(i)*fac_time
1264 if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)') &
1265 "itime",itime," Timestep scaled down to ",&
1266 d_time," ifac_time",ifac_time," itime_scal",itime_scal
1268 ! Second step of the velocity Verlet algorithm
1273 else if (lang.eq.3) then
1275 call sd_verlet2_ciccotti
1277 else if (lang.eq.1) then
1282 if (rattle) call rattle2
1284 if (d_time.ne.d_time0) then
1287 if (lang.eq.2 .or. lang.eq.3) then
1288 if (large) write (iout,'(a)') &
1289 "Restore original matrices for stochastic step"
1290 ! write (iout,*) "Calling sd_verlet_setup: 2"
1291 ! Restore the matrices of tinker integrator if the time step has been restored
1294 pfric_mat(i,j)=pfric0_mat(i,j,0)
1295 afric_mat(i,j)=afric0_mat(i,j,0)
1296 vfric_mat(i,j)=vfric0_mat(i,j,0)
1297 prand_mat(i,j)=prand0_mat(i,j,0)
1298 vrand_mat1(i,j)=vrand0_mat1(i,j,0)
1299 vrand_mat2(i,j)=vrand0_mat2(i,j,0)
1308 ! Calculate the kinetic and the total energy and the kinetic temperature
1312 ! call kinetic1(EK1)
1313 ! write (iout,*) "step",itime," EK",EK," EK1",EK1
1315 ! Couple the system to Berendsen bath if needed
1316 if (tbf .and. lang.eq.0) then
1319 kinetic_T=2.0d0/(dimen3*Rb)*EK
1320 ! Backup the coordinates, velocities, and accelerations
1324 d_t_old(j,i)=d_t(j,i)
1325 d_a_old(j,i)=d_a(j,i)
1329 if (mod(itime,ntwe).eq.0 .and. large) then
1330 write (iout,*) "Velocities, step 2"
1332 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1333 (d_t(j,i+nres),j=1,3)
1338 end subroutine velverlet_step
1339 !-----------------------------------------------------------------------------
1340 subroutine RESPA_step(itime)
1341 !-------------------------------------------------------------------------------
1342 ! Perform a single RESPA step.
1343 !-------------------------------------------------------------------------------
1344 ! implicit real*8 (a-h,o-z)
1345 ! include 'DIMENSIONS'
1349 use control, only:tcpu
1351 ! use io_conf, only:cartprint
1354 integer :: IERROR,ERRCODE
1356 ! include 'COMMON.SETUP'
1357 ! include 'COMMON.CONTROL'
1358 ! include 'COMMON.VAR'
1359 ! include 'COMMON.MD'
1361 ! include 'COMMON.LANGEVIN'
1363 ! include 'COMMON.LANGEVIN.lang0'
1365 ! include 'COMMON.CHAIN'
1366 ! include 'COMMON.DERIV'
1367 ! include 'COMMON.GEO'
1368 ! include 'COMMON.LOCAL'
1369 ! include 'COMMON.INTERACT'
1370 ! include 'COMMON.IOUNITS'
1371 ! include 'COMMON.NAMES'
1372 ! include 'COMMON.TIME1'
1373 real(kind=8),dimension(0:n_ene) :: energia_short,energia_long
1374 real(kind=8),dimension(3) :: L,vcm,incr
1375 real(kind=8),dimension(3,0:2*nres) :: dc_old0,d_t_old0,d_a_old0 !(3,0:maxres2) maxres2=2*maxres
1376 logical :: PRINT_AMTS_MSG = .false.
1377 integer :: count,rstcount !ilen,
1379 character(len=50) :: tytul
1380 integer :: maxcount_scale = 10
1381 !el common /gucio/ cm
1382 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1383 !el common /stochcalc/ stochforcvec
1384 integer :: itime,itt,i,j,itsplit
1386 !el common /cipiszcze/ itt
1388 real(kind=8) :: epdrift,tt0,epdriftmax
1391 if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres)) !(MAXRES6) maxres6=6*maxres
1395 if (large.and. mod(itime,ntwe).eq.0) then
1396 write (iout,*) "***************** RESPA itime",itime
1397 write (iout,*) "Cartesian and internal coordinates: step 0"
1399 call pdbout(0.0d0,"cipiszcze",iout)
1401 write (iout,*) "Accelerations from long-range forces"
1403 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1404 (d_a(j,i+nres),j=1,3)
1406 write (iout,*) "Velocities, step 0"
1408 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1409 (d_t(j,i+nres),j=1,3)
1414 ! Perform the initial RESPA step (increment velocities)
1415 ! write (iout,*) "*********************** RESPA ini"
1418 if (mod(itime,ntwe).eq.0 .and. large) then
1419 write (iout,*) "Velocities, end"
1421 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1422 (d_t(j,i+nres),j=1,3)
1426 ! Compute the short-range forces
1432 ! 7/2/2009 commented out
1434 ! call etotal_short(energia_short)
1437 ! 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
1440 d_a(j,i)=d_a_short(j,i)
1444 if (large.and. mod(itime,ntwe).eq.0) then
1445 write (iout,*) "energia_short",energia_short(0)
1446 write (iout,*) "Accelerations from short-range forces"
1448 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1449 (d_a(j,i+nres),j=1,3)
1454 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1456 t_enegrad=t_enegrad+tcpu()-tt0
1461 d_t_old(j,i)=d_t(j,i)
1462 d_a_old(j,i)=d_a(j,i)
1465 ! 6/30/08 A-MTS: attempt at increasing the split number
1468 dc_old0(j,i)=dc_old(j,i)
1469 d_t_old0(j,i)=d_t_old(j,i)
1470 d_a_old0(j,i)=d_a_old(j,i)
1473 if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
1474 if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
1481 ! write (iout,*) "itime",itime," ntime_split",ntime_split
1482 ! Split the time step
1483 d_time=d_time0/ntime_split
1484 ! Perform the short-range RESPA steps (velocity Verlet increments of
1485 ! positions and velocities using short-range forces)
1486 ! write (iout,*) "*********************** RESPA split"
1487 do itsplit=1,ntime_split
1490 else if (lang.eq.2 .or. lang.eq.3) then
1492 call stochastic_force(stochforcvec)
1495 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
1497 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1502 ! First step of the velocity Verlet algorithm
1507 else if (lang.eq.3) then
1509 call sd_verlet1_ciccotti
1511 else if (lang.eq.1) then
1516 ! Build the chain from the newly calculated coordinates
1517 call chainbuild_cart
1518 if (rattle) call rattle1
1520 if (large.and. mod(itime,ntwe).eq.0) then
1521 write (iout,*) "***** ITSPLIT",itsplit
1522 write (iout,*) "Cartesian and internal coordinates: step 1"
1523 call pdbout(0.0d0,"cipiszcze",iout)
1526 write (iout,*) "Velocities, step 1"
1528 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1529 (d_t(j,i+nres),j=1,3)
1538 ! Calculate energy and forces
1540 call etotal_short(energia_short)
1541 if (large.and. mod(itime,ntwe).eq.0) &
1542 call enerprint(energia_short)
1545 t_eshort=t_eshort+MPI_Wtime()-tt0
1547 t_eshort=t_eshort+tcpu()-tt0
1551 ! Get the new accelerations
1553 ! 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
1556 d_a_short(j,i)=d_a(j,i)
1560 if (large.and. mod(itime,ntwe).eq.0) then
1561 write (iout,*)"energia_short",energia_short(0)
1562 write (iout,*) "Accelerations from short-range forces"
1564 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1565 (d_a(j,i+nres),j=1,3)
1570 ! Determine maximum acceleration and scale down the timestep if needed
1572 amax=amax/ntime_split**2
1573 call predict_edrift(epdrift)
1574 if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0) &
1575 write (iout,*) "amax",amax," damax",damax,&
1576 " epdrift",epdrift," epdriftmax",epdriftmax
1577 ! Exit loop and try with increased split number if the change of
1578 ! acceleration is too big
1579 if (amax.gt.damax .or. epdrift.gt.edriftmax) then
1580 if (ntime_split.lt.maxtime_split) then
1582 ntime_split=ntime_split*2
1585 dc_old(j,i)=dc_old0(j,i)
1586 d_t_old(j,i)=d_t_old0(j,i)
1587 d_a_old(j,i)=d_a_old0(j,i)
1590 if (PRINT_AMTS_MSG) then
1591 write (iout,*) "acceleration/energy drift too large",amax,&
1592 epdrift," split increased to ",ntime_split," itime",itime,&
1598 "Uh-hu. Bumpy landscape. Maximum splitting number",&
1600 " already reached!!! Trying to carry on!"
1604 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1606 t_enegrad=t_enegrad+tcpu()-tt0
1608 ! Second step of the velocity Verlet algorithm
1613 else if (lang.eq.3) then
1615 call sd_verlet2_ciccotti
1617 else if (lang.eq.1) then
1622 if (rattle) call rattle2
1623 ! Backup the coordinates, velocities, and accelerations
1627 d_t_old(j,i)=d_t(j,i)
1628 d_a_old(j,i)=d_a(j,i)
1635 ! Restore the time step
1637 ! Compute long-range forces
1644 call etotal_long(energia_long)
1645 if (large.and. mod(itime,ntwe).eq.0) &
1646 call enerprint(energia_long)
1649 t_elong=t_elong+MPI_Wtime()-tt0
1651 t_elong=t_elong+tcpu()-tt0
1657 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1659 t_enegrad=t_enegrad+tcpu()-tt0
1661 ! Compute accelerations from long-range forces
1663 if (large.and. mod(itime,ntwe).eq.0) then
1664 write (iout,*) "energia_long",energia_long(0)
1665 write (iout,*) "Cartesian and internal coordinates: step 2"
1667 call pdbout(0.0d0,"cipiszcze",iout)
1669 write (iout,*) "Accelerations from long-range forces"
1671 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1672 (d_a(j,i+nres),j=1,3)
1674 write (iout,*) "Velocities, step 2"
1676 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1677 (d_t(j,i+nres),j=1,3)
1681 ! Compute the final RESPA step (increment velocities)
1682 ! write (iout,*) "*********************** RESPA fin"
1684 ! Compute the complete potential energy
1686 potEcomp(i)=energia_short(i)+energia_long(i)
1688 potE=potEcomp(0)-potEcomp(20)
1689 ! potE=energia_short(0)+energia_long(0)
1691 ! Calculate the kinetic and the total energy and the kinetic temperature
1694 ! Couple the system to Berendsen bath if needed
1695 if (tbf .and. lang.eq.0) then
1698 kinetic_T=2.0d0/(dimen3*Rb)*EK
1699 ! Backup the coordinates, velocities, and accelerations
1701 if (mod(itime,ntwe).eq.0 .and. large) then
1702 write (iout,*) "Velocities, end"
1704 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1705 (d_t(j,i+nres),j=1,3)
1710 end subroutine RESPA_step
1711 !-----------------------------------------------------------------------------
1712 subroutine RESPA_vel
1713 ! First and last RESPA step (incrementing velocities using long-range
1716 ! implicit real*8 (a-h,o-z)
1717 ! include 'DIMENSIONS'
1718 ! include 'COMMON.CONTROL'
1719 ! include 'COMMON.VAR'
1720 ! include 'COMMON.MD'
1721 ! include 'COMMON.CHAIN'
1722 ! include 'COMMON.DERIV'
1723 ! include 'COMMON.GEO'
1724 ! include 'COMMON.LOCAL'
1725 ! include 'COMMON.INTERACT'
1726 ! include 'COMMON.IOUNITS'
1727 ! include 'COMMON.NAMES'
1728 integer :: i,j,inres
1731 d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
1735 d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
1739 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1742 d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
1747 end subroutine RESPA_vel
1748 !-----------------------------------------------------------------------------
1750 ! Applying velocity Verlet algorithm - step 1 to coordinates
1752 ! implicit real*8 (a-h,o-z)
1753 ! include 'DIMENSIONS'
1754 ! include 'COMMON.CONTROL'
1755 ! include 'COMMON.VAR'
1756 ! include 'COMMON.MD'
1757 ! include 'COMMON.CHAIN'
1758 ! include 'COMMON.DERIV'
1759 ! include 'COMMON.GEO'
1760 ! include 'COMMON.LOCAL'
1761 ! include 'COMMON.INTERACT'
1762 ! include 'COMMON.IOUNITS'
1763 ! include 'COMMON.NAMES'
1764 real(kind=8) :: adt,adt2
1765 integer :: i,j,inres
1768 write (iout,*) "VELVERLET1 START: DC"
1770 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),&
1771 (dc(j,i+nres),j=1,3)
1775 adt=d_a_old(j,0)*d_time
1777 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1778 d_t_new(j,0)=d_t_old(j,0)+adt2
1779 d_t(j,0)=d_t_old(j,0)+adt
1783 adt=d_a_old(j,i)*d_time
1785 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1786 d_t_new(j,i)=d_t_old(j,i)+adt2
1787 d_t(j,i)=d_t_old(j,i)+adt
1791 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1794 adt=d_a_old(j,inres)*d_time
1796 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1797 d_t_new(j,inres)=d_t_old(j,inres)+adt2
1798 d_t(j,inres)=d_t_old(j,inres)+adt
1803 write (iout,*) "VELVERLET1 END: DC"
1805 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),&
1806 (dc(j,i+nres),j=1,3)
1810 end subroutine verlet1
1811 !-----------------------------------------------------------------------------
1813 ! Step 2 of the velocity Verlet algorithm: update velocities
1815 ! implicit real*8 (a-h,o-z)
1816 ! include 'DIMENSIONS'
1817 ! include 'COMMON.CONTROL'
1818 ! include 'COMMON.VAR'
1819 ! include 'COMMON.MD'
1820 ! include 'COMMON.CHAIN'
1821 ! include 'COMMON.DERIV'
1822 ! include 'COMMON.GEO'
1823 ! include 'COMMON.LOCAL'
1824 ! include 'COMMON.INTERACT'
1825 ! include 'COMMON.IOUNITS'
1826 ! include 'COMMON.NAMES'
1827 integer :: i,j,inres
1830 d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
1834 d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
1838 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1841 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
1846 end subroutine verlet2
1847 !-----------------------------------------------------------------------------
1848 subroutine sddir_precalc
1849 ! Applying velocity Verlet algorithm - step 1 to coordinates
1850 ! implicit real*8 (a-h,o-z)
1851 ! include 'DIMENSIONS'
1857 ! include 'COMMON.CONTROL'
1858 ! include 'COMMON.VAR'
1859 ! include 'COMMON.MD'
1861 ! include 'COMMON.LANGEVIN'
1863 ! include 'COMMON.LANGEVIN.lang0'
1865 ! include 'COMMON.CHAIN'
1866 ! include 'COMMON.DERIV'
1867 ! include 'COMMON.GEO'
1868 ! include 'COMMON.LOCAL'
1869 ! include 'COMMON.INTERACT'
1870 ! include 'COMMON.IOUNITS'
1871 ! include 'COMMON.NAMES'
1872 ! include 'COMMON.TIME1'
1873 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1874 !el common /stochcalc/ stochforcvec
1875 real(kind=8) :: time00
1877 ! Compute friction and stochastic forces
1882 time_fric=time_fric+MPI_Wtime()-time00
1884 call stochastic_force(stochforcvec)
1885 time_stoch=time_stoch+MPI_Wtime()-time00
1888 ! Compute the acceleration due to friction forces (d_af_work) and stochastic
1889 ! forces (d_as_work)
1891 call ginv_mult(fric_work, d_af_work)
1892 call ginv_mult(stochforcvec, d_as_work)
1894 end subroutine sddir_precalc
1895 !-----------------------------------------------------------------------------
1896 subroutine sddir_verlet1
1897 ! Applying velocity Verlet algorithm - step 1 to velocities
1900 ! implicit real*8 (a-h,o-z)
1901 ! include 'DIMENSIONS'
1902 ! include 'COMMON.CONTROL'
1903 ! include 'COMMON.VAR'
1904 ! include 'COMMON.MD'
1906 ! include 'COMMON.LANGEVIN'
1908 ! include 'COMMON.LANGEVIN.lang0'
1910 ! include 'COMMON.CHAIN'
1911 ! include 'COMMON.DERIV'
1912 ! include 'COMMON.GEO'
1913 ! include 'COMMON.LOCAL'
1914 ! include 'COMMON.INTERACT'
1915 ! include 'COMMON.IOUNITS'
1916 ! include 'COMMON.NAMES'
1917 ! Revised 3/31/05 AL: correlation between random contributions to
1918 ! position and velocity increments included.
1919 real(kind=8) :: sqrt13 = 0.57735026918962576451d0 ! 1/sqrt(3)
1920 real(kind=8) :: adt,adt2
1921 integer :: i,j,ind,inres
1923 ! Add the contribution from BOTH friction and stochastic force to the
1924 ! coordinates, but ONLY the contribution from the friction forces to velocities
1927 adt=(d_a_old(j,0)+d_af_work(j))*d_time
1928 adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
1929 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1930 d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
1931 d_t(j,0)=d_t_old(j,0)+adt
1936 adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
1937 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1938 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1939 d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
1940 d_t(j,i)=d_t_old(j,i)+adt
1945 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1948 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
1949 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1950 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1951 d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
1952 d_t(j,inres)=d_t_old(j,inres)+adt
1958 end subroutine sddir_verlet1
1959 !-----------------------------------------------------------------------------
1960 subroutine sddir_verlet2
1961 ! Calculating the adjusted velocities for accelerations
1964 ! implicit real*8 (a-h,o-z)
1965 ! include 'DIMENSIONS'
1966 ! include 'COMMON.CONTROL'
1967 ! include 'COMMON.VAR'
1968 ! include 'COMMON.MD'
1970 ! include 'COMMON.LANGEVIN'
1972 ! include 'COMMON.LANGEVIN.lang0'
1974 ! include 'COMMON.CHAIN'
1975 ! include 'COMMON.DERIV'
1976 ! include 'COMMON.GEO'
1977 ! include 'COMMON.LOCAL'
1978 ! include 'COMMON.INTERACT'
1979 ! include 'COMMON.IOUNITS'
1980 ! include 'COMMON.NAMES'
1981 real(kind=8),dimension(6*nres) :: stochforcvec,d_as_work1 !(MAXRES6) maxres6=6*maxres
1982 real(kind=8) :: cos60 = 0.5d0, sin60 = 0.86602540378443864676d0
1983 integer :: i,j,ind,inres
1984 ! Revised 3/31/05 AL: correlation between random contributions to
1985 ! position and velocity increments included.
1986 ! The correlation coefficients are calculated at low-friction limit.
1987 ! Also, friction forces are now not calculated with new velocities.
1989 ! call friction_force
1990 call stochastic_force(stochforcvec)
1992 ! Compute the acceleration due to friction forces (d_af_work) and stochastic
1993 ! forces (d_as_work)
1995 call ginv_mult(stochforcvec, d_as_work1)
2001 d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j)) &
2002 +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
2007 d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j)) &
2008 +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
2013 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2016 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres) &
2017 +d_af_work(ind+j))+sin60*d_as_work(ind+j) &
2018 +cos60*d_as_work1(ind+j))*d_time
2024 end subroutine sddir_verlet2
2025 !-----------------------------------------------------------------------------
2026 subroutine max_accel
2028 ! Find the maximum difference in the accelerations of the the sites
2029 ! at the beginning and the end of the time step.
2032 ! implicit real*8 (a-h,o-z)
2033 ! include 'DIMENSIONS'
2034 ! include 'COMMON.CONTROL'
2035 ! include 'COMMON.VAR'
2036 ! include 'COMMON.MD'
2037 ! include 'COMMON.CHAIN'
2038 ! include 'COMMON.DERIV'
2039 ! include 'COMMON.GEO'
2040 ! include 'COMMON.LOCAL'
2041 ! include 'COMMON.INTERACT'
2042 ! include 'COMMON.IOUNITS'
2043 real(kind=8),dimension(3) :: aux,accel,accel_old
2044 real(kind=8) :: dacc
2048 ! aux(j)=d_a(j,0)-d_a_old(j,0)
2049 accel_old(j)=d_a_old(j,0)
2056 ! 7/3/08 changed to asymmetric difference
2058 ! accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
2059 accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
2060 accel(j)=accel(j)+0.5d0*d_a(j,i)
2061 ! if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
2062 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
2063 dacc=dabs(accel(j)-accel_old(j))
2064 ! write (iout,*) i,dacc
2065 if (dacc.gt.amax) amax=dacc
2073 accel_old(j)=d_a_old(j,0)
2078 accel_old(j)=accel_old(j)+d_a_old(j,1)
2079 accel(j)=accel(j)+d_a(j,1)
2083 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2085 ! accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
2086 accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
2087 accel(j)=accel(j)+d_a(j,i+nres)
2091 ! if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
2092 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
2093 dacc=dabs(accel(j)-accel_old(j))
2094 ! write (iout,*) "side-chain",i,dacc
2095 if (dacc.gt.amax) amax=dacc
2099 accel_old(j)=accel_old(j)+d_a_old(j,i)
2100 accel(j)=accel(j)+d_a(j,i)
2101 ! aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
2105 end subroutine max_accel
2106 !-----------------------------------------------------------------------------
2107 subroutine predict_edrift(epdrift)
2109 ! Predict the drift of the potential energy
2112 use control_data, only: lmuca
2113 ! implicit real*8 (a-h,o-z)
2114 ! include 'DIMENSIONS'
2115 ! include 'COMMON.CONTROL'
2116 ! include 'COMMON.VAR'
2117 ! include 'COMMON.MD'
2118 ! include 'COMMON.CHAIN'
2119 ! include 'COMMON.DERIV'
2120 ! include 'COMMON.GEO'
2121 ! include 'COMMON.LOCAL'
2122 ! include 'COMMON.INTERACT'
2123 ! include 'COMMON.IOUNITS'
2124 ! include 'COMMON.MUCA'
2125 real(kind=8) :: epdrift,epdriftij
2127 ! Drift of the potential energy
2133 epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
2134 if (lmuca) epdriftij=epdriftij*factor
2135 ! write (iout,*) "back",i,j,epdriftij
2136 if (epdriftij.gt.epdrift) epdrift=epdriftij
2140 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2143 dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
2144 if (lmuca) epdriftij=epdriftij*factor
2145 ! write (iout,*) "side",i,j,epdriftij
2146 if (epdriftij.gt.epdrift) epdrift=epdriftij
2150 epdrift=0.5d0*epdrift*d_time*d_time
2151 ! write (iout,*) "epdrift",epdrift
2153 end subroutine predict_edrift
2154 !-----------------------------------------------------------------------------
2155 subroutine verlet_bath
2157 ! Coupling to the thermostat by using the Berendsen algorithm
2160 ! implicit real*8 (a-h,o-z)
2161 ! include 'DIMENSIONS'
2162 ! include 'COMMON.CONTROL'
2163 ! include 'COMMON.VAR'
2164 ! include 'COMMON.MD'
2165 ! include 'COMMON.CHAIN'
2166 ! include 'COMMON.DERIV'
2167 ! include 'COMMON.GEO'
2168 ! include 'COMMON.LOCAL'
2169 ! include 'COMMON.INTERACT'
2170 ! include 'COMMON.IOUNITS'
2171 ! include 'COMMON.NAMES'
2172 real(kind=8) :: T_half,fact
2173 integer :: i,j,inres
2175 T_half=2.0d0/(dimen3*Rb)*EK
2176 fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
2177 ! write(iout,*) "T_half", T_half
2178 ! write(iout,*) "EK", EK
2179 ! write(iout,*) "fact", fact
2181 d_t(j,0)=fact*d_t(j,0)
2185 d_t(j,i)=fact*d_t(j,i)
2189 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2192 d_t(j,inres)=fact*d_t(j,inres)
2197 end subroutine verlet_bath
2198 !-----------------------------------------------------------------------------
2200 ! Set up the initial conditions of a MD simulation
2203 use control, only:tcpu
2204 !el use io_basic, only:ilen
2207 use minimm, only:minim_dc,minimize,sc_move
2208 use io_config, only:readrst
2209 use io, only:statout
2210 ! implicit real*8 (a-h,o-z)
2211 ! include 'DIMENSIONS'
2214 character(len=16) :: form
2215 integer :: IERROR,ERRCODE
2217 ! include 'COMMON.SETUP'
2218 ! include 'COMMON.CONTROL'
2219 ! include 'COMMON.VAR'
2220 ! include 'COMMON.MD'
2222 ! include 'COMMON.LANGEVIN'
2224 ! include 'COMMON.LANGEVIN.lang0'
2226 ! include 'COMMON.CHAIN'
2227 ! include 'COMMON.DERIV'
2228 ! include 'COMMON.GEO'
2229 ! include 'COMMON.LOCAL'
2230 ! include 'COMMON.INTERACT'
2231 ! include 'COMMON.IOUNITS'
2232 ! include 'COMMON.NAMES'
2233 ! include 'COMMON.REMD'
2234 real(kind=8),dimension(0:n_ene) :: energia_long,energia_short
2235 real(kind=8),dimension(3) :: vcm,incr,L
2236 real(kind=8) :: xv,sigv,lowb,highb
2237 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
2238 character(len=256) :: qstr
2241 character(len=50) :: tytul
2242 logical :: file_exist
2243 !el common /gucio/ cm
2244 integer :: i,j,ipos,iq,iw,nft_sc,iretcode,nfun,itime,ierr
2245 real(kind=8) :: etot,tt0
2249 ! write(iout,*) "d_time", d_time
2250 ! Compute the standard deviations of stochastic forces for Langevin dynamics
2251 ! if the friction coefficients do not depend on surface area
2252 if (lang.gt.0 .and. .not.surfarea) then
2254 stdforcp(i)=stdfp*dsqrt(gamp)
2257 stdforcsc(i)=stdfsc(iabs(itype(i))) &
2258 *dsqrt(gamsc(iabs(itype(i))))
2261 ! Open the pdb file for snapshotshots
2264 if (ilen(tmpdir).gt.0) &
2265 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2266 liczba(:ilen(liczba))//".pdb")
2268 file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2272 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) &
2273 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2274 liczba(:ilen(liczba))//".x")
2275 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2278 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) &
2279 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2280 liczba(:ilen(liczba))//".cx")
2281 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2287 if (ilen(tmpdir).gt.0) &
2288 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
2289 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
2291 if (ilen(tmpdir).gt.0) &
2292 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
2293 cartname=prefix(:ilen(prefix))//"_MD.cx"
2297 write (qstr,'(256(1h ))')
2300 iq = qinfrag(i,iset)*10
2301 iw = wfrag(i,iset)/100
2303 if(me.eq.king.or..not.out1file) &
2304 write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
2305 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
2310 iq = qinpair(i,iset)*10
2311 iw = wpair(i,iset)/100
2313 if(me.eq.king.or..not.out1file) &
2314 write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
2315 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
2319 ! pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
2321 ! cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
2323 ! cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
2325 ! statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
2329 if (restart1file) then
2331 inquire(file=mremd_rst_name,exist=file_exist)
2332 write (*,*) me," Before broadcast: file_exist",file_exist
2334 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,&
2337 write (*,*) me," After broadcast: file_exist",file_exist
2338 ! inquire(file=mremd_rst_name,exist=file_exist)
2339 if(me.eq.king.or..not.out1file) &
2340 write(iout,*) "Initial state read by master and distributed"
2342 if (ilen(tmpdir).gt.0) &
2343 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_' &
2344 //liczba(:ilen(liczba))//'.rst')
2345 inquire(file=rest2name,exist=file_exist)
2348 if(.not.restart1file) then
2349 if(me.eq.king.or..not.out1file) &
2350 write(iout,*) "Initial state will be read from file ",&
2351 rest2name(:ilen(rest2name))
2354 call rescale_weights(t_bath)
2356 if(me.eq.king.or..not.out1file)then
2357 if (restart1file) then
2358 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),&
2361 write(iout,*) "File ",rest2name(:ilen(rest2name)),&
2364 write(iout,*) "Initial velocities randomly generated"
2370 ! Generate initial velocities
2371 if(me.eq.king.or..not.out1file) &
2372 write(iout,*) "Initial velocities randomly generated"
2376 ! rest2name = prefix(:ilen(prefix))//'.rst'
2377 if(me.eq.king.or..not.out1file)then
2378 write (iout,*) "Initial velocities"
2380 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2381 (d_t(j,i+nres),j=1,3)
2383 ! Zeroing the total angular momentum of the system
2384 write(iout,*) "Calling the zero-angular momentum subroutine"
2387 ! Getting the potential energy and forces and velocities and accelerations
2389 ! write (iout,*) "velocity of the center of the mass:"
2390 ! write (iout,*) (vcm(j),j=1,3)
2392 d_t(j,0)=d_t(j,0)-vcm(j)
2394 ! Removing the velocity of the center of mass
2396 if(me.eq.king.or..not.out1file)then
2397 write (iout,*) "vcm right after adjustment:"
2398 write (iout,*) (vcm(j),j=1,3)
2402 if(iranconf.ne.0) then
2404 print *, 'Calling OVERLAP_SC'
2405 call overlap_sc(fail)
2408 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
2409 print *,'SC_move',nft_sc,etot
2410 if(me.eq.king.or..not.out1file) &
2411 write(iout,*) 'SC_move',nft_sc,etot
2415 print *, 'Calling MINIM_DC'
2416 call minim_dc(etot,iretcode,nfun)
2418 call geom_to_var(nvar,varia)
2419 print *,'Calling MINIMIZE.'
2420 call minimize(etot,varia,iretcode,nfun)
2421 call var_to_geom(nvar,varia)
2423 if(me.eq.king.or..not.out1file) &
2424 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2427 call chainbuild_cart
2432 kinetic_T=2.0d0/(dimen3*Rb)*EK
2433 if(me.eq.king.or..not.out1file)then
2443 call etotal(potEcomp)
2444 if (large) call enerprint(potEcomp)
2447 t_etotal=t_etotal+MPI_Wtime()-tt0
2449 t_etotal=t_etotal+tcpu()-tt0
2456 if (amax*d_time .gt. dvmax) then
2457 d_time=d_time*dvmax/amax
2458 if(me.eq.king.or..not.out1file) write (iout,*) &
2459 "Time step reduced to",d_time,&
2460 " because of too large initial acceleration."
2462 if(me.eq.king.or..not.out1file)then
2463 write(iout,*) "Potential energy and its components"
2464 call enerprint(potEcomp)
2465 ! write(iout,*) (potEcomp(i),i=0,n_ene)
2467 potE=potEcomp(0)-potEcomp(20)
2470 if (ntwe.ne.0) call statout(itime)
2471 if(me.eq.king.or..not.out1file) &
2472 write (iout,'(/a/3(a25,1pe14.5/))') "Initial:", &
2473 " Kinetic energy",EK," Potential energy",potE, &
2474 " Total energy",totE," Maximum acceleration ", &
2477 write (iout,*) "Initial coordinates"
2479 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),&
2482 write (iout,*) "Initial dC"
2484 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),&
2485 (dc(j,i+nres),j=1,3)
2487 write (iout,*) "Initial velocities"
2488 write (iout,"(13x,' backbone ',23x,' side chain')")
2490 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2491 (d_t(j,i+nres),j=1,3)
2493 write (iout,*) "Initial accelerations"
2495 ! write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2496 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),&
2497 (d_a(j,i+nres),j=1,3)
2503 d_t_old(j,i)=d_t(j,i)
2504 d_a_old(j,i)=d_a(j,i)
2506 ! write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
2515 call etotal_short(energia_short)
2516 if (large) call enerprint(potEcomp)
2519 t_eshort=t_eshort+MPI_Wtime()-tt0
2521 t_eshort=t_eshort+tcpu()-tt0
2526 if(.not.out1file .and. large) then
2527 write (iout,*) "energia_long",energia_long(0),&
2528 " energia_short",energia_short(0),&
2529 " total",energia_long(0)+energia_short(0)
2530 write (iout,*) "Initial fast-force accelerations"
2532 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
2533 (d_a(j,i+nres),j=1,3)
2536 ! 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
2539 d_a_short(j,i)=d_a(j,i)
2548 call etotal_long(energia_long)
2549 if (large) call enerprint(potEcomp)
2552 t_elong=t_elong+MPI_Wtime()-tt0
2554 t_elong=t_elong+tcpu()-tt0
2559 if(.not.out1file .and. large) then
2560 write (iout,*) "energia_long",energia_long(0)
2561 write (iout,*) "Initial slow-force accelerations"
2563 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
2564 (d_a(j,i+nres),j=1,3)
2568 t_enegrad=t_enegrad+MPI_Wtime()-tt0
2570 t_enegrad=t_enegrad+tcpu()-tt0
2574 end subroutine init_MD
2575 !-----------------------------------------------------------------------------
2576 subroutine random_vel
2578 ! implicit real*8 (a-h,o-z)
2580 use random, only:anorm_distr
2582 ! include 'DIMENSIONS'
2583 ! include 'COMMON.CONTROL'
2584 ! include 'COMMON.VAR'
2585 ! include 'COMMON.MD'
2587 ! include 'COMMON.LANGEVIN'
2589 ! include 'COMMON.LANGEVIN.lang0'
2591 ! include 'COMMON.CHAIN'
2592 ! include 'COMMON.DERIV'
2593 ! include 'COMMON.GEO'
2594 ! include 'COMMON.LOCAL'
2595 ! include 'COMMON.INTERACT'
2596 ! include 'COMMON.IOUNITS'
2597 ! include 'COMMON.NAMES'
2598 ! include 'COMMON.TIME1'
2599 real(kind=8) :: xv,sigv,lowb,highb ,Ek1
2601 real(kind=8) ,allocatable, dimension(:) :: DDU1,DDU2,DL2,DL1,xsolv,DML,rs
2602 real(kind=8) :: sumx
2604 real(kind=8) ,allocatable, dimension(:) :: rsold
2605 real (kind=8),allocatable,dimension(:,:) :: matold
2608 integer :: i,j,ii,k,ind,mark,imark
2609 ! Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
2610 ! First generate velocities in the eigenspace of the G matrix
2611 ! write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
2618 sigv=dsqrt((Rb*t_bath)/geigen(i))
2621 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
2623 write (iout,*) "i",i," ii",ii," geigen",geigen(i),&
2624 " d_t_work_new",d_t_work_new(ii)
2635 Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
2638 write (iout,*) "Ek from eigenvectors",Ek1
2639 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
2643 ! Transform velocities to UNRES coordinate space
2644 allocate (DL1(2*nres))
2645 allocate (DDU1(2*nres))
2646 allocate (DL2(2*nres))
2647 allocate (DDU2(2*nres))
2648 allocate (xsolv(2*nres))
2649 allocate (DML(2*nres))
2650 allocate (rs(2*nres))
2652 allocate (rsold(2*nres))
2653 allocate (matold(2*nres,2*nres))
2655 matold(1,1)=DMorig(1)
2656 matold(1,2)=DU1orig(1)
2657 matold(1,3)=DU2orig(1)
2658 write (*,*) DMorig(1),DU1orig(1),DU2orig(1)
2663 matold(i,j)=DMorig(i)
2664 matold(i,j-1)=DU1orig(i-1)
2666 matold(i,j-2)=DU2orig(i-2)
2674 matold(i,j+1)=DU1orig(i)
2680 matold(i,j+2)=DU2orig(i)
2684 matold(dimen,dimen)=DMorig(dimen)
2685 matold(dimen,dimen-1)=DU1orig(dimen-1)
2686 matold(dimen,dimen-2)=DU2orig(dimen-2)
2687 write(iout,*) "old gmatrix"
2688 call matout(dimen,dimen,2*nres,2*nres,matold)
2692 ! Find the ith eigenvector of the pentadiagonal inertiq matrix
2696 DML(j)=DMorig(j)-geigen(i)
2699 DML(j-1)=DMorig(j)-geigen(i)
2704 DDU1(imark-1)=DU2orig(imark-1)
2705 do j=imark+1,dimen-1
2706 DDU1(j-1)=DU1orig(j)
2714 DDU2(j)=DU2orig(j+1)
2723 write (iout,*) "DL2,DL1,DML,DDU1,DDU2"
2724 write(iout,'(10f10.5)') (DL2(k),k=3,dimen-1)
2725 write(iout,'(10f10.5)') (DL1(k),k=2,dimen-1)
2726 write(iout,'(10f10.5)') (DML(k),k=1,dimen-1)
2727 write(iout,'(10f10.5)') (DDU1(k),k=1,dimen-2)
2728 write(iout,'(10f10.5)') (DDU2(k),k=1,dimen-3)
2731 if (imark.gt.2) rs(imark-2)=-DU2orig(imark-2)
2732 if (imark.gt.1) rs(imark-1)=-DU1orig(imark-1)
2733 if (imark.lt.dimen) rs(imark)=-DU1orig(imark)
2734 if (imark.lt.dimen-1) rs(imark+1)=-DU2orig(imark)
2738 ! write (iout,*) "Vector rs"
2740 ! write (iout,*) j,rs(j)
2743 call FDIAG(dimen-1,DL2,DL1,DML,DDU1,DDU2,rs,xsolv,mark)
2750 sumx=-geigen(i)*xsolv(j)
2752 sumx=sumx+matold(j,k)*xsolv(k)
2755 sumx=sumx+matold(j,k)*xsolv(k-1)
2757 write(iout,'(i5,3f10.5)') j,sumx,rsold(j),sumx-rsold(j)
2760 sumx=-geigen(i)*xsolv(j-1)
2762 sumx=sumx+matold(j,k)*xsolv(k)
2765 sumx=sumx+matold(j,k)*xsolv(k-1)
2767 write(iout,'(i5,3f10.5)') j-1,sumx,rsold(j-1),sumx-rsold(j-1)
2771 "Solution of equations system",i," for eigenvalue",geigen(i)
2773 write(iout,'(i5,f10.5)') j,xsolv(j)
2776 do j=dimen-1,imark,-1
2781 write (iout,*) "Un-normalized eigenvector",i," for eigenvalue",geigen(i)
2783 write(iout,'(i5,f10.5)') j,xsolv(j)
2786 ! Normalize ith eigenvector
2789 sumx=sumx+xsolv(j)**2
2793 xsolv(j)=xsolv(j)/sumx
2796 write (iout,*) "Eigenvector",i," for eigenvalue",geigen(i)
2798 write(iout,'(i5,3f10.5)') j,xsolv(j)
2801 ! All done at this point for eigenvector i, exit loop
2809 write (iout,*) "Unable to find eigenvector",i
2812 ! write (iout,*) "i=",i
2814 ! write (iout,*) "k=",k
2817 ! write(iout,*) "ind",ind," ind1",3*(i-1)+k
2818 d_t_work(ind)=d_t_work(ind) &
2819 +xsolv(j)*d_t_work_new(3*(i-1)+k)
2822 enddo ! i (loop over the eigenvectors)
2825 write (iout,*) "d_t_work"
2827 write (iout,"(i5,f10.5)") i,d_t_work(i)
2832 if (itype(i).eq.10) then
2839 ! write (iout,*) "k",k," ii+k",ii+k," ii+j+k",ii+j+k,"EK1",Ek1
2840 Ek1=Ek1+0.5d0*mp*((d_t_work(ii+j+k)+d_t_work_new(ii+k))/2)**2+&
2841 0.5d0*0.25d0*IP*(d_t_work(ii+j+k)-d_t_work_new(ii+k))**2
2844 if (itype(i).ne.10) ii=ii+3
2845 write (iout,*) "i",i," itype",itype(i)," mass",msc(itype(i))
2846 write (iout,*) "ii",ii
2849 write (iout,*) "k",k," ii",ii,"EK1",EK1
2850 if (iabs(itype(i)).ne.10) Ek1=Ek1+0.5d0*Isc(iabs(itype(i)))*(d_t_work(ii)-d_t_work(ii-3))**2
2851 Ek1=Ek1+0.5d0*msc(iabs(itype(i)))*d_t_work(ii)**2
2853 write (iout,*) "i",i," ii",ii
2855 write (iout,*) "Ek from d_t_work",Ek1
2856 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
2872 d_t(k,j)=d_t_work(ind)
2875 if (itype(j).ne.10 .and. itype(j).ne.ntyp1) then
2877 d_t(k,j+nres)=d_t_work(ind)
2883 write (iout,*) "Random velocities in the Calpha,SC space"
2885 write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3)
2888 write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3)
2895 if (itype(i).eq.10) then
2897 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
2901 d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
2902 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
2907 write (iout,*)"Random velocities in the virtual-bond-vector space"
2909 write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3)
2912 write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3)
2915 write (iout,*) "Ek from d_t_work",Ek1
2916 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
2924 d_t_work(ind)=d_t_work(ind) &
2925 +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
2927 ! write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
2931 ! Transfer to the d_t vector
2933 d_t(j,0)=d_t_work(j)
2939 d_t(j,i)=d_t_work(ind)
2943 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2946 d_t(j,i+nres)=d_t_work(ind)
2952 ! write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",&
2953 ! 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
2957 end subroutine random_vel
2958 !-----------------------------------------------------------------------------
2960 subroutine sd_verlet_p_setup
2961 ! Sets up the parameters of stochastic Verlet algorithm
2962 ! implicit real*8 (a-h,o-z)
2963 ! include 'DIMENSIONS'
2964 use control, only: tcpu
2969 ! include 'COMMON.CONTROL'
2970 ! include 'COMMON.VAR'
2971 ! include 'COMMON.MD'
2973 ! include 'COMMON.LANGEVIN'
2975 ! include 'COMMON.LANGEVIN.lang0'
2977 ! include 'COMMON.CHAIN'
2978 ! include 'COMMON.DERIV'
2979 ! include 'COMMON.GEO'
2980 ! include 'COMMON.LOCAL'
2981 ! include 'COMMON.INTERACT'
2982 ! include 'COMMON.IOUNITS'
2983 ! include 'COMMON.NAMES'
2984 ! include 'COMMON.TIME1'
2985 real(kind=8),dimension(6*nres) :: emgdt !(MAXRES6) maxres6=6*maxres
2986 real(kind=8) :: pterm,vterm,rho,rhoc,vsig
2987 real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
2988 prand_vec,vrand_vec1,vrand_vec2 !(MAXRES6) maxres6=6*maxres
2989 logical :: lprn = .false.
2990 real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
2991 real(kind=8) :: ktm,gdt,egdt,gdt2,gdt3,gdt4,gdt5,gdt6,gdt7,gdt8,&
2993 integer :: i,maxres2
3000 ! AL 8/17/04 Code adapted from tinker
3002 ! Get the frictional and random terms for stochastic dynamics in the
3003 ! eigenspace of mass-scaled UNRES friction matrix
3007 gdt = fricgam(i) * d_time
3009 ! Stochastic dynamics reduces to simple MD for zero friction
3011 if (gdt .le. zero) then
3012 pfric_vec(i) = 1.0d0
3013 vfric_vec(i) = d_time
3014 afric_vec(i) = 0.5d0 * d_time * d_time
3015 prand_vec(i) = 0.0d0
3016 vrand_vec1(i) = 0.0d0
3017 vrand_vec2(i) = 0.0d0
3019 ! Analytical expressions when friction coefficient is large
3022 if (gdt .ge. gdt_radius) then
3025 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
3026 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
3027 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
3028 vterm = 1.0d0 - egdt**2
3029 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
3031 ! Use series expansions when friction coefficient is small
3042 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0 &
3043 - gdt5/120.0d0 + gdt6/720.0d0 &
3044 - gdt7/5040.0d0 + gdt8/40320.0d0 &
3045 - gdt9/362880.0d0) / fricgam(i)**2
3046 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
3047 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
3048 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0 &
3049 + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0 &
3050 + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0 &
3051 + 127.0d0*gdt9/90720.0d0
3052 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0 &
3053 - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0 &
3054 - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0 &
3055 - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
3056 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0 &
3057 - 17.0d0*gdt2/1280.0d0 &
3058 + 17.0d0*gdt3/6144.0d0 &
3059 + 40967.0d0*gdt4/34406400.0d0 &
3060 - 57203.0d0*gdt5/275251200.0d0 &
3061 - 1429487.0d0*gdt6/13212057600.0d0)
3064 ! Compute the scaling factors of random terms for the nonzero friction case
3066 ktm = 0.5d0*d_time/fricgam(i)
3067 psig = dsqrt(ktm*pterm) / fricgam(i)
3068 vsig = dsqrt(ktm*vterm)
3069 rhoc = dsqrt(1.0d0 - rho*rho)
3071 vrand_vec1(i) = vsig * rho
3072 vrand_vec2(i) = vsig * rhoc
3077 "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",&
3080 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),&
3081 afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
3085 ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
3088 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
3089 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
3090 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
3091 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
3092 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
3093 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
3096 t_sdsetup=t_sdsetup+MPI_Wtime()
3098 t_sdsetup=t_sdsetup+tcpu()-tt0
3101 end subroutine sd_verlet_p_setup
3102 !-----------------------------------------------------------------------------
3103 subroutine eigtransf1(n,ndim,ab,d,c)
3107 real(kind=8) :: ab(ndim,ndim,n),c(ndim,n),d(ndim)
3113 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
3118 end subroutine eigtransf1
3119 !-----------------------------------------------------------------------------
3120 subroutine eigtransf(n,ndim,a,b,d,c)
3124 real(kind=8) :: a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
3130 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
3135 end subroutine eigtransf
3136 !-----------------------------------------------------------------------------
3137 subroutine sd_verlet1
3139 ! Applying stochastic velocity Verlet algorithm - step 1 to velocities
3141 ! implicit real*8 (a-h,o-z)
3142 ! include 'DIMENSIONS'
3143 ! include 'COMMON.CONTROL'
3144 ! include 'COMMON.VAR'
3145 ! include 'COMMON.MD'
3147 ! include 'COMMON.LANGEVIN'
3149 ! include 'COMMON.LANGEVIN.lang0'
3151 ! include 'COMMON.CHAIN'
3152 ! include 'COMMON.DERIV'
3153 ! include 'COMMON.GEO'
3154 ! include 'COMMON.LOCAL'
3155 ! include 'COMMON.INTERACT'
3156 ! include 'COMMON.IOUNITS'
3157 ! include 'COMMON.NAMES'
3158 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
3159 !el common /stochcalc/ stochforcvec
3160 logical :: lprn = .false.
3161 real(kind=8) :: ddt1,ddt2
3162 integer :: i,j,ind,inres
3164 ! write (iout,*) "dc_old"
3166 ! write (iout,'(i5,3f10.5,5x,3f10.5)')
3167 ! & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3170 dc_work(j)=dc_old(j,0)
3171 d_t_work(j)=d_t_old(j,0)
3172 d_a_work(j)=d_a_old(j,0)
3177 dc_work(ind+j)=dc_old(j,i)
3178 d_t_work(ind+j)=d_t_old(j,i)
3179 d_a_work(ind+j)=d_a_old(j,i)
3184 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3186 dc_work(ind+j)=dc_old(j,i+nres)
3187 d_t_work(ind+j)=d_t_old(j,i+nres)
3188 d_a_work(ind+j)=d_a_old(j,i+nres)
3196 "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",&
3200 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
3201 vfric_mat(i,j),afric_mat(i,j),&
3202 prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3210 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) &
3211 +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3212 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3213 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3215 d_t_work_new(i)=ddt1+0.5d0*ddt2
3216 d_t_work(i)=ddt1+ddt2
3221 d_t(j,0)=d_t_work(j)
3226 dc(j,i)=dc_work(ind+j)
3227 d_t(j,i)=d_t_work(ind+j)
3232 if (itype(i).ne.10) then
3235 dc(j,inres)=dc_work(ind+j)
3236 d_t(j,inres)=d_t_work(ind+j)
3242 end subroutine sd_verlet1
3243 !-----------------------------------------------------------------------------
3244 subroutine sd_verlet2
3246 ! Calculating the adjusted velocities for accelerations
3248 ! implicit real*8 (a-h,o-z)
3249 ! include 'DIMENSIONS'
3250 ! include 'COMMON.CONTROL'
3251 ! include 'COMMON.VAR'
3252 ! include 'COMMON.MD'
3254 ! include 'COMMON.LANGEVIN'
3256 ! include 'COMMON.LANGEVIN.lang0'
3258 ! include 'COMMON.CHAIN'
3259 ! include 'COMMON.DERIV'
3260 ! include 'COMMON.GEO'
3261 ! include 'COMMON.LOCAL'
3262 ! include 'COMMON.INTERACT'
3263 ! include 'COMMON.IOUNITS'
3264 ! include 'COMMON.NAMES'
3265 !el real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV !(MAXRES6) maxres6=6*maxres
3266 real(kind=8),dimension(6*nres) :: stochforcvecV !(MAXRES6) maxres6=6*maxres
3267 !el common /stochcalc/ stochforcvec
3269 real(kind=8) :: ddt1,ddt2
3270 integer :: i,j,ind,inres
3271 ! Compute the stochastic forces which contribute to velocity change
3273 call stochastic_force(stochforcvecV)
3280 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3281 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+ &
3282 vrand_mat2(i,j)*stochforcvecV(j)
3284 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
3288 d_t(j,0)=d_t_work(j)
3293 d_t(j,i)=d_t_work(ind+j)
3298 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3301 d_t(j,inres)=d_t_work(ind+j)
3307 end subroutine sd_verlet2
3308 !-----------------------------------------------------------------------------
3309 subroutine sd_verlet_ciccotti_setup
3311 ! Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
3313 ! implicit real*8 (a-h,o-z)
3314 ! include 'DIMENSIONS'
3315 use control, only: tcpu
3320 ! include 'COMMON.CONTROL'
3321 ! include 'COMMON.VAR'
3322 ! include 'COMMON.MD'
3324 ! include 'COMMON.LANGEVIN'
3326 ! include 'COMMON.LANGEVIN.lang0'
3328 ! include 'COMMON.CHAIN'
3329 ! include 'COMMON.DERIV'
3330 ! include 'COMMON.GEO'
3331 ! include 'COMMON.LOCAL'
3332 ! include 'COMMON.INTERACT'
3333 ! include 'COMMON.IOUNITS'
3334 ! include 'COMMON.NAMES'
3335 ! include 'COMMON.TIME1'
3336 real(kind=8),dimension(6*nres) :: emgdt !(MAXRES6) maxres6=6*maxres
3337 real(kind=8) :: pterm,vterm,rho,rhoc,vsig
3338 real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
3339 prand_vec,vrand_vec1,vrand_vec2 !(MAXRES6) maxres6=6*maxres
3340 logical :: lprn = .false.
3341 real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
3342 real(kind=8) :: ktm,gdt,egdt,tt0
3343 integer :: i,maxres2
3350 ! AL 8/17/04 Code adapted from tinker
3352 ! Get the frictional and random terms for stochastic dynamics in the
3353 ! eigenspace of mass-scaled UNRES friction matrix
3357 write (iout,*) "i",i," fricgam",fricgam(i)
3358 gdt = fricgam(i) * d_time
3360 ! Stochastic dynamics reduces to simple MD for zero friction
3362 if (gdt .le. zero) then
3363 pfric_vec(i) = 1.0d0
3364 vfric_vec(i) = d_time
3365 afric_vec(i) = 0.5d0*d_time*d_time
3366 prand_vec(i) = afric_vec(i)
3367 vrand_vec2(i) = vfric_vec(i)
3369 ! Analytical expressions when friction coefficient is large
3374 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
3375 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
3376 prand_vec(i) = afric_vec(i)
3377 vrand_vec2(i) = vfric_vec(i)
3379 ! Compute the scaling factors of random terms for the nonzero friction case
3381 ! ktm = 0.5d0*d_time/fricgam(i)
3382 ! psig = dsqrt(ktm*pterm) / fricgam(i)
3383 ! vsig = dsqrt(ktm*vterm)
3384 ! prand_vec(i) = psig*afric_vec(i)
3385 ! vrand_vec2(i) = vsig*vfric_vec(i)
3390 "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",&
3393 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),&
3394 afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
3398 ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
3400 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
3401 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
3402 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
3403 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
3404 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
3406 t_sdsetup=t_sdsetup+MPI_Wtime()
3408 t_sdsetup=t_sdsetup+tcpu()-tt0
3411 end subroutine sd_verlet_ciccotti_setup
3412 !-----------------------------------------------------------------------------
3413 subroutine sd_verlet1_ciccotti
3415 ! Applying stochastic velocity Verlet algorithm - step 1 to velocities
3416 ! implicit real*8 (a-h,o-z)
3418 ! include 'DIMENSIONS'
3422 ! include 'COMMON.CONTROL'
3423 ! include 'COMMON.VAR'
3424 ! include 'COMMON.MD'
3426 ! include 'COMMON.LANGEVIN'
3428 ! include 'COMMON.LANGEVIN.lang0'
3430 ! include 'COMMON.CHAIN'
3431 ! include 'COMMON.DERIV'
3432 ! include 'COMMON.GEO'
3433 ! include 'COMMON.LOCAL'
3434 ! include 'COMMON.INTERACT'
3435 ! include 'COMMON.IOUNITS'
3436 ! include 'COMMON.NAMES'
3437 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
3438 !el common /stochcalc/ stochforcvec
3439 logical :: lprn = .false.
3440 real(kind=8) :: ddt1,ddt2
3441 integer :: i,j,ind,inres
3442 ! write (iout,*) "dc_old"
3444 ! write (iout,'(i5,3f10.5,5x,3f10.5)')
3445 ! & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3448 dc_work(j)=dc_old(j,0)
3449 d_t_work(j)=d_t_old(j,0)
3450 d_a_work(j)=d_a_old(j,0)
3455 dc_work(ind+j)=dc_old(j,i)
3456 d_t_work(ind+j)=d_t_old(j,i)
3457 d_a_work(ind+j)=d_a_old(j,i)
3462 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3464 dc_work(ind+j)=dc_old(j,i+nres)
3465 d_t_work(ind+j)=d_t_old(j,i+nres)
3466 d_a_work(ind+j)=d_a_old(j,i+nres)
3475 "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",&
3479 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
3480 vfric_mat(i,j),afric_mat(i,j),&
3481 prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3489 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) &
3490 +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3491 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3492 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3494 d_t_work_new(i)=ddt1+0.5d0*ddt2
3495 d_t_work(i)=ddt1+ddt2
3500 d_t(j,0)=d_t_work(j)
3505 dc(j,i)=dc_work(ind+j)
3506 d_t(j,i)=d_t_work(ind+j)
3511 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3514 dc(j,inres)=dc_work(ind+j)
3515 d_t(j,inres)=d_t_work(ind+j)
3521 end subroutine sd_verlet1_ciccotti
3522 !-----------------------------------------------------------------------------
3523 subroutine sd_verlet2_ciccotti
3525 ! Calculating the adjusted velocities for accelerations
3527 ! implicit real*8 (a-h,o-z)
3528 ! include 'DIMENSIONS'
3529 ! include 'COMMON.CONTROL'
3530 ! include 'COMMON.VAR'
3531 ! include 'COMMON.MD'
3533 ! include 'COMMON.LANGEVIN'
3535 ! include 'COMMON.LANGEVIN.lang0'
3537 ! include 'COMMON.CHAIN'
3538 ! include 'COMMON.DERIV'
3539 ! include 'COMMON.GEO'
3540 ! include 'COMMON.LOCAL'
3541 ! include 'COMMON.INTERACT'
3542 ! include 'COMMON.IOUNITS'
3543 ! include 'COMMON.NAMES'
3544 !el real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV !(MAXRES6) maxres6=6*maxres
3545 real(kind=8),dimension(6*nres) :: stochforcvecV !(MAXRES6) maxres6=6*maxres
3546 !el common /stochcalc/ stochforcvec
3547 real(kind=8) :: ddt1,ddt2
3548 integer :: i,j,ind,inres
3550 ! Compute the stochastic forces which contribute to velocity change
3552 call stochastic_force(stochforcvecV)
3559 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3560 ! ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
3561 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
3563 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
3567 d_t(j,0)=d_t_work(j)
3572 d_t(j,i)=d_t_work(ind+j)
3577 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3580 d_t(j,inres)=d_t_work(ind+j)
3586 end subroutine sd_verlet2_ciccotti
3588 !-----------------------------------------------------------------------------
3590 !-----------------------------------------------------------------------------
3591 subroutine inertia_tensor
3593 ! Calculating the intertia tensor for the entire protein in order to
3594 ! remove the perpendicular components of velocity matrix which cause
3595 ! the molecule to rotate.
3598 ! implicit real*8 (a-h,o-z)
3599 ! include 'DIMENSIONS'
3600 ! include 'COMMON.CONTROL'
3601 ! include 'COMMON.VAR'
3602 ! include 'COMMON.MD'
3603 ! include 'COMMON.CHAIN'
3604 ! include 'COMMON.DERIV'
3605 ! include 'COMMON.GEO'
3606 ! include 'COMMON.LOCAL'
3607 ! include 'COMMON.INTERACT'
3608 ! include 'COMMON.IOUNITS'
3609 ! include 'COMMON.NAMES'
3611 real(kind=8),dimension(3,3) :: Im,Imcp,eigvec,Id
3612 real(kind=8),dimension(3) :: pr,eigval,L,vp,vrot
3613 real(kind=8) :: M_SC,mag,mag2
3614 real(kind=8),dimension(3,0:nres) :: vpp !(3,0:MAXRES)
3615 real(kind=8),dimension(3) :: vs_p,pp,incr,v
3616 real(kind=8),dimension(3,3) :: pr1,pr2
3618 !el common /gucio/ cm
3619 integer :: iti,inres,i,j,k
3630 ! calculating the center of the mass of the protein
3633 cm(j)=cm(j)+c(j,i)+0.5d0*dc(j,i)
3642 M_SC=M_SC+msc(iabs(iti))
3645 cm(j)=cm(j)+msc(iabs(iti))*c(j,inres)
3649 cm(j)=cm(j)/(M_SC+(nct-nnt)*mp)
3654 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
3656 Im(1,1)=Im(1,1)+mp*(pr(2)*pr(2)+pr(3)*pr(3))
3657 Im(1,2)=Im(1,2)-mp*pr(1)*pr(2)
3658 Im(1,3)=Im(1,3)-mp*pr(1)*pr(3)
3659 Im(2,3)=Im(2,3)-mp*pr(2)*pr(3)
3660 Im(2,2)=Im(2,2)+mp*(pr(3)*pr(3)+pr(1)*pr(1))
3661 Im(3,3)=Im(3,3)+mp*(pr(1)*pr(1)+pr(2)*pr(2))
3668 pr(j)=c(j,inres)-cm(j)
3670 Im(1,1)=Im(1,1)+msc(iabs(iti))*(pr(2)*pr(2)+pr(3)*pr(3))
3671 Im(1,2)=Im(1,2)-msc(iabs(iti))*pr(1)*pr(2)
3672 Im(1,3)=Im(1,3)-msc(iabs(iti))*pr(1)*pr(3)
3673 Im(2,3)=Im(2,3)-msc(iabs(iti))*pr(2)*pr(3)
3674 Im(2,2)=Im(2,2)+msc(iabs(iti))*(pr(3)*pr(3)+pr(1)*pr(1))
3675 Im(3,3)=Im(3,3)+msc(iabs(iti))*(pr(1)*pr(1)+pr(2)*pr(2))
3679 Im(1,1)=Im(1,1)+Ip*(1-dc_norm(1,i)*dc_norm(1,i))* &
3680 vbld(i+1)*vbld(i+1)*0.25d0
3681 Im(1,2)=Im(1,2)+Ip*(-dc_norm(1,i)*dc_norm(2,i))* &
3682 vbld(i+1)*vbld(i+1)*0.25d0
3683 Im(1,3)=Im(1,3)+Ip*(-dc_norm(1,i)*dc_norm(3,i))* &
3684 vbld(i+1)*vbld(i+1)*0.25d0
3685 Im(2,3)=Im(2,3)+Ip*(-dc_norm(2,i)*dc_norm(3,i))* &
3686 vbld(i+1)*vbld(i+1)*0.25d0
3687 Im(2,2)=Im(2,2)+Ip*(1-dc_norm(2,i)*dc_norm(2,i))* &
3688 vbld(i+1)*vbld(i+1)*0.25d0
3689 Im(3,3)=Im(3,3)+Ip*(1-dc_norm(3,i)*dc_norm(3,i))* &
3690 vbld(i+1)*vbld(i+1)*0.25d0
3695 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3698 Im(1,1)=Im(1,1)+Isc(iti)*(1-dc_norm(1,inres)* &
3699 dc_norm(1,inres))*vbld(inres)*vbld(inres)
3700 Im(1,2)=Im(1,2)-Isc(iti)*(dc_norm(1,inres)* &
3701 dc_norm(2,inres))*vbld(inres)*vbld(inres)
3702 Im(1,3)=Im(1,3)-Isc(iti)*(dc_norm(1,inres)* &
3703 dc_norm(3,inres))*vbld(inres)*vbld(inres)
3704 Im(2,3)=Im(2,3)-Isc(iti)*(dc_norm(2,inres)* &
3705 dc_norm(3,inres))*vbld(inres)*vbld(inres)
3706 Im(2,2)=Im(2,2)+Isc(iti)*(1-dc_norm(2,inres)* &
3707 dc_norm(2,inres))*vbld(inres)*vbld(inres)
3708 Im(3,3)=Im(3,3)+Isc(iti)*(1-dc_norm(3,inres)* &
3709 dc_norm(3,inres))*vbld(inres)*vbld(inres)
3714 ! write(iout,*) "The angular momentum before adjustment:"
3715 ! write(iout,*) (L(j),j=1,3)
3721 ! Copying the Im matrix for the djacob subroutine
3729 ! Finding the eigenvectors and eignvalues of the inertia tensor
3730 call djacob(3,3,10000,1.0d-10,Imcp,eigvec,eigval)
3731 ! write (iout,*) "Eigenvalues & Eigenvectors"
3732 ! write (iout,'(5x,3f10.5)') (eigval(i),i=1,3)
3735 ! write (iout,'(i5,3f10.5)') i,(eigvec(i,j),j=1,3)
3737 ! Constructing the diagonalized matrix
3739 if (dabs(eigval(i)).gt.1.0d-15) then
3740 Id(i,i)=1.0d0/eigval(i)
3747 Imcp(i,j)=eigvec(j,i)
3753 pr1(i,j)=pr1(i,j)+Id(i,k)*Imcp(k,j)
3760 pr2(i,j)=pr2(i,j)+eigvec(i,k)*pr1(k,j)
3764 ! Calculating the total rotational velocity of the molecule
3767 vrot(i)=vrot(i)+pr2(i,j)*L(j)
3770 ! Resetting the velocities
3772 call vecpr(vrot(1),dc(1,i),vp)
3774 d_t(j,i)=d_t(j,i)-vp(j)
3778 if(itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3780 call vecpr(vrot(1),dc(1,inres),vp)
3782 d_t(j,inres)=d_t(j,inres)-vp(j)
3787 ! write(iout,*) "The angular momentum after adjustment:"
3788 ! write(iout,*) (L(j),j=1,3)
3791 end subroutine inertia_tensor
3792 !-----------------------------------------------------------------------------
3793 subroutine angmom(cm,L)
3796 ! implicit real*8 (a-h,o-z)
3797 ! include 'DIMENSIONS'
3798 ! include 'COMMON.CONTROL'
3799 ! include 'COMMON.VAR'
3800 ! include 'COMMON.MD'
3801 ! include 'COMMON.CHAIN'
3802 ! include 'COMMON.DERIV'
3803 ! include 'COMMON.GEO'
3804 ! include 'COMMON.LOCAL'
3805 ! include 'COMMON.INTERACT'
3806 ! include 'COMMON.IOUNITS'
3807 ! include 'COMMON.NAMES'
3809 real(kind=8),dimension(3) :: L,cm,pr,vp,vrot,incr,v,pp
3810 integer :: iti,inres,i,j
3811 ! Calculate the angular momentum
3820 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
3823 v(j)=incr(j)+0.5d0*d_t(j,i)
3826 incr(j)=incr(j)+d_t(j,i)
3828 call vecpr(pr(1),v(1),vp)
3834 pp(j)=0.5d0*d_t(j,i)
3836 call vecpr(pr(1),pp(1),vp)
3848 pr(j)=c(j,inres)-cm(j)
3850 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3852 v(j)=incr(j)+d_t(j,inres)
3859 call vecpr(pr(1),v(1),vp)
3860 ! write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),
3861 ! & " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
3863 L(j)=L(j)+msc(iabs(iti))*vp(j)
3865 ! write (iout,*) "L",(l(j),j=1,3)
3866 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3868 v(j)=incr(j)+d_t(j,inres)
3870 call vecpr(dc(1,inres),d_t(1,inres),vp)
3872 L(j)=L(j)+Isc(iti)*vp(j)
3876 incr(j)=incr(j)+d_t(j,i)
3880 end subroutine angmom
3881 !-----------------------------------------------------------------------------
3882 subroutine vcm_vel(vcm)
3885 ! implicit real*8 (a-h,o-z)
3886 ! include 'DIMENSIONS'
3887 ! include 'COMMON.VAR'
3888 ! include 'COMMON.MD'
3889 ! include 'COMMON.CHAIN'
3890 ! include 'COMMON.DERIV'
3891 ! include 'COMMON.GEO'
3892 ! include 'COMMON.LOCAL'
3893 ! include 'COMMON.INTERACT'
3894 ! include 'COMMON.IOUNITS'
3895 real(kind=8),dimension(3) :: vcm,vv
3896 real(kind=8) :: summas,amas
3908 vcm(j)=vcm(j)+mp*(vv(j)+0.5d0*d_t(j,i))
3911 amas=msc(iabs(itype(i)))
3913 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3915 vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres))
3919 vcm(j)=vcm(j)+amas*vv(j)
3923 vv(j)=vv(j)+d_t(j,i)
3926 ! write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas
3928 vcm(j)=vcm(j)/summas
3931 end subroutine vcm_vel
3932 !-----------------------------------------------------------------------------
3934 !-----------------------------------------------------------------------------
3936 ! RATTLE algorithm for velocity Verlet - step 1, UNRES
3940 ! implicit real*8 (a-h,o-z)
3941 ! include 'DIMENSIONS'
3943 ! include 'COMMON.CONTROL'
3944 ! include 'COMMON.VAR'
3945 ! include 'COMMON.MD'
3947 ! include 'COMMON.LANGEVIN'
3949 ! include 'COMMON.LANGEVIN.lang0'
3951 ! include 'COMMON.CHAIN'
3952 ! include 'COMMON.DERIV'
3953 ! include 'COMMON.GEO'
3954 ! include 'COMMON.LOCAL'
3955 ! include 'COMMON.INTERACT'
3956 ! include 'COMMON.IOUNITS'
3957 ! include 'COMMON.NAMES'
3958 ! include 'COMMON.TIME1'
3959 !el real(kind=8) :: gginv(2*nres,2*nres),&
3960 !el gdc(3,2*nres,2*nres)
3961 real(kind=8) :: dC_uncor(3,2*nres) !,&
3962 !el real(kind=8) :: Cmat(2*nres,2*nres)
3963 real(kind=8) :: x(2*nres),xcorr(3,2*nres) !maxres2=2*maxres
3964 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
3965 !el common /przechowalnia/ nbond
3966 integer :: max_rattle = 5
3967 logical :: lprn = .false., lprn1 = .false., not_done
3968 real(kind=8) :: tol_rattle = 1.0d-5
3970 integer :: ii,i,j,jj,l,ind,ind1,nres2
3973 !el /common/ przechowalnia
3975 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
3976 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
3977 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
3979 if (lprn) write (iout,*) "RATTLE1"
3982 if (itype(i).ne.10) nbond=nbond+1
3984 ! Make a folded form of the Ginv-matrix
3997 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4001 if (itype(k).ne.10) then
4005 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4012 if (itype(i).ne.10) then
4022 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4026 if (itype(k).ne.10) then
4030 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4038 write (iout,*) "Matrix GGinv"
4039 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
4045 if (iter.gt.max_rattle) then
4046 write (iout,*) "Error - too many iterations in RATTLE."
4049 ! Calculate the matrix C = GG**(-1) dC_old o dC
4054 dC_uncor(j,ind1)=dC(j,i)
4058 if (itype(i).ne.10) then
4061 dC_uncor(j,ind1)=dC(j,i+nres)
4070 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
4074 if (itype(k).ne.10) then
4077 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
4082 ! Calculate deviations from standard virtual-bond lengths
4086 x(ind)=vbld(i+1)**2-vbl**2
4089 if (itype(i).ne.10) then
4091 x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i))**2
4095 write (iout,*) "Coordinates and violations"
4097 write(iout,'(i5,3f10.5,5x,e15.5)') &
4098 i,(dC_uncor(j,i),j=1,3),x(i)
4100 write (iout,*) "Velocities and violations"
4104 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4105 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4108 if (itype(i).ne.10) then
4110 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4111 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4112 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4115 ! write (iout,*) "gdc"
4117 ! write (iout,*) "i",i
4119 ! write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4125 if (dabs(x(i)).gt.xmax) then
4129 if (xmax.lt.tol_rattle) then
4133 ! Calculate the matrix of the system of equations
4138 Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
4143 write (iout,*) "Matrix Cmat"
4144 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4146 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4147 ! Add constraint term to positions
4154 xx = xx+x(ii)*gdc(j,ind,ii)
4158 d_t_new(j,i)=d_t_new(j,i)-xx/d_time
4162 if (itype(i).ne.10) then
4167 xx = xx+x(ii)*gdc(j,ind,ii)
4170 dC(j,i+nres)=dC(j,i+nres)-xx
4171 d_t_new(j,i+nres)=d_t_new(j,i+nres)-xx/d_time
4175 ! Rebuild the chain using the new coordinates
4176 call chainbuild_cart
4178 write (iout,*) "New coordinates, Lagrange multipliers,",&
4179 " and differences between actual and standard bond lengths"
4183 xx=vbld(i+1)**2-vbl**2
4184 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4185 i,(dC(j,i),j=1,3),x(ind),xx
4188 if (itype(i).ne.10) then
4190 xx=vbld(i+nres)**2-vbldsc0(1,itype(i))**2
4191 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4192 i,(dC(j,i+nres),j=1,3),x(ind),xx
4195 write (iout,*) "Velocities and violations"
4199 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4200 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4203 if (itype(i).ne.10) then
4205 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4206 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4207 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4214 10 write (iout,*) "Error - singularity in solving the system",&
4215 " of equations for Lagrange multipliers."
4219 "RATTLE inactive; use -DRATTLE switch at compile time."
4222 end subroutine rattle1
4223 !-----------------------------------------------------------------------------
4225 ! RATTLE algorithm for velocity Verlet - step 2, UNRES
4229 ! implicit real*8 (a-h,o-z)
4230 ! include 'DIMENSIONS'
4232 ! include 'COMMON.CONTROL'
4233 ! include 'COMMON.VAR'
4234 ! include 'COMMON.MD'
4236 ! include 'COMMON.LANGEVIN'
4238 ! include 'COMMON.LANGEVIN.lang0'
4240 ! include 'COMMON.CHAIN'
4241 ! include 'COMMON.DERIV'
4242 ! include 'COMMON.GEO'
4243 ! include 'COMMON.LOCAL'
4244 ! include 'COMMON.INTERACT'
4245 ! include 'COMMON.IOUNITS'
4246 ! include 'COMMON.NAMES'
4247 ! include 'COMMON.TIME1'
4248 !el real(kind=8) :: gginv(2*nres,2*nres),&
4249 !el gdc(3,2*nres,2*nres)
4250 real(kind=8) :: dC_uncor(3,2*nres) !,&
4251 !el Cmat(2*nres,2*nres)
4252 real(kind=8) :: x(2*nres) !maxres2=2*maxres
4253 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
4254 !el common /przechowalnia/ nbond
4255 integer :: max_rattle = 5
4256 logical :: lprn = .false., lprn1 = .false., not_done
4257 real(kind=8) :: tol_rattle = 1.0d-5
4261 !el /common/ przechowalnia
4262 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4263 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4264 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4266 if (lprn) write (iout,*) "RATTLE2"
4267 if (lprn) write (iout,*) "Velocity correction"
4268 ! Calculate the matrix G dC
4274 gdc(j,i,ind)=GGinv(i,ind)*dC(j,k)
4278 if (itype(k).ne.10) then
4281 gdc(j,i,ind)=GGinv(i,ind)*dC(j,k+nres)
4287 ! write (iout,*) "gdc"
4289 ! write (iout,*) "i",i
4291 ! write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4295 ! Calculate the matrix of the system of equations
4302 Cmat(ind,j)=Cmat(ind,j)+dC(k,i)*gdc(k,ind,j)
4307 if (itype(i).ne.10) then
4312 Cmat(ind,j)=Cmat(ind,j)+dC(k,i+nres)*gdc(k,ind,j)
4317 ! Calculate the scalar product dC o d_t_new
4321 x(ind)=scalar(d_t(1,i),dC(1,i))
4324 if (itype(i).ne.10) then
4326 x(ind)=scalar(d_t(1,i+nres),dC(1,i+nres))
4330 write (iout,*) "Velocities and violations"
4334 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4335 i,ind,(d_t(j,i),j=1,3),x(ind)
4338 if (itype(i).ne.10) then
4340 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4341 i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind)
4347 if (dabs(x(i)).gt.xmax) then
4351 if (xmax.lt.tol_rattle) then
4356 write (iout,*) "Matrix Cmat"
4357 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4359 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4360 ! Add constraint term to velocities
4367 xx = xx+x(ii)*gdc(j,ind,ii)
4369 d_t(j,i)=d_t(j,i)-xx
4373 if (itype(i).ne.10) then
4378 xx = xx+x(ii)*gdc(j,ind,ii)
4380 d_t(j,i+nres)=d_t(j,i+nres)-xx
4386 "New velocities, Lagrange multipliers violations"
4390 if (lprn) write (iout,'(2i5,3f10.5,5x,2e15.5)') &
4391 i,ind,(d_t(j,i),j=1,3),x(ind),scalar(d_t(1,i),dC(1,i))
4394 if (itype(i).ne.10) then
4396 write (iout,'(2i5,3f10.5,5x,2e15.5)') &
4397 i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind),&
4398 scalar(d_t(1,i+nres),dC(1,i+nres))
4404 10 write (iout,*) "Error - singularity in solving the system",&
4405 " of equations for Lagrange multipliers."
4409 "RATTLE inactive; use -DRATTLE option at compile time."
4412 end subroutine rattle2
4413 !-----------------------------------------------------------------------------
4414 subroutine rattle_brown
4415 ! RATTLE/LINCS algorithm for Brownian dynamics, UNRES
4419 ! implicit real*8 (a-h,o-z)
4420 ! include 'DIMENSIONS'
4422 ! include 'COMMON.CONTROL'
4423 ! include 'COMMON.VAR'
4424 ! include 'COMMON.MD'
4426 ! include 'COMMON.LANGEVIN'
4428 ! include 'COMMON.LANGEVIN.lang0'
4430 ! include 'COMMON.CHAIN'
4431 ! include 'COMMON.DERIV'
4432 ! include 'COMMON.GEO'
4433 ! include 'COMMON.LOCAL'
4434 ! include 'COMMON.INTERACT'
4435 ! include 'COMMON.IOUNITS'
4436 ! include 'COMMON.NAMES'
4437 ! include 'COMMON.TIME1'
4438 !el real(kind=8) :: gginv(2*nres,2*nres),&
4439 !el gdc(3,2*nres,2*nres)
4440 real(kind=8) :: dC_uncor(3,2*nres) !,&
4441 !el real(kind=8) :: Cmat(2*nres,2*nres)
4442 real(kind=8) :: x(2*nres) !maxres2=2*maxres
4443 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
4444 !el common /przechowalnia/ nbond
4445 integer :: max_rattle = 5
4446 logical :: lprn = .true., lprn1 = .true., not_done
4447 real(kind=8) :: tol_rattle = 1.0d-5
4451 !el /common/ przechowalnia
4452 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4453 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4454 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4457 if (lprn) write (iout,*) "RATTLE_BROWN"
4460 if (itype(i).ne.10) nbond=nbond+1
4462 ! Make a folded form of the Ginv-matrix
4475 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4479 if (itype(k).ne.10) then
4483 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4490 if (itype(i).ne.10) then
4500 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4504 if (itype(k).ne.10) then
4508 if (j.eq.1 .and. l.eq.1)GGinv(ii,jj)=fricmat(ind,ind1)
4516 write (iout,*) "Matrix GGinv"
4517 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
4523 if (iter.gt.max_rattle) then
4524 write (iout,*) "Error - too many iterations in RATTLE."
4527 ! Calculate the matrix C = GG**(-1) dC_old o dC
4532 dC_uncor(j,ind1)=dC(j,i)
4536 if (itype(i).ne.10) then
4539 dC_uncor(j,ind1)=dC(j,i+nres)
4548 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
4552 if (itype(k).ne.10) then
4555 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
4560 ! Calculate deviations from standard virtual-bond lengths
4564 x(ind)=vbld(i+1)**2-vbl**2
4567 if (itype(i).ne.10) then
4569 x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i))**2
4573 write (iout,*) "Coordinates and violations"
4575 write(iout,'(i5,3f10.5,5x,e15.5)') &
4576 i,(dC_uncor(j,i),j=1,3),x(i)
4578 write (iout,*) "Velocities and violations"
4582 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4583 i,ind,(d_t(j,i),j=1,3),scalar(d_t(1,i),dC_old(1,i))
4586 if (itype(i).ne.10) then
4588 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4589 i+nres,ind,(d_t(j,i+nres),j=1,3),&
4590 scalar(d_t(1,i+nres),dC_old(1,i+nres))
4593 write (iout,*) "gdc"
4595 write (iout,*) "i",i
4597 write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4603 if (dabs(x(i)).gt.xmax) then
4607 if (xmax.lt.tol_rattle) then
4611 ! Calculate the matrix of the system of equations
4616 Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
4621 write (iout,*) "Matrix Cmat"
4622 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4624 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4625 ! Add constraint term to positions
4632 xx = xx+x(ii)*gdc(j,ind,ii)
4635 d_t(j,i)=d_t(j,i)+xx/d_time
4640 if (itype(i).ne.10) then
4645 xx = xx+x(ii)*gdc(j,ind,ii)
4648 d_t(j,i+nres)=d_t(j,i+nres)+xx/d_time
4649 dC(j,i+nres)=dC(j,i+nres)+xx
4653 ! Rebuild the chain using the new coordinates
4654 call chainbuild_cart
4656 write (iout,*) "New coordinates, Lagrange multipliers,",&
4657 " and differences between actual and standard bond lengths"
4661 xx=vbld(i+1)**2-vbl**2
4662 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4663 i,(dC(j,i),j=1,3),x(ind),xx
4666 if (itype(i).ne.10) then
4668 xx=vbld(i+nres)**2-vbldsc0(1,itype(i))**2
4669 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4670 i,(dC(j,i+nres),j=1,3),x(ind),xx
4673 write (iout,*) "Velocities and violations"
4677 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4678 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4681 if (itype(i).ne.10) then
4683 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4684 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4685 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4692 10 write (iout,*) "Error - singularity in solving the system",&
4693 " of equations for Lagrange multipliers."
4697 "RATTLE inactive; use -DRATTLE option at compile time"
4700 end subroutine rattle_brown
4701 !-----------------------------------------------------------------------------
4703 !-----------------------------------------------------------------------------
4704 subroutine friction_force
4709 ! implicit real*8 (a-h,o-z)
4710 ! include 'DIMENSIONS'
4711 ! include 'COMMON.VAR'
4712 ! include 'COMMON.CHAIN'
4713 ! include 'COMMON.DERIV'
4714 ! include 'COMMON.GEO'
4715 ! include 'COMMON.LOCAL'
4716 ! include 'COMMON.INTERACT'
4717 ! include 'COMMON.MD'
4719 ! include 'COMMON.LANGEVIN'
4721 ! include 'COMMON.LANGEVIN.lang0'
4723 ! include 'COMMON.IOUNITS'
4724 !el real(kind=8),dimension(6*nres) :: gamvec !(MAXRES6) maxres6=6*maxres
4725 !el common /syfek/ gamvec
4726 real(kind=8) :: vv(3),vvtot(3,nres),v_work(6*nres) !,&
4727 !el ginvfric(2*nres,2*nres) !maxres2=2*maxres
4728 !el common /przechowalnia/ ginvfric
4730 logical :: lprn = .false., checkmode = .false.
4731 integer :: i,j,ind,k,nres2,nres6
4735 if(.not.allocated(gamvec)) allocate(gamvec(nres6)) !(MAXRES6)
4736 if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
4744 d_t_work(j)=d_t(j,0)
4749 d_t_work(ind+j)=d_t(j,i)
4754 if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
4756 d_t_work(ind+j)=d_t(j,i+nres)
4762 call fricmat_mult(d_t_work,fric_work)
4764 if (.not.checkmode) return
4767 write (iout,*) "d_t_work and fric_work"
4769 write (iout,'(i3,2e15.5)') i,d_t_work(i),fric_work(i)
4773 friction(j,0)=fric_work(j)
4778 friction(j,i)=fric_work(ind+j)
4783 if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
4785 friction(j,i+nres)=fric_work(ind+j)
4791 write(iout,*) "Friction backbone"
4793 write(iout,'(i5,3e15.5,5x,3e15.5)') &
4794 i,(friction(j,i),j=1,3),(d_t(j,i),j=1,3)
4796 write(iout,*) "Friction side chain"
4798 write(iout,'(i5,3e15.5,5x,3e15.5)') &
4799 i,(friction(j,i+nres),j=1,3),(d_t(j,i+nres),j=1,3)
4808 vvtot(j,i)=vv(j)+0.5d0*d_t(j,i)
4809 vvtot(j,i+nres)=vv(j)+d_t(j,i+nres)
4810 vv(j)=vv(j)+d_t(j,i)
4813 write (iout,*) "vvtot backbone and sidechain"
4815 write (iout,'(i5,3e15.5,5x,3e15.5)') i,(vvtot(j,i),j=1,3),&
4816 (vvtot(j,i+nres),j=1,3)
4821 v_work(ind+j)=vvtot(j,i)
4827 v_work(ind+j)=vvtot(j,i+nres)
4831 write (iout,*) "v_work gamvec and site-based friction forces"
4833 write (iout,'(i5,3e15.5)') i,v_work(i),gamvec(i),&
4837 ! fric_work1(i)=0.0d0
4839 ! fric_work1(i)=fric_work1(i)-A(j,i)*gamvec(j)*v_work(j)
4842 ! write (iout,*) "fric_work and fric_work1"
4844 ! write (iout,'(i5,2e15.5)') i,fric_work(i),fric_work1(i)
4850 ginvfric(i,j)=ginvfric(i,j)+ginv(i,k)*fricmat(k,j)
4854 write (iout,*) "ginvfric"
4856 write (iout,'(i5,100f8.3)') i,(ginvfric(i,j),j=1,dimen)
4858 write (iout,*) "symmetry check"
4861 write (iout,*) i,j,ginvfric(i,j)-ginvfric(j,i)
4866 end subroutine friction_force
4867 !-----------------------------------------------------------------------------
4868 subroutine setup_fricmat
4872 use control_data, only:time_Bcast
4873 use control, only:tcpu
4875 ! implicit real*8 (a-h,o-z)
4879 real(kind=8) :: time00
4881 ! include 'DIMENSIONS'
4882 ! include 'COMMON.VAR'
4883 ! include 'COMMON.CHAIN'
4884 ! include 'COMMON.DERIV'
4885 ! include 'COMMON.GEO'
4886 ! include 'COMMON.LOCAL'
4887 ! include 'COMMON.INTERACT'
4888 ! include 'COMMON.MD'
4889 ! include 'COMMON.SETUP'
4890 ! include 'COMMON.TIME1'
4891 ! integer licznik /0/
4894 ! include 'COMMON.LANGEVIN'
4896 ! include 'COMMON.LANGEVIN.lang0'
4898 ! include 'COMMON.IOUNITS'
4900 integer :: i,j,ind,ind1,m
4901 logical :: lprn = .false.
4902 real(kind=8) :: dtdi !el ,gamvec(2*nres)
4903 !el real(kind=8),dimension(2*nres,2*nres) :: ginvfric,fcopy
4904 real(kind=8),dimension(2*nres,2*nres) :: fcopy
4905 !el real(kind=8),dimension(2*nres*(2*nres+1)/2) :: Ghalf !(mmaxres2) (mmaxres2=(maxres2*(maxres2+1)/2))
4906 !el common /syfek/ gamvec
4907 real(kind=8) :: work(8*2*nres)
4908 integer :: iwork(2*nres)
4909 !el common /przechowalnia/ ginvfric,Ghalf,fcopy
4910 integer :: ii,iti,k,l,nzero,nres2,nres6,ierr
4912 if (fg_rank.ne.king) goto 10
4917 if(.not.allocated(gamvec)) allocate(gamvec(nres2)) !(MAXRES2)
4918 if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
4919 !el if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
4920 !el allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
4921 if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !maxres2=2*maxres
4923 !el if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
4924 ! Zeroing out fricmat
4930 ! Load the friction coefficients corresponding to peptide groups
4936 ! Load the friction coefficients corresponding to side chains
4944 gamvec(ii)=gamsc(iabs(iti))
4946 if (surfarea) call sdarea(gamvec)
4948 ! write (iout,*) "Matrix A and vector gamma"
4950 ! write (iout,'(i2,$)') i
4952 ! write (iout,'(f4.1,$)') A(i,j)
4954 ! write (iout,'(f8.3)') gamvec(i)
4958 write (iout,*) "Vector gamvec"
4960 write (iout,'(i5,f10.5)') i, gamvec(i)
4964 ! The friction matrix
4969 dtdi=dtdi+A(j,k)*A(j,i)*gamvec(j)
4976 write (iout,'(//a)') "Matrix fricmat"
4977 call matout2(dimen,dimen,nres2,nres2,fricmat)
4979 if (lang.eq.2 .or. lang.eq.3) then
4980 ! Mass-scale the friction matrix if non-direct integration will be performed
4986 Ginvfric(i,j)=Ginvfric(i,j)+ &
4987 Gsqrm(i,k)*Gsqrm(l,j)*fricmat(k,l)
4992 ! Diagonalize the friction matrix
4997 Ghalf(ind)=Ginvfric(i,j)
5000 call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
5003 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
5004 " mass-scaled friction matrix"
5005 call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
5007 ! Precompute matrices for tinker stochastic integrator
5014 mt1(i,j)=mt1(i,j)+fricvec(k,i)*gsqrm(k,j)
5015 mt2(i,j)=mt2(i,j)+fricvec(k,i)*gsqrp(k,j)
5021 else if (lang.eq.4) then
5022 ! Diagonalize the friction matrix
5027 Ghalf(ind)=fricmat(i,j)
5030 call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
5033 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
5035 call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
5037 ! Determine the number of zero eigenvalues of the friction matrix
5038 nzero=max0(dimen-dimen1,0)
5039 ! do while (fricgam(nzero+1).le.1.0d-5 .and. nzero.lt.dimen)
5042 write (iout,*) "Number of zero eigenvalues:",nzero
5047 fricmat(i,j)=fricmat(i,j) &
5048 +fricvec(i,k)*fricvec(j,k)/fricgam(k)
5053 write (iout,'(//a)') "Generalized inverse of fricmat"
5054 call matout(dimen,dimen,nres6,nres6,fricmat)
5059 if (nfgtasks.gt.1) then
5060 if (fg_rank.eq.0) then
5061 ! The matching BROADCAST for fg processors is called in ERGASTULUM
5067 call MPI_Bcast(10,1,MPI_INTEGER,king,FG_COMM,IERROR)
5069 time_Bcast=time_Bcast+MPI_Wtime()-time00
5071 time_Bcast=time_Bcast+tcpu()-time00
5073 ! print *,"Processor",myrank,
5074 ! & " BROADCAST iorder in SETUP_FRICMAT"
5077 write (iout,*) "setup_fricmat licznik"!,licznik !sp
5083 ! Scatter the friction matrix
5084 call MPI_Scatterv(fricmat(1,1),nginv_counts(0),&
5085 nginv_start(0),MPI_DOUBLE_PRECISION,fcopy(1,1),&
5086 myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
5089 time_scatter=time_scatter+MPI_Wtime()-time00
5090 time_scatter_fmat=time_scatter_fmat+MPI_Wtime()-time00
5092 time_scatter=time_scatter+tcpu()-time00
5093 time_scatter_fmat=time_scatter_fmat+tcpu()-time00
5097 do j=1,2*my_ng_count
5098 fricmat(j,i)=fcopy(i,j)
5101 ! write (iout,*) "My chunk of fricmat"
5102 ! call MATOUT2(my_ng_count,dimen,maxres2,maxres2,fcopy)
5106 end subroutine setup_fricmat
5107 !-----------------------------------------------------------------------------
5108 subroutine stochastic_force(stochforcvec)
5111 use random, only:anorm_distr
5112 ! implicit real*8 (a-h,o-z)
5113 ! include 'DIMENSIONS'
5114 use control, only: tcpu
5119 ! include 'COMMON.VAR'
5120 ! include 'COMMON.CHAIN'
5121 ! include 'COMMON.DERIV'
5122 ! include 'COMMON.GEO'
5123 ! include 'COMMON.LOCAL'
5124 ! include 'COMMON.INTERACT'
5125 ! include 'COMMON.MD'
5126 ! include 'COMMON.TIME1'
5128 ! include 'COMMON.LANGEVIN'
5130 ! include 'COMMON.LANGEVIN.lang0'
5132 ! include 'COMMON.IOUNITS'
5134 real(kind=8) :: x,sig,lowb,highb
5135 real(kind=8) :: ff(3),force(3,0:2*nres),zeta2,lowb2
5136 real(kind=8) :: highb2,sig2,forcvec(6*nres),stochforcvec(6*nres)
5137 real(kind=8) :: time00
5138 logical :: lprn = .false.
5143 stochforc(j,i)=0.0d0
5153 ! Compute the stochastic forces acting on bodies. Store in force.
5159 force(j,i)=anorm_distr(x,sig,lowb,highb)
5167 force(j,i+nres)=anorm_distr(x,sig2,lowb2,highb2)
5171 time_fsample=time_fsample+MPI_Wtime()-time00
5173 time_fsample=time_fsample+tcpu()-time00
5175 ! Compute the stochastic forces acting on virtual-bond vectors.
5181 stochforc(j,i)=ff(j)+0.5d0*force(j,i)
5184 ff(j)=ff(j)+force(j,i)
5186 if (itype(i+1).ne.ntyp1) then
5188 stochforc(j,i)=stochforc(j,i)+force(j,i+nres+1)
5189 ff(j)=ff(j)+force(j,i+nres+1)
5194 stochforc(j,0)=ff(j)+force(j,nnt+nres)
5197 if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
5199 stochforc(j,i+nres)=force(j,i+nres)
5205 stochforcvec(j)=stochforc(j,0)
5210 stochforcvec(ind+j)=stochforc(j,i)
5215 if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
5217 stochforcvec(ind+j)=stochforc(j,i+nres)
5223 write (iout,*) "stochforcvec"
5225 write(iout,'(i5,e15.5)') i,stochforcvec(i)
5227 write(iout,*) "Stochastic forces backbone"
5229 write(iout,'(i5,3e15.5)') i,(stochforc(j,i),j=1,3)
5231 write(iout,*) "Stochastic forces side chain"
5233 write(iout,'(i5,3e15.5)') &
5234 i,(stochforc(j,i+nres),j=1,3)
5242 write (iout,*) i,ind
5244 forcvec(ind+j)=force(j,i)
5249 write (iout,*) i,ind
5251 forcvec(j+ind)=force(j,i+nres)
5256 write (iout,*) "forcvec"
5260 write (iout,'(2i3,2f10.5)') i,j,force(j,i),&
5267 write (iout,'(2i3,2f10.5)') i,j,force(j,i+nres),&
5276 end subroutine stochastic_force
5277 !-----------------------------------------------------------------------------
5278 subroutine sdarea(gamvec)
5280 ! Scale the friction coefficients according to solvent accessible surface areas
5281 ! Code adapted from TINKER
5285 ! implicit real*8 (a-h,o-z)
5286 ! include 'DIMENSIONS'
5287 ! include 'COMMON.CONTROL'
5288 ! include 'COMMON.VAR'
5289 ! include 'COMMON.MD'
5291 ! include 'COMMON.LANGEVIN'
5293 ! include 'COMMON.LANGEVIN.lang0'
5295 ! include 'COMMON.CHAIN'
5296 ! include 'COMMON.DERIV'
5297 ! include 'COMMON.GEO'
5298 ! include 'COMMON.LOCAL'
5299 ! include 'COMMON.INTERACT'
5300 ! include 'COMMON.IOUNITS'
5301 ! include 'COMMON.NAMES'
5302 real(kind=8),dimension(2*nres) :: radius,gamvec !(maxres2)
5303 real(kind=8),parameter :: twosix = 1.122462048309372981d0
5304 logical :: lprn = .false.
5305 real(kind=8) :: probe,area,ratio
5306 integer :: i,j,ind,iti
5308 ! determine new friction coefficients every few SD steps
5310 ! set the atomic radii to estimates of sigma values
5312 ! print *,"Entered sdarea"
5318 ! Load peptide group radii
5322 ! Load side chain radii
5325 radius(i+nres)=restok(iti)
5328 ! write (iout,*) "i",i," radius",radius(i)
5331 radius(i) = radius(i) / twosix
5332 if (radius(i) .ne. 0.0d0) radius(i) = radius(i) + probe
5335 ! scale atomic friction coefficients by accessible area
5337 if (lprn) write (iout,*) &
5338 "Original gammas, surface areas, scaling factors, new gammas, ",&
5339 "std's of stochastic forces"
5342 if (radius(i).gt.0.0d0) then
5343 call surfatom (i,area,radius)
5344 ratio = dmax1(area/(4.0d0*pi*radius(i)**2),1.0d-1)
5345 if (lprn) write (iout,'(i5,3f10.5,$)') &
5346 i,gamvec(ind+1),area,ratio
5349 gamvec(ind) = ratio * gamvec(ind)
5351 stdforcp(i)=stdfp*dsqrt(gamvec(ind))
5352 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcp(i)
5356 if (radius(i+nres).gt.0.0d0) then
5357 call surfatom (i+nres,area,radius)
5358 ratio = dmax1(area/(4.0d0*pi*radius(i+nres)**2),1.0d-1)
5359 if (lprn) write (iout,'(i5,3f10.5,$)') &
5360 i,gamvec(ind+1),area,ratio
5363 gamvec(ind) = ratio * gamvec(ind)
5365 stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamvec(ind))
5366 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcsc(i)
5371 end subroutine sdarea
5372 !-----------------------------------------------------------------------------
5374 !-----------------------------------------------------------------------------
5377 ! ###################################################
5378 ! ## COPYRIGHT (C) 1996 by Jay William Ponder ##
5379 ! ## All Rights Reserved ##
5380 ! ###################################################
5382 ! ################################################################
5384 ! ## subroutine surfatom -- exposed surface area of an atom ##
5386 ! ################################################################
5389 ! "surfatom" performs an analytical computation of the surface
5390 ! area of a specified atom; a simplified version of "surface"
5392 ! literature references:
5394 ! T. J. Richmond, "Solvent Accessible Surface Area and
5395 ! Excluded Volume in Proteins", Journal of Molecular Biology,
5398 ! L. Wesson and D. Eisenberg, "Atomic Solvation Parameters
5399 ! Applied to Molecular Dynamics of Proteins in Solution",
5400 ! Protein Science, 1, 227-235 (1992)
5402 ! variables and parameters:
5404 ! ir number of atom for which area is desired
5405 ! area accessible surface area of the atom
5406 ! radius radii of each of the individual atoms
5409 subroutine surfatom(ir,area,radius)
5411 ! implicit real*8 (a-h,o-z)
5412 ! include 'DIMENSIONS'
5414 ! include 'COMMON.GEO'
5415 ! include 'COMMON.IOUNITS'
5417 integer :: nsup,nstart_sup
5418 ! double precision c,dc,dc_old,d_c_work,xloc,xrot,dc_norm
5419 ! common /chain/ c(3,maxres2+2),dc(3,0:maxres2),dc_old(3,0:maxres2),
5420 ! & xloc(3,maxres),xrot(3,maxres),dc_norm(3,0:maxres2),
5421 ! & dc_work(MAXRES6),nres,nres0
5422 integer,parameter :: maxarc=300
5426 integer :: mi,ni,narc
5427 integer :: key(maxarc)
5428 integer :: intag(maxarc)
5429 integer :: intag1(maxarc)
5430 real(kind=8) :: area,arcsum
5431 real(kind=8) :: arclen,exang
5432 real(kind=8) :: delta,delta2
5433 real(kind=8) :: eps,rmove
5434 real(kind=8) :: xr,yr,zr
5435 real(kind=8) :: rr,rrsq
5436 real(kind=8) :: rplus,rminus
5437 real(kind=8) :: axx,axy,axz
5438 real(kind=8) :: ayx,ayy
5439 real(kind=8) :: azx,azy,azz
5440 real(kind=8) :: uxj,uyj,uzj
5441 real(kind=8) :: tx,ty,tz
5442 real(kind=8) :: txb,tyb,td
5443 real(kind=8) :: tr2,tr,txr,tyr
5444 real(kind=8) :: tk1,tk2
5445 real(kind=8) :: thec,the,t,tb
5446 real(kind=8) :: txk,tyk,tzk
5447 real(kind=8) :: t1,ti,tf,tt
5448 real(kind=8) :: txj,tyj,tzj
5449 real(kind=8) :: ccsq,cc,xysq
5450 real(kind=8) :: bsqk,bk,cosine
5451 real(kind=8) :: dsqj,gi,pix2
5452 real(kind=8) :: therk,dk,gk
5453 real(kind=8) :: risqk,rik
5454 real(kind=8) :: radius(2*nres) !(maxatm) (maxatm=maxres2)
5455 real(kind=8) :: ri(maxarc),risq(maxarc)
5456 real(kind=8) :: ux(maxarc),uy(maxarc),uz(maxarc)
5457 real(kind=8) :: xc(maxarc),yc(maxarc),zc(maxarc)
5458 real(kind=8) :: xc1(maxarc),yc1(maxarc),zc1(maxarc)
5459 real(kind=8) :: dsq(maxarc),bsq(maxarc)
5460 real(kind=8) :: dsq1(maxarc),bsq1(maxarc)
5461 real(kind=8) :: arci(maxarc),arcf(maxarc)
5462 real(kind=8) :: ex(maxarc),lt(maxarc),gr(maxarc)
5463 real(kind=8) :: b(maxarc),b1(maxarc),bg(maxarc)
5464 real(kind=8) :: kent(maxarc),kout(maxarc)
5465 real(kind=8) :: ther(maxarc)
5466 logical :: moved,top
5467 logical :: omit(maxarc)
5470 maxatm = 2*nres !maxres2 maxres2=2*maxres
5476 ! zero out the surface area for the sphere of interest
5479 ! write (2,*) "ir",ir," radius",radius(ir)
5480 if (radius(ir) .eq. 0.0d0) return
5482 ! set the overlap significance and connectivity shift
5486 delta2 = delta * delta
5491 ! store coordinates and radius of the sphere of interest
5499 ! initialize values of some counters and summations
5508 ! test each sphere to see if it overlaps the sphere of interest
5511 if (i.eq.ir .or. radius(i).eq.0.0d0) goto 30
5512 rplus = rr + radius(i)
5514 if (abs(tx) .ge. rplus) goto 30
5516 if (abs(ty) .ge. rplus) goto 30
5518 if (abs(tz) .ge. rplus) goto 30
5520 ! check for sphere overlap by testing distance against radii
5522 xysq = tx*tx + ty*ty
5523 if (xysq .lt. delta2) then
5530 if (rplus-cc .le. delta) goto 30
5531 rminus = rr - radius(i)
5533 ! check to see if sphere of interest is completely buried
5535 if (cc-abs(rminus) .le. delta) then
5536 if (rminus .le. 0.0d0) goto 170
5540 ! check for too many overlaps with sphere of interest
5542 if (io .ge. maxarc) then
5544 20 format (/,' SURFATOM -- Increase the Value of MAXARC')
5548 ! get overlap between current sphere and sphere of interest
5557 gr(io) = (ccsq+rplus*rminus) / (2.0d0*rr*b1(io))
5563 ! case where no other spheres overlap the sphere of interest
5566 area = 4.0d0 * pi * rrsq
5570 ! case where only one sphere overlaps the sphere of interest
5573 area = pix2 * (1.0d0 + gr(1))
5574 area = mod(area,4.0d0*pi) * rrsq
5578 ! case where many spheres intersect the sphere of interest;
5579 ! sort the intersecting spheres by their degree of overlap
5581 call sort2 (io,gr,key)
5584 intag(i) = intag1(k)
5593 ! get radius of each overlap circle on surface of the sphere
5598 risq(i) = rrsq - gi*gi
5599 ri(i) = sqrt(risq(i))
5600 ther(i) = 0.5d0*pi - asin(min(1.0d0,max(-1.0d0,gr(i))))
5603 ! find boundary of inaccessible area on sphere of interest
5606 if (.not. omit(k)) then
5613 ! check to see if J circle is intersecting K circle;
5614 ! get distance between circle centers and sum of radii
5617 if (omit(j)) goto 60
5618 cc = (txk*xc(j)+tyk*yc(j)+tzk*zc(j))/(bk*b(j))
5619 cc = acos(min(1.0d0,max(-1.0d0,cc)))
5620 td = therk + ther(j)
5622 ! check to see if circles enclose separate regions
5624 if (cc .ge. td) goto 60
5626 ! check for circle J completely inside circle K
5628 if (cc+ther(j) .lt. therk) goto 40
5630 ! check for circles that are essentially parallel
5632 if (cc .gt. delta) goto 50
5637 ! check to see if sphere of interest is completely buried
5640 if (pix2-cc .le. td) goto 170
5646 ! find T value of circle intersections
5649 if (omit(k)) goto 110
5664 ! rotation matrix elements
5676 if (.not. omit(j)) then
5681 ! rotate spheres so K vector colinear with z-axis
5683 uxj = txj*axx + tyj*axy - tzj*axz
5684 uyj = tyj*ayy - txj*ayx
5685 uzj = txj*azx + tyj*azy + tzj*azz
5686 cosine = min(1.0d0,max(-1.0d0,uzj/b(j)))
5687 if (acos(cosine) .lt. therk+ther(j)) then
5688 dsqj = uxj*uxj + uyj*uyj
5693 tr2 = risqk*dsqj - tb*tb
5699 ! get T values of intersection for K circle
5702 tb = min(1.0d0,max(-1.0d0,tb))
5704 if (tyb-txr .lt. 0.0d0) tk1 = pix2 - tk1
5706 tb = min(1.0d0,max(-1.0d0,tb))
5708 if (tyb+txr .lt. 0.0d0) tk2 = pix2 - tk2
5709 thec = (rrsq*uzj-gk*bg(j)) / (rik*ri(j)*b(j))
5710 if (abs(thec) .lt. 1.0d0) then
5712 else if (thec .ge. 1.0d0) then
5714 else if (thec .le. -1.0d0) then
5718 ! see if "tk1" is entry or exit point; check t=0 point;
5719 ! "ti" is exit point, "tf" is entry point
5721 cosine = min(1.0d0,max(-1.0d0, &
5722 (uzj*gk-uxj*rik)/(b(j)*rr)))
5723 if ((acos(cosine)-ther(j))*(tk2-tk1) .le. 0.0d0) then
5731 if (narc .ge. maxarc) then
5733 70 format (/,' SURFATOM -- Increase the Value',&
5737 if (tf .le. ti) then
5758 ! special case; K circle without intersections
5760 if (narc .le. 0) goto 90
5762 ! general case; sum up arclength and set connectivity code
5764 call sort2 (narc,arci,key)
5769 if (narc .gt. 1) then
5772 if (t .lt. arci(j)) then
5773 arcsum = arcsum + arci(j) - t
5774 exang = exang + ex(ni)
5776 if (jb .ge. maxarc) then
5778 80 format (/,' SURFATOM -- Increase the Value',&
5783 kent(jb) = maxarc*i + k
5785 kout(jb) = maxarc*k + i
5794 arcsum = arcsum + pix2 - t
5796 exang = exang + ex(ni)
5799 kent(jb) = maxarc*i + k
5801 kout(jb) = maxarc*k + i
5808 arclen = arclen + gr(k)*arcsum
5811 if (arclen .eq. 0.0d0) goto 170
5812 if (jb .eq. 0) goto 150
5814 ! find number of independent boundaries and check connectivity
5818 if (kout(k) .ne. 0) then
5825 if (m .eq. kent(ii)) then
5828 if (j .eq. jb) goto 150
5840 ! attempt to fix connectivity error by moving atom slightly
5844 140 format (/,' SURFATOM -- Connectivity Error at Atom',i6)
5853 ! compute the exposed surface area for the sphere of interest
5856 area = ib*pix2 + exang + arclen
5857 area = mod(area,4.0d0*pi) * rrsq
5859 ! attempt to fix negative area by moving atom slightly
5861 if (area .lt. 0.0d0) then
5864 160 format (/,' SURFATOM -- Negative Area at Atom',i6)
5875 end subroutine surfatom
5876 !----------------------------------------------------------------
5877 !----------------------------------------------------------------
5878 subroutine alloc_MD_arrays
5879 !EL Allocation of arrays used by MD module
5881 integer :: nres2,nres6
5884 !----------------------
5888 allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
5889 allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
5890 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
5891 allocate(fricvec(nres2,nres2))
5892 allocate(pfric_mat(nres2,nres2),vfric_mat(nres2,nres2))
5893 allocate(afric_mat(nres2,nres2),prand_mat(nres2,nres2))
5894 allocate(vrand_mat1(nres2,nres2),vrand_mat2(nres2,nres2)) !(MAXRES2,MAXRES2)
5895 allocate(pfric0_mat(nres2,nres2,0:maxflag_stoch))
5896 allocate(afric0_mat(nres2,nres2,0:maxflag_stoch))
5897 allocate(vfric0_mat(nres2,nres2,0:maxflag_stoch))
5898 allocate(prand0_mat(nres2,nres2,0:maxflag_stoch))
5899 allocate(vrand0_mat1(nres2,nres2,0:maxflag_stoch))
5900 allocate(vrand0_mat2(nres2,nres2,0:maxflag_stoch)) !(MAXRES2,MAXRES2,0:maxflag_stoch)
5901 allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
5903 allocate(mt1(nres2,nres2),mt2(nres2,nres2),mt3(nres2,nres2)) !(maxres2,maxres2)
5904 !----------------------
5906 ! commom.langevin.lang0
5908 allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
5909 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
5910 allocate(fricvec(nres2,nres2)) !(MAXRES2,MAXRES2)
5911 allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
5912 allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
5915 !el if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2))
5916 !----------------------
5917 ! commom.hairpin in CSA module
5918 !----------------------
5919 ! common.mce in MCM_MD module
5920 !----------------------
5922 ! common /mdgrad/ in module.energy
5923 ! common /back_constr/ in module.energy
5924 ! common /qmeas/ in module.energy
5927 allocate(potEcomp(0:n_ene+4)) !(0:n_ene+4)
5929 allocate(d_t(3,0:nres2),d_a(3,0:nres2),d_t_old(3,0:nres2)) !(3,0:MAXRES2)
5930 allocate(d_a_work(nres6)) !(6*MAXRES)
5932 allocate(DM(nres2),DU1(nres2),DU2(nres2))
5933 allocate(DMorig(nres2),DU1orig(nres2),DU2orig(nres2))
5935 allocate(Gmat(nres2,nres2),A(nres2,nres2))
5936 if(.not.allocated(Ginv)) allocate(Ginv(nres2,nres2)) !in control: ergastulum
5937 allocate(Gsqrp(nres2,nres2),Gsqrm(nres2,nres2),Gvec(nres2,nres2)) !(maxres2,maxres2)
5939 allocate(Geigen(nres2)) !(maxres2)
5940 if(.not.allocated(vtot)) allocate(vtot(nres2)) !(maxres2)
5941 ! common /inertia/ in io_conf: parmread
5942 ! real(kind=8),dimension(:),allocatable :: ISC,msc !(ntyp+1)
5943 ! common /langevin/in io read_MDpar
5944 ! real(kind=8),dimension(:),allocatable :: gamsc !(ntyp1)
5945 ! real(kind=8),dimension(:),allocatable :: stdfsc !(ntyp)
5946 ! in io_conf: parmread
5947 ! real(kind=8),dimension(:),allocatable :: restok !(ntyp+1)
5948 ! common /mdpmpi/ in control: ergastulum
5949 if(.not.allocated(ng_start)) allocate(ng_start(0:nfgtasks-1))
5950 if(.not.allocated(ng_counts)) allocate(ng_counts(0:nfgtasks-1))
5951 if(.not.allocated(nginv_counts)) allocate(nginv_counts(0:nfgtasks-1)) !(0:MaxProcs-1)
5952 if(.not.allocated(nginv_start)) allocate(nginv_start(0:nfgtasks)) !(0:MaxProcs)
5953 !----------------------
5954 ! common.muca in read_muca
5955 ! common /double_muca/
5956 ! real(kind=8) :: elow,ehigh,factor,hbin,factor_min
5957 ! real(kind=8),dimension(:),allocatable :: emuca,nemuca,&
5958 ! nemuca2,hist !(4*maxres)
5959 ! common /integer_muca/
5960 ! integer :: nmuca,imtime,muca_smooth
5962 ! real(kind=8),dimension(:),allocatable :: elowi,ehighi !(maxprocs)
5963 !----------------------
5965 ! common /mdgrad/ in module.energy
5966 ! common /back_constr/ in module.energy
5967 ! common /qmeas/ in module.energy
5971 allocate(d_t_work(nres6),d_t_work_new(nres6),d_af_work(nres6))
5972 allocate(d_as_work(nres6),kinetic_force(nres6)) !(MAXRES6)
5973 allocate(d_t_new(3,0:nres2),d_a_old(3,0:nres2),d_a_short(3,0:nres2)) !,d_a !(3,0:MAXRES2)
5974 allocate(stdforcp(nres),stdforcsc(nres)) !(MAXRES)
5975 !----------------------
5977 allocate(D_ban(nres6)) !(MAXRES6) maxres6=6*maxres
5978 ! common /stochcalc/ stochforcvec
5979 allocate(stochforcvec(nres6)) !(MAXRES6) maxres6=6*maxres
5980 !----------------------
5982 end subroutine alloc_MD_arrays
5983 !-----------------------------------------------------------------------------
5984 !-----------------------------------------------------------------------------