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 ! AL 4/17/17: Reduce the steps if NaNs occurred.
1167 if (potEcomp(0).gt.0.99e20 .or. isnan(potEcomp(0)).gt.0) then
1172 if (large.and. mod(itime,ntwe).eq.0) &
1173 call enerprint(potEcomp)
1176 t_etotal=t_etotal+MPI_Wtime()-tt0
1178 t_etotal=t_etotal+tcpu()-tt0
1181 potE=potEcomp(0)-potEcomp(20)
1183 ! Get the new accelerations
1186 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1188 t_enegrad=t_enegrad+tcpu()-tt0
1190 ! Determine maximum acceleration and scale down the timestep if needed
1192 amax=amax/(itime_scal+1)**2
1193 call predict_edrift(epdrift)
1194 if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then
1195 ! Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
1197 ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax)) &
1199 itime_scal=itime_scal+ifac_time
1200 ! fac_time=dmin1(damax/amax,0.5d0)
1201 fac_time=0.5d0**ifac_time
1202 d_time=d_time*fac_time
1203 if (lang.eq.2 .or. lang.eq.3) then
1205 ! write (iout,*) "Calling sd_verlet_setup: 1"
1206 ! Rescale the stochastic forces and recalculate or restore
1207 ! the matrices of tinker integrator
1208 if (itime_scal.gt.maxflag_stoch) then
1209 if (large) write (iout,'(a,i5,a)') &
1210 "Calculate matrices for stochastic step;",&
1211 " itime_scal ",itime_scal
1213 call sd_verlet_p_setup
1215 call sd_verlet_ciccotti_setup
1217 write (iout,'(2a,i3,a,i3,1h.)') &
1218 "Warning: cannot store matrices for stochastic",&
1219 " integration because the index",itime_scal,&
1220 " is greater than",maxflag_stoch
1221 write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",&
1222 " integration Langevin algorithm for better efficiency."
1223 else if (flag_stoch(itime_scal)) then
1224 if (large) write (iout,'(a,i5,a,l1)') &
1225 "Restore matrices for stochastic step; itime_scal ",&
1226 itime_scal," flag ",flag_stoch(itime_scal)
1229 pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
1230 afric_mat(i,j)=afric0_mat(i,j,itime_scal)
1231 vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
1232 prand_mat(i,j)=prand0_mat(i,j,itime_scal)
1233 vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
1234 vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
1238 if (large) write (iout,'(2a,i5,a,l1)') &
1239 "Calculate & store matrices for stochastic step;",&
1240 " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
1242 call sd_verlet_p_setup
1244 call sd_verlet_ciccotti_setup
1246 flag_stoch(ifac_time)=.true.
1249 pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
1250 afric0_mat(i,j,itime_scal)=afric_mat(i,j)
1251 vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
1252 prand0_mat(i,j,itime_scal)=prand_mat(i,j)
1253 vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
1254 vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
1258 fac_time=1.0d0/dsqrt(fac_time)
1260 stochforcvec(i)=fac_time*stochforcvec(i)
1263 else if (lang.eq.1) then
1264 ! Rescale the accelerations due to stochastic forces
1265 fac_time=1.0d0/dsqrt(fac_time)
1267 d_as_work(i)=d_as_work(i)*fac_time
1270 if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)') &
1271 "itime",itime," Timestep scaled down to ",&
1272 d_time," ifac_time",ifac_time," itime_scal",itime_scal
1274 ! Second step of the velocity Verlet algorithm
1279 else if (lang.eq.3) then
1281 call sd_verlet2_ciccotti
1283 else if (lang.eq.1) then
1288 if (rattle) call rattle2
1290 if (d_time.ne.d_time0) then
1293 if (lang.eq.2 .or. lang.eq.3) then
1294 if (large) write (iout,'(a)') &
1295 "Restore original matrices for stochastic step"
1296 ! write (iout,*) "Calling sd_verlet_setup: 2"
1297 ! Restore the matrices of tinker integrator if the time step has been restored
1300 pfric_mat(i,j)=pfric0_mat(i,j,0)
1301 afric_mat(i,j)=afric0_mat(i,j,0)
1302 vfric_mat(i,j)=vfric0_mat(i,j,0)
1303 prand_mat(i,j)=prand0_mat(i,j,0)
1304 vrand_mat1(i,j)=vrand0_mat1(i,j,0)
1305 vrand_mat2(i,j)=vrand0_mat2(i,j,0)
1314 ! Calculate the kinetic and the total energy and the kinetic temperature
1318 ! call kinetic1(EK1)
1319 ! write (iout,*) "step",itime," EK",EK," EK1",EK1
1321 ! Couple the system to Berendsen bath if needed
1322 if (tbf .and. lang.eq.0) then
1325 kinetic_T=2.0d0/(dimen3*Rb)*EK
1326 ! Backup the coordinates, velocities, and accelerations
1330 d_t_old(j,i)=d_t(j,i)
1331 d_a_old(j,i)=d_a(j,i)
1335 if (mod(itime,ntwe).eq.0 .and. large) then
1336 write (iout,*) "Velocities, step 2"
1338 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1339 (d_t(j,i+nres),j=1,3)
1344 end subroutine velverlet_step
1345 !-----------------------------------------------------------------------------
1346 subroutine RESPA_step(itime)
1347 !-------------------------------------------------------------------------------
1348 ! Perform a single RESPA step.
1349 !-------------------------------------------------------------------------------
1350 ! implicit real*8 (a-h,o-z)
1351 ! include 'DIMENSIONS'
1355 use control, only:tcpu
1357 ! use io_conf, only:cartprint
1360 integer :: IERROR,ERRCODE
1362 ! include 'COMMON.SETUP'
1363 ! include 'COMMON.CONTROL'
1364 ! include 'COMMON.VAR'
1365 ! include 'COMMON.MD'
1367 ! include 'COMMON.LANGEVIN'
1369 ! include 'COMMON.LANGEVIN.lang0'
1371 ! include 'COMMON.CHAIN'
1372 ! include 'COMMON.DERIV'
1373 ! include 'COMMON.GEO'
1374 ! include 'COMMON.LOCAL'
1375 ! include 'COMMON.INTERACT'
1376 ! include 'COMMON.IOUNITS'
1377 ! include 'COMMON.NAMES'
1378 ! include 'COMMON.TIME1'
1379 real(kind=8),dimension(0:n_ene) :: energia_short,energia_long
1380 real(kind=8),dimension(3) :: L,vcm,incr
1381 real(kind=8),dimension(3,0:2*nres) :: dc_old0,d_t_old0,d_a_old0 !(3,0:maxres2) maxres2=2*maxres
1382 logical :: PRINT_AMTS_MSG = .false.
1383 integer :: count,rstcount !ilen,
1385 character(len=50) :: tytul
1386 integer :: maxcount_scale = 10
1387 !el common /gucio/ cm
1388 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1389 !el common /stochcalc/ stochforcvec
1390 integer :: itime,itt,i,j,itsplit
1392 !el common /cipiszcze/ itt
1394 real(kind=8) :: epdrift,tt0,epdriftmax
1397 if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres)) !(MAXRES6) maxres6=6*maxres
1401 if (large.and. mod(itime,ntwe).eq.0) then
1402 write (iout,*) "***************** RESPA itime",itime
1403 write (iout,*) "Cartesian and internal coordinates: step 0"
1405 call pdbout(0.0d0,"cipiszcze",iout)
1407 write (iout,*) "Accelerations from long-range forces"
1409 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1410 (d_a(j,i+nres),j=1,3)
1412 write (iout,*) "Velocities, step 0"
1414 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1415 (d_t(j,i+nres),j=1,3)
1420 ! Perform the initial RESPA step (increment velocities)
1421 ! write (iout,*) "*********************** RESPA ini"
1424 if (mod(itime,ntwe).eq.0 .and. large) then
1425 write (iout,*) "Velocities, end"
1427 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1428 (d_t(j,i+nres),j=1,3)
1432 ! Compute the short-range forces
1438 ! 7/2/2009 commented out
1440 ! call etotal_short(energia_short)
1443 ! 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
1446 d_a(j,i)=d_a_short(j,i)
1450 if (large.and. mod(itime,ntwe).eq.0) then
1451 write (iout,*) "energia_short",energia_short(0)
1452 write (iout,*) "Accelerations from short-range forces"
1454 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1455 (d_a(j,i+nres),j=1,3)
1460 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1462 t_enegrad=t_enegrad+tcpu()-tt0
1467 d_t_old(j,i)=d_t(j,i)
1468 d_a_old(j,i)=d_a(j,i)
1471 ! 6/30/08 A-MTS: attempt at increasing the split number
1474 dc_old0(j,i)=dc_old(j,i)
1475 d_t_old0(j,i)=d_t_old(j,i)
1476 d_a_old0(j,i)=d_a_old(j,i)
1479 if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
1480 if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
1487 ! write (iout,*) "itime",itime," ntime_split",ntime_split
1488 ! Split the time step
1489 d_time=d_time0/ntime_split
1490 ! Perform the short-range RESPA steps (velocity Verlet increments of
1491 ! positions and velocities using short-range forces)
1492 ! write (iout,*) "*********************** RESPA split"
1493 do itsplit=1,ntime_split
1496 else if (lang.eq.2 .or. lang.eq.3) then
1498 call stochastic_force(stochforcvec)
1501 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
1503 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1508 ! First step of the velocity Verlet algorithm
1513 else if (lang.eq.3) then
1515 call sd_verlet1_ciccotti
1517 else if (lang.eq.1) then
1522 ! Build the chain from the newly calculated coordinates
1523 call chainbuild_cart
1524 if (rattle) call rattle1
1526 if (large.and. mod(itime,ntwe).eq.0) then
1527 write (iout,*) "***** ITSPLIT",itsplit
1528 write (iout,*) "Cartesian and internal coordinates: step 1"
1529 call pdbout(0.0d0,"cipiszcze",iout)
1532 write (iout,*) "Velocities, step 1"
1534 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1535 (d_t(j,i+nres),j=1,3)
1544 ! Calculate energy and forces
1546 call etotal_short(energia_short)
1547 ! AL 4/17/17: Exit itime_split loop when energy goes infinite
1548 if (energia_short(0).gt.0.99e20 .or. isnan(energia_short(0)) ) then
1549 if (PRINT_AMTS_MSG) &
1550 write (iout,*) "Infinities/NaNs in energia_short",energia_short(0),"; increasing ntime_split to",ntime_split
1551 ntime_split=ntime_split*2
1552 if (ntime_split.gt.maxtime_split) then
1555 "Cannot rescue the run; aborting job. Retry with a smaller time step"
1557 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1560 "Cannot rescue the run; terminating. Retry with a smaller time step"
1566 if (large.and. mod(itime,ntwe).eq.0) &
1567 call enerprint(energia_short)
1570 t_eshort=t_eshort+MPI_Wtime()-tt0
1572 t_eshort=t_eshort+tcpu()-tt0
1576 ! Get the new accelerations
1578 ! 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
1581 d_a_short(j,i)=d_a(j,i)
1585 if (large.and. mod(itime,ntwe).eq.0) then
1586 write (iout,*)"energia_short",energia_short(0)
1587 write (iout,*) "Accelerations from short-range forces"
1589 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1590 (d_a(j,i+nres),j=1,3)
1595 ! Determine maximum acceleration and scale down the timestep if needed
1597 amax=amax/ntime_split**2
1598 call predict_edrift(epdrift)
1599 if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0) &
1600 write (iout,*) "amax",amax," damax",damax,&
1601 " epdrift",epdrift," epdriftmax",epdriftmax
1602 ! Exit loop and try with increased split number if the change of
1603 ! acceleration is too big
1604 if (amax.gt.damax .or. epdrift.gt.edriftmax) then
1605 if (ntime_split.lt.maxtime_split) then
1607 ntime_split=ntime_split*2
1608 ! AL 4/17/17: We should exit the itime_split loop when acceleration change is too big
1612 dc_old(j,i)=dc_old0(j,i)
1613 d_t_old(j,i)=d_t_old0(j,i)
1614 d_a_old(j,i)=d_a_old0(j,i)
1617 if (PRINT_AMTS_MSG) then
1618 write (iout,*) "acceleration/energy drift too large",amax,&
1619 epdrift," split increased to ",ntime_split," itime",itime,&
1625 "Uh-hu. Bumpy landscape. Maximum splitting number",&
1627 " already reached!!! Trying to carry on!"
1631 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1633 t_enegrad=t_enegrad+tcpu()-tt0
1635 ! Second step of the velocity Verlet algorithm
1640 else if (lang.eq.3) then
1642 call sd_verlet2_ciccotti
1644 else if (lang.eq.1) then
1649 if (rattle) call rattle2
1650 ! Backup the coordinates, velocities, and accelerations
1654 d_t_old(j,i)=d_t(j,i)
1655 d_a_old(j,i)=d_a(j,i)
1662 ! Restore the time step
1664 ! Compute long-range forces
1671 call etotal_long(energia_long)
1672 if (energia_long(0).gt.0.99e20 .or. isnan(energia_long(0))) then
1675 "Infinitied/NaNs in energia_long, Aborting MPI job."
1677 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1679 write (iout,*) "Infinitied/NaNs in energia_long, terminating."
1683 if (large.and. mod(itime,ntwe).eq.0) &
1684 call enerprint(energia_long)
1687 t_elong=t_elong+MPI_Wtime()-tt0
1689 t_elong=t_elong+tcpu()-tt0
1695 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1697 t_enegrad=t_enegrad+tcpu()-tt0
1699 ! Compute accelerations from long-range forces
1701 if (large.and. mod(itime,ntwe).eq.0) then
1702 write (iout,*) "energia_long",energia_long(0)
1703 write (iout,*) "Cartesian and internal coordinates: step 2"
1705 call pdbout(0.0d0,"cipiszcze",iout)
1707 write (iout,*) "Accelerations from long-range forces"
1709 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1710 (d_a(j,i+nres),j=1,3)
1712 write (iout,*) "Velocities, step 2"
1714 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1715 (d_t(j,i+nres),j=1,3)
1719 ! Compute the final RESPA step (increment velocities)
1720 ! write (iout,*) "*********************** RESPA fin"
1722 ! Compute the complete potential energy
1724 potEcomp(i)=energia_short(i)+energia_long(i)
1726 potE=potEcomp(0)-potEcomp(20)
1727 ! potE=energia_short(0)+energia_long(0)
1729 ! Calculate the kinetic and the total energy and the kinetic temperature
1732 ! Couple the system to Berendsen bath if needed
1733 if (tbf .and. lang.eq.0) then
1736 kinetic_T=2.0d0/(dimen3*Rb)*EK
1737 ! Backup the coordinates, velocities, and accelerations
1739 if (mod(itime,ntwe).eq.0 .and. large) then
1740 write (iout,*) "Velocities, end"
1742 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1743 (d_t(j,i+nres),j=1,3)
1748 end subroutine RESPA_step
1749 !-----------------------------------------------------------------------------
1750 subroutine RESPA_vel
1751 ! First and last RESPA step (incrementing velocities using long-range
1754 ! implicit real*8 (a-h,o-z)
1755 ! include 'DIMENSIONS'
1756 ! include 'COMMON.CONTROL'
1757 ! include 'COMMON.VAR'
1758 ! include 'COMMON.MD'
1759 ! include 'COMMON.CHAIN'
1760 ! include 'COMMON.DERIV'
1761 ! include 'COMMON.GEO'
1762 ! include 'COMMON.LOCAL'
1763 ! include 'COMMON.INTERACT'
1764 ! include 'COMMON.IOUNITS'
1765 ! include 'COMMON.NAMES'
1766 integer :: i,j,inres
1769 d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
1773 d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
1777 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1780 d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
1785 end subroutine RESPA_vel
1786 !-----------------------------------------------------------------------------
1788 ! Applying velocity Verlet algorithm - step 1 to coordinates
1790 ! implicit real*8 (a-h,o-z)
1791 ! include 'DIMENSIONS'
1792 ! include 'COMMON.CONTROL'
1793 ! include 'COMMON.VAR'
1794 ! include 'COMMON.MD'
1795 ! include 'COMMON.CHAIN'
1796 ! include 'COMMON.DERIV'
1797 ! include 'COMMON.GEO'
1798 ! include 'COMMON.LOCAL'
1799 ! include 'COMMON.INTERACT'
1800 ! include 'COMMON.IOUNITS'
1801 ! include 'COMMON.NAMES'
1802 real(kind=8) :: adt,adt2
1803 integer :: i,j,inres
1806 write (iout,*) "VELVERLET1 START: DC"
1808 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),&
1809 (dc(j,i+nres),j=1,3)
1813 adt=d_a_old(j,0)*d_time
1815 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1816 d_t_new(j,0)=d_t_old(j,0)+adt2
1817 d_t(j,0)=d_t_old(j,0)+adt
1821 adt=d_a_old(j,i)*d_time
1823 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1824 d_t_new(j,i)=d_t_old(j,i)+adt2
1825 d_t(j,i)=d_t_old(j,i)+adt
1829 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1832 adt=d_a_old(j,inres)*d_time
1834 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1835 d_t_new(j,inres)=d_t_old(j,inres)+adt2
1836 d_t(j,inres)=d_t_old(j,inres)+adt
1841 write (iout,*) "VELVERLET1 END: DC"
1843 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),&
1844 (dc(j,i+nres),j=1,3)
1848 end subroutine verlet1
1849 !-----------------------------------------------------------------------------
1851 ! Step 2 of the velocity Verlet algorithm: update velocities
1853 ! implicit real*8 (a-h,o-z)
1854 ! include 'DIMENSIONS'
1855 ! include 'COMMON.CONTROL'
1856 ! include 'COMMON.VAR'
1857 ! include 'COMMON.MD'
1858 ! include 'COMMON.CHAIN'
1859 ! include 'COMMON.DERIV'
1860 ! include 'COMMON.GEO'
1861 ! include 'COMMON.LOCAL'
1862 ! include 'COMMON.INTERACT'
1863 ! include 'COMMON.IOUNITS'
1864 ! include 'COMMON.NAMES'
1865 integer :: i,j,inres
1868 d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
1872 d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
1876 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1879 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
1884 end subroutine verlet2
1885 !-----------------------------------------------------------------------------
1886 subroutine sddir_precalc
1887 ! Applying velocity Verlet algorithm - step 1 to coordinates
1888 ! implicit real*8 (a-h,o-z)
1889 ! include 'DIMENSIONS'
1895 ! include 'COMMON.CONTROL'
1896 ! include 'COMMON.VAR'
1897 ! include 'COMMON.MD'
1899 ! include 'COMMON.LANGEVIN'
1901 ! include 'COMMON.LANGEVIN.lang0'
1903 ! include 'COMMON.CHAIN'
1904 ! include 'COMMON.DERIV'
1905 ! include 'COMMON.GEO'
1906 ! include 'COMMON.LOCAL'
1907 ! include 'COMMON.INTERACT'
1908 ! include 'COMMON.IOUNITS'
1909 ! include 'COMMON.NAMES'
1910 ! include 'COMMON.TIME1'
1911 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1912 !el common /stochcalc/ stochforcvec
1913 real(kind=8) :: time00
1915 ! Compute friction and stochastic forces
1920 time_fric=time_fric+MPI_Wtime()-time00
1922 call stochastic_force(stochforcvec)
1923 time_stoch=time_stoch+MPI_Wtime()-time00
1926 ! Compute the acceleration due to friction forces (d_af_work) and stochastic
1927 ! forces (d_as_work)
1929 call ginv_mult(fric_work, d_af_work)
1930 call ginv_mult(stochforcvec, d_as_work)
1932 end subroutine sddir_precalc
1933 !-----------------------------------------------------------------------------
1934 subroutine sddir_verlet1
1935 ! Applying velocity Verlet algorithm - step 1 to velocities
1938 ! implicit real*8 (a-h,o-z)
1939 ! include 'DIMENSIONS'
1940 ! include 'COMMON.CONTROL'
1941 ! include 'COMMON.VAR'
1942 ! include 'COMMON.MD'
1944 ! include 'COMMON.LANGEVIN'
1946 ! include 'COMMON.LANGEVIN.lang0'
1948 ! include 'COMMON.CHAIN'
1949 ! include 'COMMON.DERIV'
1950 ! include 'COMMON.GEO'
1951 ! include 'COMMON.LOCAL'
1952 ! include 'COMMON.INTERACT'
1953 ! include 'COMMON.IOUNITS'
1954 ! include 'COMMON.NAMES'
1955 ! Revised 3/31/05 AL: correlation between random contributions to
1956 ! position and velocity increments included.
1957 real(kind=8) :: sqrt13 = 0.57735026918962576451d0 ! 1/sqrt(3)
1958 real(kind=8) :: adt,adt2
1959 integer :: i,j,ind,inres
1961 ! Add the contribution from BOTH friction and stochastic force to the
1962 ! coordinates, but ONLY the contribution from the friction forces to velocities
1965 adt=(d_a_old(j,0)+d_af_work(j))*d_time
1966 adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
1967 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1968 d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
1969 d_t(j,0)=d_t_old(j,0)+adt
1974 adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
1975 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1976 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1977 d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
1978 d_t(j,i)=d_t_old(j,i)+adt
1983 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1986 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
1987 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1988 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1989 d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
1990 d_t(j,inres)=d_t_old(j,inres)+adt
1996 end subroutine sddir_verlet1
1997 !-----------------------------------------------------------------------------
1998 subroutine sddir_verlet2
1999 ! Calculating the adjusted velocities for accelerations
2002 ! implicit real*8 (a-h,o-z)
2003 ! include 'DIMENSIONS'
2004 ! include 'COMMON.CONTROL'
2005 ! include 'COMMON.VAR'
2006 ! include 'COMMON.MD'
2008 ! include 'COMMON.LANGEVIN'
2010 ! include 'COMMON.LANGEVIN.lang0'
2012 ! include 'COMMON.CHAIN'
2013 ! include 'COMMON.DERIV'
2014 ! include 'COMMON.GEO'
2015 ! include 'COMMON.LOCAL'
2016 ! include 'COMMON.INTERACT'
2017 ! include 'COMMON.IOUNITS'
2018 ! include 'COMMON.NAMES'
2019 real(kind=8),dimension(6*nres) :: stochforcvec,d_as_work1 !(MAXRES6) maxres6=6*maxres
2020 real(kind=8) :: cos60 = 0.5d0, sin60 = 0.86602540378443864676d0
2021 integer :: i,j,ind,inres
2022 ! Revised 3/31/05 AL: correlation between random contributions to
2023 ! position and velocity increments included.
2024 ! The correlation coefficients are calculated at low-friction limit.
2025 ! Also, friction forces are now not calculated with new velocities.
2027 ! call friction_force
2028 call stochastic_force(stochforcvec)
2030 ! Compute the acceleration due to friction forces (d_af_work) and stochastic
2031 ! forces (d_as_work)
2033 call ginv_mult(stochforcvec, d_as_work1)
2039 d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j)) &
2040 +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
2045 d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j)) &
2046 +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
2051 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2054 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres) &
2055 +d_af_work(ind+j))+sin60*d_as_work(ind+j) &
2056 +cos60*d_as_work1(ind+j))*d_time
2062 end subroutine sddir_verlet2
2063 !-----------------------------------------------------------------------------
2064 subroutine max_accel
2066 ! Find the maximum difference in the accelerations of the the sites
2067 ! at the beginning and the end of the time step.
2070 ! implicit real*8 (a-h,o-z)
2071 ! include 'DIMENSIONS'
2072 ! include 'COMMON.CONTROL'
2073 ! include 'COMMON.VAR'
2074 ! include 'COMMON.MD'
2075 ! include 'COMMON.CHAIN'
2076 ! include 'COMMON.DERIV'
2077 ! include 'COMMON.GEO'
2078 ! include 'COMMON.LOCAL'
2079 ! include 'COMMON.INTERACT'
2080 ! include 'COMMON.IOUNITS'
2081 real(kind=8),dimension(3) :: aux,accel,accel_old
2082 real(kind=8) :: dacc
2086 ! aux(j)=d_a(j,0)-d_a_old(j,0)
2087 accel_old(j)=d_a_old(j,0)
2094 ! 7/3/08 changed to asymmetric difference
2096 ! accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
2097 accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
2098 accel(j)=accel(j)+0.5d0*d_a(j,i)
2099 ! if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
2100 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
2101 dacc=dabs(accel(j)-accel_old(j))
2102 ! write (iout,*) i,dacc
2103 if (dacc.gt.amax) amax=dacc
2111 accel_old(j)=d_a_old(j,0)
2116 accel_old(j)=accel_old(j)+d_a_old(j,1)
2117 accel(j)=accel(j)+d_a(j,1)
2121 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2123 ! accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
2124 accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
2125 accel(j)=accel(j)+d_a(j,i+nres)
2129 ! if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
2130 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
2131 dacc=dabs(accel(j)-accel_old(j))
2132 ! write (iout,*) "side-chain",i,dacc
2133 if (dacc.gt.amax) amax=dacc
2137 accel_old(j)=accel_old(j)+d_a_old(j,i)
2138 accel(j)=accel(j)+d_a(j,i)
2139 ! aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
2143 end subroutine max_accel
2144 !-----------------------------------------------------------------------------
2145 subroutine predict_edrift(epdrift)
2147 ! Predict the drift of the potential energy
2150 use control_data, only: lmuca
2151 ! implicit real*8 (a-h,o-z)
2152 ! include 'DIMENSIONS'
2153 ! include 'COMMON.CONTROL'
2154 ! include 'COMMON.VAR'
2155 ! include 'COMMON.MD'
2156 ! include 'COMMON.CHAIN'
2157 ! include 'COMMON.DERIV'
2158 ! include 'COMMON.GEO'
2159 ! include 'COMMON.LOCAL'
2160 ! include 'COMMON.INTERACT'
2161 ! include 'COMMON.IOUNITS'
2162 ! include 'COMMON.MUCA'
2163 real(kind=8) :: epdrift,epdriftij
2165 ! Drift of the potential energy
2171 epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
2172 if (lmuca) epdriftij=epdriftij*factor
2173 ! write (iout,*) "back",i,j,epdriftij
2174 if (epdriftij.gt.epdrift) epdrift=epdriftij
2178 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2181 dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
2182 if (lmuca) epdriftij=epdriftij*factor
2183 ! write (iout,*) "side",i,j,epdriftij
2184 if (epdriftij.gt.epdrift) epdrift=epdriftij
2188 epdrift=0.5d0*epdrift*d_time*d_time
2189 ! write (iout,*) "epdrift",epdrift
2191 end subroutine predict_edrift
2192 !-----------------------------------------------------------------------------
2193 subroutine verlet_bath
2195 ! Coupling to the thermostat by using the Berendsen algorithm
2198 ! implicit real*8 (a-h,o-z)
2199 ! include 'DIMENSIONS'
2200 ! include 'COMMON.CONTROL'
2201 ! include 'COMMON.VAR'
2202 ! include 'COMMON.MD'
2203 ! include 'COMMON.CHAIN'
2204 ! include 'COMMON.DERIV'
2205 ! include 'COMMON.GEO'
2206 ! include 'COMMON.LOCAL'
2207 ! include 'COMMON.INTERACT'
2208 ! include 'COMMON.IOUNITS'
2209 ! include 'COMMON.NAMES'
2210 real(kind=8) :: T_half,fact
2211 integer :: i,j,inres
2213 T_half=2.0d0/(dimen3*Rb)*EK
2214 fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
2215 ! write(iout,*) "T_half", T_half
2216 ! write(iout,*) "EK", EK
2217 ! write(iout,*) "fact", fact
2219 d_t(j,0)=fact*d_t(j,0)
2223 d_t(j,i)=fact*d_t(j,i)
2227 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2230 d_t(j,inres)=fact*d_t(j,inres)
2235 end subroutine verlet_bath
2236 !-----------------------------------------------------------------------------
2238 ! Set up the initial conditions of a MD simulation
2241 use control, only:tcpu
2242 !el use io_basic, only:ilen
2245 use minimm, only:minim_dc,minimize,sc_move
2246 use io_config, only:readrst
2247 use io, only:statout
2248 ! implicit real*8 (a-h,o-z)
2249 ! include 'DIMENSIONS'
2252 character(len=16) :: form
2253 integer :: IERROR,ERRCODE
2255 ! include 'COMMON.SETUP'
2256 ! include 'COMMON.CONTROL'
2257 ! include 'COMMON.VAR'
2258 ! include 'COMMON.MD'
2260 ! include 'COMMON.LANGEVIN'
2262 ! include 'COMMON.LANGEVIN.lang0'
2264 ! include 'COMMON.CHAIN'
2265 ! include 'COMMON.DERIV'
2266 ! include 'COMMON.GEO'
2267 ! include 'COMMON.LOCAL'
2268 ! include 'COMMON.INTERACT'
2269 ! include 'COMMON.IOUNITS'
2270 ! include 'COMMON.NAMES'
2271 ! include 'COMMON.REMD'
2272 real(kind=8),dimension(0:n_ene) :: energia_long,energia_short
2273 real(kind=8),dimension(3) :: vcm,incr,L
2274 real(kind=8) :: xv,sigv,lowb,highb
2275 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
2276 character(len=256) :: qstr
2279 character(len=50) :: tytul
2280 logical :: file_exist
2281 !el common /gucio/ cm
2282 integer :: i,j,ipos,iq,iw,nft_sc,iretcode,nfun,itime,ierr
2283 real(kind=8) :: etot,tt0
2287 ! write(iout,*) "d_time", d_time
2288 ! Compute the standard deviations of stochastic forces for Langevin dynamics
2289 ! if the friction coefficients do not depend on surface area
2290 if (lang.gt.0 .and. .not.surfarea) then
2292 stdforcp(i)=stdfp*dsqrt(gamp)
2295 stdforcsc(i)=stdfsc(iabs(itype(i))) &
2296 *dsqrt(gamsc(iabs(itype(i))))
2299 ! Open the pdb file for snapshotshots
2302 if (ilen(tmpdir).gt.0) &
2303 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2304 liczba(:ilen(liczba))//".pdb")
2306 file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2310 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) &
2311 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2312 liczba(:ilen(liczba))//".x")
2313 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2316 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) &
2317 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2318 liczba(:ilen(liczba))//".cx")
2319 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2325 if (ilen(tmpdir).gt.0) &
2326 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
2327 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
2329 if (ilen(tmpdir).gt.0) &
2330 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
2331 cartname=prefix(:ilen(prefix))//"_MD.cx"
2335 write (qstr,'(256(1h ))')
2338 iq = qinfrag(i,iset)*10
2339 iw = wfrag(i,iset)/100
2341 if(me.eq.king.or..not.out1file) &
2342 write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
2343 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
2348 iq = qinpair(i,iset)*10
2349 iw = wpair(i,iset)/100
2351 if(me.eq.king.or..not.out1file) &
2352 write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
2353 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
2357 ! pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
2359 ! cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
2361 ! cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
2363 ! statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
2367 if (restart1file) then
2369 inquire(file=mremd_rst_name,exist=file_exist)
2370 write (*,*) me," Before broadcast: file_exist",file_exist
2372 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,&
2375 write (*,*) me," After broadcast: file_exist",file_exist
2376 ! inquire(file=mremd_rst_name,exist=file_exist)
2377 if(me.eq.king.or..not.out1file) &
2378 write(iout,*) "Initial state read by master and distributed"
2380 if (ilen(tmpdir).gt.0) &
2381 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_' &
2382 //liczba(:ilen(liczba))//'.rst')
2383 inquire(file=rest2name,exist=file_exist)
2386 if(.not.restart1file) then
2387 if(me.eq.king.or..not.out1file) &
2388 write(iout,*) "Initial state will be read from file ",&
2389 rest2name(:ilen(rest2name))
2392 call rescale_weights(t_bath)
2394 if(me.eq.king.or..not.out1file)then
2395 if (restart1file) then
2396 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),&
2399 write(iout,*) "File ",rest2name(:ilen(rest2name)),&
2402 write(iout,*) "Initial velocities randomly generated"
2408 ! Generate initial velocities
2409 if(me.eq.king.or..not.out1file) &
2410 write(iout,*) "Initial velocities randomly generated"
2414 ! rest2name = prefix(:ilen(prefix))//'.rst'
2415 if(me.eq.king.or..not.out1file)then
2416 write (iout,*) "Initial velocities"
2418 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2419 (d_t(j,i+nres),j=1,3)
2421 ! Zeroing the total angular momentum of the system
2422 write(iout,*) "Calling the zero-angular momentum subroutine"
2425 ! Getting the potential energy and forces and velocities and accelerations
2427 ! write (iout,*) "velocity of the center of the mass:"
2428 ! write (iout,*) (vcm(j),j=1,3)
2430 d_t(j,0)=d_t(j,0)-vcm(j)
2432 ! Removing the velocity of the center of mass
2434 if(me.eq.king.or..not.out1file)then
2435 write (iout,*) "vcm right after adjustment:"
2436 write (iout,*) (vcm(j),j=1,3)
2438 if ((.not.rest).and.(indpdb.eq.0)) then
2440 if(iranconf.ne.0) then
2442 print *, 'Calling OVERLAP_SC'
2443 call overlap_sc(fail)
2446 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
2447 print *,'SC_move',nft_sc,etot
2448 if(me.eq.king.or..not.out1file) &
2449 write(iout,*) 'SC_move',nft_sc,etot
2453 print *, 'Calling MINIM_DC'
2454 call minim_dc(etot,iretcode,nfun)
2456 call geom_to_var(nvar,varia)
2457 print *,'Calling MINIMIZE.'
2458 call minimize(etot,varia,iretcode,nfun)
2459 call var_to_geom(nvar,varia)
2461 if(me.eq.king.or..not.out1file) &
2462 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2465 call chainbuild_cart
2470 kinetic_T=2.0d0/(dimen3*Rb)*EK
2471 if(me.eq.king.or..not.out1file)then
2481 call etotal(potEcomp)
2482 if (large) call enerprint(potEcomp)
2485 t_etotal=t_etotal+MPI_Wtime()-tt0
2487 t_etotal=t_etotal+tcpu()-tt0
2494 if (amax*d_time .gt. dvmax) then
2495 d_time=d_time*dvmax/amax
2496 if(me.eq.king.or..not.out1file) write (iout,*) &
2497 "Time step reduced to",d_time,&
2498 " because of too large initial acceleration."
2500 if(me.eq.king.or..not.out1file)then
2501 write(iout,*) "Potential energy and its components"
2502 call enerprint(potEcomp)
2503 ! write(iout,*) (potEcomp(i),i=0,n_ene)
2505 potE=potEcomp(0)-potEcomp(20)
2508 if (ntwe.ne.0) call statout(itime)
2509 if(me.eq.king.or..not.out1file) &
2510 write (iout,'(/a/3(a25,1pe14.5/))') "Initial:", &
2511 " Kinetic energy",EK," Potential energy",potE, &
2512 " Total energy",totE," Maximum acceleration ", &
2515 write (iout,*) "Initial coordinates"
2517 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),&
2520 write (iout,*) "Initial dC"
2522 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),&
2523 (dc(j,i+nres),j=1,3)
2525 write (iout,*) "Initial velocities"
2526 write (iout,"(13x,' backbone ',23x,' side chain')")
2528 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2529 (d_t(j,i+nres),j=1,3)
2531 write (iout,*) "Initial accelerations"
2533 ! write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2534 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),&
2535 (d_a(j,i+nres),j=1,3)
2541 d_t_old(j,i)=d_t(j,i)
2542 d_a_old(j,i)=d_a(j,i)
2544 ! write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
2553 call etotal_short(energia_short)
2554 if (large) call enerprint(potEcomp)
2557 t_eshort=t_eshort+MPI_Wtime()-tt0
2559 t_eshort=t_eshort+tcpu()-tt0
2564 if(.not.out1file .and. large) then
2565 write (iout,*) "energia_long",energia_long(0),&
2566 " energia_short",energia_short(0),&
2567 " total",energia_long(0)+energia_short(0)
2568 write (iout,*) "Initial fast-force accelerations"
2570 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
2571 (d_a(j,i+nres),j=1,3)
2574 ! 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
2577 d_a_short(j,i)=d_a(j,i)
2586 call etotal_long(energia_long)
2587 if (large) call enerprint(potEcomp)
2590 t_elong=t_elong+MPI_Wtime()-tt0
2592 t_elong=t_elong+tcpu()-tt0
2597 if(.not.out1file .and. large) then
2598 write (iout,*) "energia_long",energia_long(0)
2599 write (iout,*) "Initial slow-force accelerations"
2601 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
2602 (d_a(j,i+nres),j=1,3)
2606 t_enegrad=t_enegrad+MPI_Wtime()-tt0
2608 t_enegrad=t_enegrad+tcpu()-tt0
2612 end subroutine init_MD
2613 !-----------------------------------------------------------------------------
2614 subroutine random_vel
2616 ! implicit real*8 (a-h,o-z)
2618 use random, only:anorm_distr
2620 ! include 'DIMENSIONS'
2621 ! include 'COMMON.CONTROL'
2622 ! include 'COMMON.VAR'
2623 ! include 'COMMON.MD'
2625 ! include 'COMMON.LANGEVIN'
2627 ! include 'COMMON.LANGEVIN.lang0'
2629 ! include 'COMMON.CHAIN'
2630 ! include 'COMMON.DERIV'
2631 ! include 'COMMON.GEO'
2632 ! include 'COMMON.LOCAL'
2633 ! include 'COMMON.INTERACT'
2634 ! include 'COMMON.IOUNITS'
2635 ! include 'COMMON.NAMES'
2636 ! include 'COMMON.TIME1'
2637 real(kind=8) :: xv,sigv,lowb,highb ,Ek1
2639 real(kind=8) ,allocatable, dimension(:) :: DDU1,DDU2,DL2,DL1,xsolv,DML,rs
2640 real(kind=8) :: sumx
2642 real(kind=8) ,allocatable, dimension(:) :: rsold
2643 real (kind=8),allocatable,dimension(:,:) :: matold
2646 integer :: i,j,ii,k,ind,mark,imark
2647 ! Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
2648 ! First generate velocities in the eigenspace of the G matrix
2649 ! write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
2656 sigv=dsqrt((Rb*t_bath)/geigen(i))
2659 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
2661 write (iout,*) "i",i," ii",ii," geigen",geigen(i),&
2662 " d_t_work_new",d_t_work_new(ii)
2673 Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
2676 write (iout,*) "Ek from eigenvectors",Ek1
2677 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
2681 ! Transform velocities to UNRES coordinate space
2682 allocate (DL1(2*nres))
2683 allocate (DDU1(2*nres))
2684 allocate (DL2(2*nres))
2685 allocate (DDU2(2*nres))
2686 allocate (xsolv(2*nres))
2687 allocate (DML(2*nres))
2688 allocate (rs(2*nres))
2690 allocate (rsold(2*nres))
2691 allocate (matold(2*nres,2*nres))
2693 matold(1,1)=DMorig(1)
2694 matold(1,2)=DU1orig(1)
2695 matold(1,3)=DU2orig(1)
2696 write (*,*) DMorig(1),DU1orig(1),DU2orig(1)
2701 matold(i,j)=DMorig(i)
2702 matold(i,j-1)=DU1orig(i-1)
2704 matold(i,j-2)=DU2orig(i-2)
2712 matold(i,j+1)=DU1orig(i)
2718 matold(i,j+2)=DU2orig(i)
2722 matold(dimen,dimen)=DMorig(dimen)
2723 matold(dimen,dimen-1)=DU1orig(dimen-1)
2724 matold(dimen,dimen-2)=DU2orig(dimen-2)
2725 write(iout,*) "old gmatrix"
2726 call matout(dimen,dimen,2*nres,2*nres,matold)
2730 ! Find the ith eigenvector of the pentadiagonal inertiq matrix
2734 DML(j)=DMorig(j)-geigen(i)
2737 DML(j-1)=DMorig(j)-geigen(i)
2742 DDU1(imark-1)=DU2orig(imark-1)
2743 do j=imark+1,dimen-1
2744 DDU1(j-1)=DU1orig(j)
2752 DDU2(j)=DU2orig(j+1)
2761 write (iout,*) "DL2,DL1,DML,DDU1,DDU2"
2762 write(iout,'(10f10.5)') (DL2(k),k=3,dimen-1)
2763 write(iout,'(10f10.5)') (DL1(k),k=2,dimen-1)
2764 write(iout,'(10f10.5)') (DML(k),k=1,dimen-1)
2765 write(iout,'(10f10.5)') (DDU1(k),k=1,dimen-2)
2766 write(iout,'(10f10.5)') (DDU2(k),k=1,dimen-3)
2769 if (imark.gt.2) rs(imark-2)=-DU2orig(imark-2)
2770 if (imark.gt.1) rs(imark-1)=-DU1orig(imark-1)
2771 if (imark.lt.dimen) rs(imark)=-DU1orig(imark)
2772 if (imark.lt.dimen-1) rs(imark+1)=-DU2orig(imark)
2776 ! write (iout,*) "Vector rs"
2778 ! write (iout,*) j,rs(j)
2781 call FDIAG(dimen-1,DL2,DL1,DML,DDU1,DDU2,rs,xsolv,mark)
2788 sumx=-geigen(i)*xsolv(j)
2790 sumx=sumx+matold(j,k)*xsolv(k)
2793 sumx=sumx+matold(j,k)*xsolv(k-1)
2795 write(iout,'(i5,3f10.5)') j,sumx,rsold(j),sumx-rsold(j)
2798 sumx=-geigen(i)*xsolv(j-1)
2800 sumx=sumx+matold(j,k)*xsolv(k)
2803 sumx=sumx+matold(j,k)*xsolv(k-1)
2805 write(iout,'(i5,3f10.5)') j-1,sumx,rsold(j-1),sumx-rsold(j-1)
2809 "Solution of equations system",i," for eigenvalue",geigen(i)
2811 write(iout,'(i5,f10.5)') j,xsolv(j)
2814 do j=dimen-1,imark,-1
2819 write (iout,*) "Un-normalized eigenvector",i," for eigenvalue",geigen(i)
2821 write(iout,'(i5,f10.5)') j,xsolv(j)
2824 ! Normalize ith eigenvector
2827 sumx=sumx+xsolv(j)**2
2831 xsolv(j)=xsolv(j)/sumx
2834 write (iout,*) "Eigenvector",i," for eigenvalue",geigen(i)
2836 write(iout,'(i5,3f10.5)') j,xsolv(j)
2839 ! All done at this point for eigenvector i, exit loop
2847 write (iout,*) "Unable to find eigenvector",i
2850 ! write (iout,*) "i=",i
2852 ! write (iout,*) "k=",k
2855 ! write(iout,*) "ind",ind," ind1",3*(i-1)+k
2856 d_t_work(ind)=d_t_work(ind) &
2857 +xsolv(j)*d_t_work_new(3*(i-1)+k)
2860 enddo ! i (loop over the eigenvectors)
2863 write (iout,*) "d_t_work"
2865 write (iout,"(i5,f10.5)") i,d_t_work(i)
2870 if (itype(i).eq.10) then
2877 ! write (iout,*) "k",k," ii+k",ii+k," ii+j+k",ii+j+k,"EK1",Ek1
2878 Ek1=Ek1+0.5d0*mp*((d_t_work(ii+j+k)+d_t_work_new(ii+k))/2)**2+&
2879 0.5d0*0.25d0*IP*(d_t_work(ii+j+k)-d_t_work_new(ii+k))**2
2882 if (itype(i).ne.10) ii=ii+3
2883 write (iout,*) "i",i," itype",itype(i)," mass",msc(itype(i))
2884 write (iout,*) "ii",ii
2887 write (iout,*) "k",k," ii",ii,"EK1",EK1
2888 if (iabs(itype(i)).ne.10) Ek1=Ek1+0.5d0*Isc(iabs(itype(i)))*(d_t_work(ii)-d_t_work(ii-3))**2
2889 Ek1=Ek1+0.5d0*msc(iabs(itype(i)))*d_t_work(ii)**2
2891 write (iout,*) "i",i," ii",ii
2893 write (iout,*) "Ek from d_t_work",Ek1
2894 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
2910 d_t(k,j)=d_t_work(ind)
2913 if (itype(j).ne.10 .and. itype(j).ne.ntyp1) then
2915 d_t(k,j+nres)=d_t_work(ind)
2921 write (iout,*) "Random velocities in the Calpha,SC space"
2923 write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3)
2926 write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3)
2933 if (itype(i).eq.10) then
2935 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
2939 d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
2940 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
2945 write (iout,*)"Random velocities in the virtual-bond-vector space"
2947 write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3)
2950 write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3)
2953 write (iout,*) "Ek from d_t_work",Ek1
2954 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
2962 d_t_work(ind)=d_t_work(ind) &
2963 +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
2965 ! write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
2969 ! Transfer to the d_t vector
2971 d_t(j,0)=d_t_work(j)
2977 d_t(j,i)=d_t_work(ind)
2981 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
2984 d_t(j,i+nres)=d_t_work(ind)
2990 ! write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",&
2991 ! 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
2995 end subroutine random_vel
2996 !-----------------------------------------------------------------------------
2998 subroutine sd_verlet_p_setup
2999 ! Sets up the parameters of stochastic Verlet algorithm
3000 ! implicit real*8 (a-h,o-z)
3001 ! include 'DIMENSIONS'
3002 use control, only: tcpu
3007 ! include 'COMMON.CONTROL'
3008 ! include 'COMMON.VAR'
3009 ! include 'COMMON.MD'
3011 ! include 'COMMON.LANGEVIN'
3013 ! include 'COMMON.LANGEVIN.lang0'
3015 ! include 'COMMON.CHAIN'
3016 ! include 'COMMON.DERIV'
3017 ! include 'COMMON.GEO'
3018 ! include 'COMMON.LOCAL'
3019 ! include 'COMMON.INTERACT'
3020 ! include 'COMMON.IOUNITS'
3021 ! include 'COMMON.NAMES'
3022 ! include 'COMMON.TIME1'
3023 real(kind=8),dimension(6*nres) :: emgdt !(MAXRES6) maxres6=6*maxres
3024 real(kind=8) :: pterm,vterm,rho,rhoc,vsig
3025 real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
3026 prand_vec,vrand_vec1,vrand_vec2 !(MAXRES6) maxres6=6*maxres
3027 logical :: lprn = .false.
3028 real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
3029 real(kind=8) :: ktm,gdt,egdt,gdt2,gdt3,gdt4,gdt5,gdt6,gdt7,gdt8,&
3031 integer :: i,maxres2
3038 ! AL 8/17/04 Code adapted from tinker
3040 ! Get the frictional and random terms for stochastic dynamics in the
3041 ! eigenspace of mass-scaled UNRES friction matrix
3045 gdt = fricgam(i) * d_time
3047 ! Stochastic dynamics reduces to simple MD for zero friction
3049 if (gdt .le. zero) then
3050 pfric_vec(i) = 1.0d0
3051 vfric_vec(i) = d_time
3052 afric_vec(i) = 0.5d0 * d_time * d_time
3053 prand_vec(i) = 0.0d0
3054 vrand_vec1(i) = 0.0d0
3055 vrand_vec2(i) = 0.0d0
3057 ! Analytical expressions when friction coefficient is large
3060 if (gdt .ge. gdt_radius) then
3063 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
3064 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
3065 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
3066 vterm = 1.0d0 - egdt**2
3067 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
3069 ! Use series expansions when friction coefficient is small
3080 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0 &
3081 - gdt5/120.0d0 + gdt6/720.0d0 &
3082 - gdt7/5040.0d0 + gdt8/40320.0d0 &
3083 - gdt9/362880.0d0) / fricgam(i)**2
3084 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
3085 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
3086 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0 &
3087 + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0 &
3088 + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0 &
3089 + 127.0d0*gdt9/90720.0d0
3090 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0 &
3091 - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0 &
3092 - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0 &
3093 - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
3094 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0 &
3095 - 17.0d0*gdt2/1280.0d0 &
3096 + 17.0d0*gdt3/6144.0d0 &
3097 + 40967.0d0*gdt4/34406400.0d0 &
3098 - 57203.0d0*gdt5/275251200.0d0 &
3099 - 1429487.0d0*gdt6/13212057600.0d0)
3102 ! Compute the scaling factors of random terms for the nonzero friction case
3104 ktm = 0.5d0*d_time/fricgam(i)
3105 psig = dsqrt(ktm*pterm) / fricgam(i)
3106 vsig = dsqrt(ktm*vterm)
3107 rhoc = dsqrt(1.0d0 - rho*rho)
3109 vrand_vec1(i) = vsig * rho
3110 vrand_vec2(i) = vsig * rhoc
3115 "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",&
3118 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),&
3119 afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
3123 ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
3126 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
3127 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
3128 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
3129 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
3130 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
3131 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
3134 t_sdsetup=t_sdsetup+MPI_Wtime()
3136 t_sdsetup=t_sdsetup+tcpu()-tt0
3139 end subroutine sd_verlet_p_setup
3140 !-----------------------------------------------------------------------------
3141 subroutine eigtransf1(n,ndim,ab,d,c)
3145 real(kind=8) :: ab(ndim,ndim,n),c(ndim,n),d(ndim)
3151 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
3156 end subroutine eigtransf1
3157 !-----------------------------------------------------------------------------
3158 subroutine eigtransf(n,ndim,a,b,d,c)
3162 real(kind=8) :: a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
3168 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
3173 end subroutine eigtransf
3174 !-----------------------------------------------------------------------------
3175 subroutine sd_verlet1
3177 ! Applying stochastic velocity Verlet algorithm - step 1 to velocities
3179 ! implicit real*8 (a-h,o-z)
3180 ! include 'DIMENSIONS'
3181 ! include 'COMMON.CONTROL'
3182 ! include 'COMMON.VAR'
3183 ! include 'COMMON.MD'
3185 ! include 'COMMON.LANGEVIN'
3187 ! include 'COMMON.LANGEVIN.lang0'
3189 ! include 'COMMON.CHAIN'
3190 ! include 'COMMON.DERIV'
3191 ! include 'COMMON.GEO'
3192 ! include 'COMMON.LOCAL'
3193 ! include 'COMMON.INTERACT'
3194 ! include 'COMMON.IOUNITS'
3195 ! include 'COMMON.NAMES'
3196 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
3197 !el common /stochcalc/ stochforcvec
3198 logical :: lprn = .false.
3199 real(kind=8) :: ddt1,ddt2
3200 integer :: i,j,ind,inres
3202 ! write (iout,*) "dc_old"
3204 ! write (iout,'(i5,3f10.5,5x,3f10.5)')
3205 ! & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3208 dc_work(j)=dc_old(j,0)
3209 d_t_work(j)=d_t_old(j,0)
3210 d_a_work(j)=d_a_old(j,0)
3215 dc_work(ind+j)=dc_old(j,i)
3216 d_t_work(ind+j)=d_t_old(j,i)
3217 d_a_work(ind+j)=d_a_old(j,i)
3222 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3224 dc_work(ind+j)=dc_old(j,i+nres)
3225 d_t_work(ind+j)=d_t_old(j,i+nres)
3226 d_a_work(ind+j)=d_a_old(j,i+nres)
3234 "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",&
3238 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
3239 vfric_mat(i,j),afric_mat(i,j),&
3240 prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3248 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) &
3249 +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3250 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3251 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3253 d_t_work_new(i)=ddt1+0.5d0*ddt2
3254 d_t_work(i)=ddt1+ddt2
3259 d_t(j,0)=d_t_work(j)
3264 dc(j,i)=dc_work(ind+j)
3265 d_t(j,i)=d_t_work(ind+j)
3270 if (itype(i).ne.10) then
3273 dc(j,inres)=dc_work(ind+j)
3274 d_t(j,inres)=d_t_work(ind+j)
3280 end subroutine sd_verlet1
3281 !-----------------------------------------------------------------------------
3282 subroutine sd_verlet2
3284 ! Calculating the adjusted velocities for accelerations
3286 ! implicit real*8 (a-h,o-z)
3287 ! include 'DIMENSIONS'
3288 ! include 'COMMON.CONTROL'
3289 ! include 'COMMON.VAR'
3290 ! include 'COMMON.MD'
3292 ! include 'COMMON.LANGEVIN'
3294 ! include 'COMMON.LANGEVIN.lang0'
3296 ! include 'COMMON.CHAIN'
3297 ! include 'COMMON.DERIV'
3298 ! include 'COMMON.GEO'
3299 ! include 'COMMON.LOCAL'
3300 ! include 'COMMON.INTERACT'
3301 ! include 'COMMON.IOUNITS'
3302 ! include 'COMMON.NAMES'
3303 !el real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV !(MAXRES6) maxres6=6*maxres
3304 real(kind=8),dimension(6*nres) :: stochforcvecV !(MAXRES6) maxres6=6*maxres
3305 !el common /stochcalc/ stochforcvec
3307 real(kind=8) :: ddt1,ddt2
3308 integer :: i,j,ind,inres
3309 ! Compute the stochastic forces which contribute to velocity change
3311 call stochastic_force(stochforcvecV)
3318 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3319 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+ &
3320 vrand_mat2(i,j)*stochforcvecV(j)
3322 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
3326 d_t(j,0)=d_t_work(j)
3331 d_t(j,i)=d_t_work(ind+j)
3336 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3339 d_t(j,inres)=d_t_work(ind+j)
3345 end subroutine sd_verlet2
3346 !-----------------------------------------------------------------------------
3347 subroutine sd_verlet_ciccotti_setup
3349 ! Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
3351 ! implicit real*8 (a-h,o-z)
3352 ! include 'DIMENSIONS'
3353 use control, only: tcpu
3358 ! include 'COMMON.CONTROL'
3359 ! include 'COMMON.VAR'
3360 ! include 'COMMON.MD'
3362 ! include 'COMMON.LANGEVIN'
3364 ! include 'COMMON.LANGEVIN.lang0'
3366 ! include 'COMMON.CHAIN'
3367 ! include 'COMMON.DERIV'
3368 ! include 'COMMON.GEO'
3369 ! include 'COMMON.LOCAL'
3370 ! include 'COMMON.INTERACT'
3371 ! include 'COMMON.IOUNITS'
3372 ! include 'COMMON.NAMES'
3373 ! include 'COMMON.TIME1'
3374 real(kind=8),dimension(6*nres) :: emgdt !(MAXRES6) maxres6=6*maxres
3375 real(kind=8) :: pterm,vterm,rho,rhoc,vsig
3376 real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
3377 prand_vec,vrand_vec1,vrand_vec2 !(MAXRES6) maxres6=6*maxres
3378 logical :: lprn = .false.
3379 real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
3380 real(kind=8) :: ktm,gdt,egdt,tt0
3381 integer :: i,maxres2
3388 ! AL 8/17/04 Code adapted from tinker
3390 ! Get the frictional and random terms for stochastic dynamics in the
3391 ! eigenspace of mass-scaled UNRES friction matrix
3395 write (iout,*) "i",i," fricgam",fricgam(i)
3396 gdt = fricgam(i) * d_time
3398 ! Stochastic dynamics reduces to simple MD for zero friction
3400 if (gdt .le. zero) then
3401 pfric_vec(i) = 1.0d0
3402 vfric_vec(i) = d_time
3403 afric_vec(i) = 0.5d0*d_time*d_time
3404 prand_vec(i) = afric_vec(i)
3405 vrand_vec2(i) = vfric_vec(i)
3407 ! Analytical expressions when friction coefficient is large
3412 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
3413 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
3414 prand_vec(i) = afric_vec(i)
3415 vrand_vec2(i) = vfric_vec(i)
3417 ! Compute the scaling factors of random terms for the nonzero friction case
3419 ! ktm = 0.5d0*d_time/fricgam(i)
3420 ! psig = dsqrt(ktm*pterm) / fricgam(i)
3421 ! vsig = dsqrt(ktm*vterm)
3422 ! prand_vec(i) = psig*afric_vec(i)
3423 ! vrand_vec2(i) = vsig*vfric_vec(i)
3428 "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",&
3431 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),&
3432 afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
3436 ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
3438 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
3439 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
3440 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
3441 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
3442 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
3444 t_sdsetup=t_sdsetup+MPI_Wtime()
3446 t_sdsetup=t_sdsetup+tcpu()-tt0
3449 end subroutine sd_verlet_ciccotti_setup
3450 !-----------------------------------------------------------------------------
3451 subroutine sd_verlet1_ciccotti
3453 ! Applying stochastic velocity Verlet algorithm - step 1 to velocities
3454 ! implicit real*8 (a-h,o-z)
3456 ! include 'DIMENSIONS'
3460 ! include 'COMMON.CONTROL'
3461 ! include 'COMMON.VAR'
3462 ! include 'COMMON.MD'
3464 ! include 'COMMON.LANGEVIN'
3466 ! include 'COMMON.LANGEVIN.lang0'
3468 ! include 'COMMON.CHAIN'
3469 ! include 'COMMON.DERIV'
3470 ! include 'COMMON.GEO'
3471 ! include 'COMMON.LOCAL'
3472 ! include 'COMMON.INTERACT'
3473 ! include 'COMMON.IOUNITS'
3474 ! include 'COMMON.NAMES'
3475 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
3476 !el common /stochcalc/ stochforcvec
3477 logical :: lprn = .false.
3478 real(kind=8) :: ddt1,ddt2
3479 integer :: i,j,ind,inres
3480 ! write (iout,*) "dc_old"
3482 ! write (iout,'(i5,3f10.5,5x,3f10.5)')
3483 ! & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3486 dc_work(j)=dc_old(j,0)
3487 d_t_work(j)=d_t_old(j,0)
3488 d_a_work(j)=d_a_old(j,0)
3493 dc_work(ind+j)=dc_old(j,i)
3494 d_t_work(ind+j)=d_t_old(j,i)
3495 d_a_work(ind+j)=d_a_old(j,i)
3500 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3502 dc_work(ind+j)=dc_old(j,i+nres)
3503 d_t_work(ind+j)=d_t_old(j,i+nres)
3504 d_a_work(ind+j)=d_a_old(j,i+nres)
3513 "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",&
3517 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
3518 vfric_mat(i,j),afric_mat(i,j),&
3519 prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3527 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) &
3528 +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3529 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3530 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3532 d_t_work_new(i)=ddt1+0.5d0*ddt2
3533 d_t_work(i)=ddt1+ddt2
3538 d_t(j,0)=d_t_work(j)
3543 dc(j,i)=dc_work(ind+j)
3544 d_t(j,i)=d_t_work(ind+j)
3549 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3552 dc(j,inres)=dc_work(ind+j)
3553 d_t(j,inres)=d_t_work(ind+j)
3559 end subroutine sd_verlet1_ciccotti
3560 !-----------------------------------------------------------------------------
3561 subroutine sd_verlet2_ciccotti
3563 ! Calculating the adjusted velocities for accelerations
3565 ! implicit real*8 (a-h,o-z)
3566 ! include 'DIMENSIONS'
3567 ! include 'COMMON.CONTROL'
3568 ! include 'COMMON.VAR'
3569 ! include 'COMMON.MD'
3571 ! include 'COMMON.LANGEVIN'
3573 ! include 'COMMON.LANGEVIN.lang0'
3575 ! include 'COMMON.CHAIN'
3576 ! include 'COMMON.DERIV'
3577 ! include 'COMMON.GEO'
3578 ! include 'COMMON.LOCAL'
3579 ! include 'COMMON.INTERACT'
3580 ! include 'COMMON.IOUNITS'
3581 ! include 'COMMON.NAMES'
3582 !el real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV !(MAXRES6) maxres6=6*maxres
3583 real(kind=8),dimension(6*nres) :: stochforcvecV !(MAXRES6) maxres6=6*maxres
3584 !el common /stochcalc/ stochforcvec
3585 real(kind=8) :: ddt1,ddt2
3586 integer :: i,j,ind,inres
3588 ! Compute the stochastic forces which contribute to velocity change
3590 call stochastic_force(stochforcvecV)
3597 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3598 ! ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
3599 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
3601 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
3605 d_t(j,0)=d_t_work(j)
3610 d_t(j,i)=d_t_work(ind+j)
3615 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3618 d_t(j,inres)=d_t_work(ind+j)
3624 end subroutine sd_verlet2_ciccotti
3626 !-----------------------------------------------------------------------------
3628 !-----------------------------------------------------------------------------
3629 subroutine inertia_tensor
3631 ! Calculating the intertia tensor for the entire protein in order to
3632 ! remove the perpendicular components of velocity matrix which cause
3633 ! the molecule to rotate.
3636 ! implicit real*8 (a-h,o-z)
3637 ! include 'DIMENSIONS'
3638 ! include 'COMMON.CONTROL'
3639 ! include 'COMMON.VAR'
3640 ! include 'COMMON.MD'
3641 ! include 'COMMON.CHAIN'
3642 ! include 'COMMON.DERIV'
3643 ! include 'COMMON.GEO'
3644 ! include 'COMMON.LOCAL'
3645 ! include 'COMMON.INTERACT'
3646 ! include 'COMMON.IOUNITS'
3647 ! include 'COMMON.NAMES'
3649 real(kind=8),dimension(3,3) :: Im,Imcp,eigvec,Id
3650 real(kind=8),dimension(3) :: pr,eigval,L,vp,vrot
3651 real(kind=8) :: M_SC,mag,mag2
3652 real(kind=8),dimension(3,0:nres) :: vpp !(3,0:MAXRES)
3653 real(kind=8),dimension(3) :: vs_p,pp,incr,v
3654 real(kind=8),dimension(3,3) :: pr1,pr2
3656 !el common /gucio/ cm
3657 integer :: iti,inres,i,j,k
3668 ! calculating the center of the mass of the protein
3671 cm(j)=cm(j)+c(j,i)+0.5d0*dc(j,i)
3680 M_SC=M_SC+msc(iabs(iti))
3683 cm(j)=cm(j)+msc(iabs(iti))*c(j,inres)
3687 cm(j)=cm(j)/(M_SC+(nct-nnt)*mp)
3692 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
3694 Im(1,1)=Im(1,1)+mp*(pr(2)*pr(2)+pr(3)*pr(3))
3695 Im(1,2)=Im(1,2)-mp*pr(1)*pr(2)
3696 Im(1,3)=Im(1,3)-mp*pr(1)*pr(3)
3697 Im(2,3)=Im(2,3)-mp*pr(2)*pr(3)
3698 Im(2,2)=Im(2,2)+mp*(pr(3)*pr(3)+pr(1)*pr(1))
3699 Im(3,3)=Im(3,3)+mp*(pr(1)*pr(1)+pr(2)*pr(2))
3706 pr(j)=c(j,inres)-cm(j)
3708 Im(1,1)=Im(1,1)+msc(iabs(iti))*(pr(2)*pr(2)+pr(3)*pr(3))
3709 Im(1,2)=Im(1,2)-msc(iabs(iti))*pr(1)*pr(2)
3710 Im(1,3)=Im(1,3)-msc(iabs(iti))*pr(1)*pr(3)
3711 Im(2,3)=Im(2,3)-msc(iabs(iti))*pr(2)*pr(3)
3712 Im(2,2)=Im(2,2)+msc(iabs(iti))*(pr(3)*pr(3)+pr(1)*pr(1))
3713 Im(3,3)=Im(3,3)+msc(iabs(iti))*(pr(1)*pr(1)+pr(2)*pr(2))
3717 Im(1,1)=Im(1,1)+Ip*(1-dc_norm(1,i)*dc_norm(1,i))* &
3718 vbld(i+1)*vbld(i+1)*0.25d0
3719 Im(1,2)=Im(1,2)+Ip*(-dc_norm(1,i)*dc_norm(2,i))* &
3720 vbld(i+1)*vbld(i+1)*0.25d0
3721 Im(1,3)=Im(1,3)+Ip*(-dc_norm(1,i)*dc_norm(3,i))* &
3722 vbld(i+1)*vbld(i+1)*0.25d0
3723 Im(2,3)=Im(2,3)+Ip*(-dc_norm(2,i)*dc_norm(3,i))* &
3724 vbld(i+1)*vbld(i+1)*0.25d0
3725 Im(2,2)=Im(2,2)+Ip*(1-dc_norm(2,i)*dc_norm(2,i))* &
3726 vbld(i+1)*vbld(i+1)*0.25d0
3727 Im(3,3)=Im(3,3)+Ip*(1-dc_norm(3,i)*dc_norm(3,i))* &
3728 vbld(i+1)*vbld(i+1)*0.25d0
3733 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3736 Im(1,1)=Im(1,1)+Isc(iti)*(1-dc_norm(1,inres)* &
3737 dc_norm(1,inres))*vbld(inres)*vbld(inres)
3738 Im(1,2)=Im(1,2)-Isc(iti)*(dc_norm(1,inres)* &
3739 dc_norm(2,inres))*vbld(inres)*vbld(inres)
3740 Im(1,3)=Im(1,3)-Isc(iti)*(dc_norm(1,inres)* &
3741 dc_norm(3,inres))*vbld(inres)*vbld(inres)
3742 Im(2,3)=Im(2,3)-Isc(iti)*(dc_norm(2,inres)* &
3743 dc_norm(3,inres))*vbld(inres)*vbld(inres)
3744 Im(2,2)=Im(2,2)+Isc(iti)*(1-dc_norm(2,inres)* &
3745 dc_norm(2,inres))*vbld(inres)*vbld(inres)
3746 Im(3,3)=Im(3,3)+Isc(iti)*(1-dc_norm(3,inres)* &
3747 dc_norm(3,inres))*vbld(inres)*vbld(inres)
3752 ! write(iout,*) "The angular momentum before adjustment:"
3753 ! write(iout,*) (L(j),j=1,3)
3759 ! Copying the Im matrix for the djacob subroutine
3767 ! Finding the eigenvectors and eignvalues of the inertia tensor
3768 call djacob(3,3,10000,1.0d-10,Imcp,eigvec,eigval)
3769 ! write (iout,*) "Eigenvalues & Eigenvectors"
3770 ! write (iout,'(5x,3f10.5)') (eigval(i),i=1,3)
3773 ! write (iout,'(i5,3f10.5)') i,(eigvec(i,j),j=1,3)
3775 ! Constructing the diagonalized matrix
3777 if (dabs(eigval(i)).gt.1.0d-15) then
3778 Id(i,i)=1.0d0/eigval(i)
3785 Imcp(i,j)=eigvec(j,i)
3791 pr1(i,j)=pr1(i,j)+Id(i,k)*Imcp(k,j)
3798 pr2(i,j)=pr2(i,j)+eigvec(i,k)*pr1(k,j)
3802 ! Calculating the total rotational velocity of the molecule
3805 vrot(i)=vrot(i)+pr2(i,j)*L(j)
3808 ! Resetting the velocities
3810 call vecpr(vrot(1),dc(1,i),vp)
3812 d_t(j,i)=d_t(j,i)-vp(j)
3816 if(itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3818 call vecpr(vrot(1),dc(1,inres),vp)
3820 d_t(j,inres)=d_t(j,inres)-vp(j)
3825 ! write(iout,*) "The angular momentum after adjustment:"
3826 ! write(iout,*) (L(j),j=1,3)
3829 end subroutine inertia_tensor
3830 !-----------------------------------------------------------------------------
3831 subroutine angmom(cm,L)
3834 ! implicit real*8 (a-h,o-z)
3835 ! include 'DIMENSIONS'
3836 ! include 'COMMON.CONTROL'
3837 ! include 'COMMON.VAR'
3838 ! include 'COMMON.MD'
3839 ! include 'COMMON.CHAIN'
3840 ! include 'COMMON.DERIV'
3841 ! include 'COMMON.GEO'
3842 ! include 'COMMON.LOCAL'
3843 ! include 'COMMON.INTERACT'
3844 ! include 'COMMON.IOUNITS'
3845 ! include 'COMMON.NAMES'
3847 real(kind=8),dimension(3) :: L,cm,pr,vp,vrot,incr,v,pp
3848 integer :: iti,inres,i,j
3849 ! Calculate the angular momentum
3858 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
3861 v(j)=incr(j)+0.5d0*d_t(j,i)
3864 incr(j)=incr(j)+d_t(j,i)
3866 call vecpr(pr(1),v(1),vp)
3872 pp(j)=0.5d0*d_t(j,i)
3874 call vecpr(pr(1),pp(1),vp)
3886 pr(j)=c(j,inres)-cm(j)
3888 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3890 v(j)=incr(j)+d_t(j,inres)
3897 call vecpr(pr(1),v(1),vp)
3898 ! write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),
3899 ! & " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
3901 L(j)=L(j)+msc(iabs(iti))*vp(j)
3903 ! write (iout,*) "L",(l(j),j=1,3)
3904 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3906 v(j)=incr(j)+d_t(j,inres)
3908 call vecpr(dc(1,inres),d_t(1,inres),vp)
3910 L(j)=L(j)+Isc(iti)*vp(j)
3914 incr(j)=incr(j)+d_t(j,i)
3918 end subroutine angmom
3919 !-----------------------------------------------------------------------------
3920 subroutine vcm_vel(vcm)
3923 ! implicit real*8 (a-h,o-z)
3924 ! include 'DIMENSIONS'
3925 ! include 'COMMON.VAR'
3926 ! include 'COMMON.MD'
3927 ! include 'COMMON.CHAIN'
3928 ! include 'COMMON.DERIV'
3929 ! include 'COMMON.GEO'
3930 ! include 'COMMON.LOCAL'
3931 ! include 'COMMON.INTERACT'
3932 ! include 'COMMON.IOUNITS'
3933 real(kind=8),dimension(3) :: vcm,vv
3934 real(kind=8) :: summas,amas
3946 vcm(j)=vcm(j)+mp*(vv(j)+0.5d0*d_t(j,i))
3949 amas=msc(iabs(itype(i)))
3951 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
3953 vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres))
3957 vcm(j)=vcm(j)+amas*vv(j)
3961 vv(j)=vv(j)+d_t(j,i)
3964 ! write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas
3966 vcm(j)=vcm(j)/summas
3969 end subroutine vcm_vel
3970 !-----------------------------------------------------------------------------
3972 !-----------------------------------------------------------------------------
3974 ! RATTLE algorithm for velocity Verlet - step 1, UNRES
3978 ! implicit real*8 (a-h,o-z)
3979 ! include 'DIMENSIONS'
3981 ! include 'COMMON.CONTROL'
3982 ! include 'COMMON.VAR'
3983 ! include 'COMMON.MD'
3985 ! include 'COMMON.LANGEVIN'
3987 ! include 'COMMON.LANGEVIN.lang0'
3989 ! include 'COMMON.CHAIN'
3990 ! include 'COMMON.DERIV'
3991 ! include 'COMMON.GEO'
3992 ! include 'COMMON.LOCAL'
3993 ! include 'COMMON.INTERACT'
3994 ! include 'COMMON.IOUNITS'
3995 ! include 'COMMON.NAMES'
3996 ! include 'COMMON.TIME1'
3997 !el real(kind=8) :: gginv(2*nres,2*nres),&
3998 !el gdc(3,2*nres,2*nres)
3999 real(kind=8) :: dC_uncor(3,2*nres) !,&
4000 !el real(kind=8) :: Cmat(2*nres,2*nres)
4001 real(kind=8) :: x(2*nres),xcorr(3,2*nres) !maxres2=2*maxres
4002 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
4003 !el common /przechowalnia/ nbond
4004 integer :: max_rattle = 5
4005 logical :: lprn = .false., lprn1 = .false., not_done
4006 real(kind=8) :: tol_rattle = 1.0d-5
4008 integer :: ii,i,j,jj,l,ind,ind1,nres2
4011 !el /common/ przechowalnia
4013 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4014 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4015 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4017 if (lprn) write (iout,*) "RATTLE1"
4020 if (itype(i).ne.10) nbond=nbond+1
4022 ! Make a folded form of the Ginv-matrix
4035 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4039 if (itype(k).ne.10) then
4043 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4050 if (itype(i).ne.10) then
4060 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4064 if (itype(k).ne.10) then
4068 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4076 write (iout,*) "Matrix GGinv"
4077 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
4083 if (iter.gt.max_rattle) then
4084 write (iout,*) "Error - too many iterations in RATTLE."
4087 ! Calculate the matrix C = GG**(-1) dC_old o dC
4092 dC_uncor(j,ind1)=dC(j,i)
4096 if (itype(i).ne.10) then
4099 dC_uncor(j,ind1)=dC(j,i+nres)
4108 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
4112 if (itype(k).ne.10) then
4115 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
4120 ! Calculate deviations from standard virtual-bond lengths
4124 x(ind)=vbld(i+1)**2-vbl**2
4127 if (itype(i).ne.10) then
4129 x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i))**2
4133 write (iout,*) "Coordinates and violations"
4135 write(iout,'(i5,3f10.5,5x,e15.5)') &
4136 i,(dC_uncor(j,i),j=1,3),x(i)
4138 write (iout,*) "Velocities and violations"
4142 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4143 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4146 if (itype(i).ne.10) then
4148 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4149 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4150 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4153 ! write (iout,*) "gdc"
4155 ! write (iout,*) "i",i
4157 ! write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4163 if (dabs(x(i)).gt.xmax) then
4167 if (xmax.lt.tol_rattle) then
4171 ! Calculate the matrix of the system of equations
4176 Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
4181 write (iout,*) "Matrix Cmat"
4182 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4184 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4185 ! Add constraint term to positions
4192 xx = xx+x(ii)*gdc(j,ind,ii)
4196 d_t_new(j,i)=d_t_new(j,i)-xx/d_time
4200 if (itype(i).ne.10) then
4205 xx = xx+x(ii)*gdc(j,ind,ii)
4208 dC(j,i+nres)=dC(j,i+nres)-xx
4209 d_t_new(j,i+nres)=d_t_new(j,i+nres)-xx/d_time
4213 ! Rebuild the chain using the new coordinates
4214 call chainbuild_cart
4216 write (iout,*) "New coordinates, Lagrange multipliers,",&
4217 " and differences between actual and standard bond lengths"
4221 xx=vbld(i+1)**2-vbl**2
4222 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4223 i,(dC(j,i),j=1,3),x(ind),xx
4226 if (itype(i).ne.10) then
4228 xx=vbld(i+nres)**2-vbldsc0(1,itype(i))**2
4229 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4230 i,(dC(j,i+nres),j=1,3),x(ind),xx
4233 write (iout,*) "Velocities and violations"
4237 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4238 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4241 if (itype(i).ne.10) then
4243 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4244 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4245 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4252 10 write (iout,*) "Error - singularity in solving the system",&
4253 " of equations for Lagrange multipliers."
4257 "RATTLE inactive; use -DRATTLE switch at compile time."
4260 end subroutine rattle1
4261 !-----------------------------------------------------------------------------
4263 ! RATTLE algorithm for velocity Verlet - step 2, UNRES
4267 ! implicit real*8 (a-h,o-z)
4268 ! include 'DIMENSIONS'
4270 ! include 'COMMON.CONTROL'
4271 ! include 'COMMON.VAR'
4272 ! include 'COMMON.MD'
4274 ! include 'COMMON.LANGEVIN'
4276 ! include 'COMMON.LANGEVIN.lang0'
4278 ! include 'COMMON.CHAIN'
4279 ! include 'COMMON.DERIV'
4280 ! include 'COMMON.GEO'
4281 ! include 'COMMON.LOCAL'
4282 ! include 'COMMON.INTERACT'
4283 ! include 'COMMON.IOUNITS'
4284 ! include 'COMMON.NAMES'
4285 ! include 'COMMON.TIME1'
4286 !el real(kind=8) :: gginv(2*nres,2*nres),&
4287 !el gdc(3,2*nres,2*nres)
4288 real(kind=8) :: dC_uncor(3,2*nres) !,&
4289 !el Cmat(2*nres,2*nres)
4290 real(kind=8) :: x(2*nres) !maxres2=2*maxres
4291 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
4292 !el common /przechowalnia/ nbond
4293 integer :: max_rattle = 5
4294 logical :: lprn = .false., lprn1 = .false., not_done
4295 real(kind=8) :: tol_rattle = 1.0d-5
4299 !el /common/ przechowalnia
4300 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4301 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4302 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4304 if (lprn) write (iout,*) "RATTLE2"
4305 if (lprn) write (iout,*) "Velocity correction"
4306 ! Calculate the matrix G dC
4312 gdc(j,i,ind)=GGinv(i,ind)*dC(j,k)
4316 if (itype(k).ne.10) then
4319 gdc(j,i,ind)=GGinv(i,ind)*dC(j,k+nres)
4325 ! write (iout,*) "gdc"
4327 ! write (iout,*) "i",i
4329 ! write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4333 ! Calculate the matrix of the system of equations
4340 Cmat(ind,j)=Cmat(ind,j)+dC(k,i)*gdc(k,ind,j)
4345 if (itype(i).ne.10) then
4350 Cmat(ind,j)=Cmat(ind,j)+dC(k,i+nres)*gdc(k,ind,j)
4355 ! Calculate the scalar product dC o d_t_new
4359 x(ind)=scalar(d_t(1,i),dC(1,i))
4362 if (itype(i).ne.10) then
4364 x(ind)=scalar(d_t(1,i+nres),dC(1,i+nres))
4368 write (iout,*) "Velocities and violations"
4372 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4373 i,ind,(d_t(j,i),j=1,3),x(ind)
4376 if (itype(i).ne.10) then
4378 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4379 i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind)
4385 if (dabs(x(i)).gt.xmax) then
4389 if (xmax.lt.tol_rattle) then
4394 write (iout,*) "Matrix Cmat"
4395 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4397 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4398 ! Add constraint term to velocities
4405 xx = xx+x(ii)*gdc(j,ind,ii)
4407 d_t(j,i)=d_t(j,i)-xx
4411 if (itype(i).ne.10) then
4416 xx = xx+x(ii)*gdc(j,ind,ii)
4418 d_t(j,i+nres)=d_t(j,i+nres)-xx
4424 "New velocities, Lagrange multipliers violations"
4428 if (lprn) write (iout,'(2i5,3f10.5,5x,2e15.5)') &
4429 i,ind,(d_t(j,i),j=1,3),x(ind),scalar(d_t(1,i),dC(1,i))
4432 if (itype(i).ne.10) then
4434 write (iout,'(2i5,3f10.5,5x,2e15.5)') &
4435 i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind),&
4436 scalar(d_t(1,i+nres),dC(1,i+nres))
4442 10 write (iout,*) "Error - singularity in solving the system",&
4443 " of equations for Lagrange multipliers."
4447 "RATTLE inactive; use -DRATTLE option at compile time."
4450 end subroutine rattle2
4451 !-----------------------------------------------------------------------------
4452 subroutine rattle_brown
4453 ! RATTLE/LINCS algorithm for Brownian dynamics, UNRES
4457 ! implicit real*8 (a-h,o-z)
4458 ! include 'DIMENSIONS'
4460 ! include 'COMMON.CONTROL'
4461 ! include 'COMMON.VAR'
4462 ! include 'COMMON.MD'
4464 ! include 'COMMON.LANGEVIN'
4466 ! include 'COMMON.LANGEVIN.lang0'
4468 ! include 'COMMON.CHAIN'
4469 ! include 'COMMON.DERIV'
4470 ! include 'COMMON.GEO'
4471 ! include 'COMMON.LOCAL'
4472 ! include 'COMMON.INTERACT'
4473 ! include 'COMMON.IOUNITS'
4474 ! include 'COMMON.NAMES'
4475 ! include 'COMMON.TIME1'
4476 !el real(kind=8) :: gginv(2*nres,2*nres),&
4477 !el gdc(3,2*nres,2*nres)
4478 real(kind=8) :: dC_uncor(3,2*nres) !,&
4479 !el real(kind=8) :: Cmat(2*nres,2*nres)
4480 real(kind=8) :: x(2*nres) !maxres2=2*maxres
4481 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
4482 !el common /przechowalnia/ nbond
4483 integer :: max_rattle = 5
4484 logical :: lprn = .true., lprn1 = .true., not_done
4485 real(kind=8) :: tol_rattle = 1.0d-5
4489 !el /common/ przechowalnia
4490 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4491 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4492 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4495 if (lprn) write (iout,*) "RATTLE_BROWN"
4498 if (itype(i).ne.10) nbond=nbond+1
4500 ! Make a folded form of the Ginv-matrix
4513 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4517 if (itype(k).ne.10) then
4521 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4528 if (itype(i).ne.10) then
4538 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4542 if (itype(k).ne.10) then
4546 if (j.eq.1 .and. l.eq.1)GGinv(ii,jj)=fricmat(ind,ind1)
4554 write (iout,*) "Matrix GGinv"
4555 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
4561 if (iter.gt.max_rattle) then
4562 write (iout,*) "Error - too many iterations in RATTLE."
4565 ! Calculate the matrix C = GG**(-1) dC_old o dC
4570 dC_uncor(j,ind1)=dC(j,i)
4574 if (itype(i).ne.10) then
4577 dC_uncor(j,ind1)=dC(j,i+nres)
4586 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
4590 if (itype(k).ne.10) then
4593 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
4598 ! Calculate deviations from standard virtual-bond lengths
4602 x(ind)=vbld(i+1)**2-vbl**2
4605 if (itype(i).ne.10) then
4607 x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i))**2
4611 write (iout,*) "Coordinates and violations"
4613 write(iout,'(i5,3f10.5,5x,e15.5)') &
4614 i,(dC_uncor(j,i),j=1,3),x(i)
4616 write (iout,*) "Velocities and violations"
4620 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4621 i,ind,(d_t(j,i),j=1,3),scalar(d_t(1,i),dC_old(1,i))
4624 if (itype(i).ne.10) then
4626 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4627 i+nres,ind,(d_t(j,i+nres),j=1,3),&
4628 scalar(d_t(1,i+nres),dC_old(1,i+nres))
4631 write (iout,*) "gdc"
4633 write (iout,*) "i",i
4635 write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4641 if (dabs(x(i)).gt.xmax) then
4645 if (xmax.lt.tol_rattle) then
4649 ! Calculate the matrix of the system of equations
4654 Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
4659 write (iout,*) "Matrix Cmat"
4660 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4662 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4663 ! Add constraint term to positions
4670 xx = xx+x(ii)*gdc(j,ind,ii)
4673 d_t(j,i)=d_t(j,i)+xx/d_time
4678 if (itype(i).ne.10) then
4683 xx = xx+x(ii)*gdc(j,ind,ii)
4686 d_t(j,i+nres)=d_t(j,i+nres)+xx/d_time
4687 dC(j,i+nres)=dC(j,i+nres)+xx
4691 ! Rebuild the chain using the new coordinates
4692 call chainbuild_cart
4694 write (iout,*) "New coordinates, Lagrange multipliers,",&
4695 " and differences between actual and standard bond lengths"
4699 xx=vbld(i+1)**2-vbl**2
4700 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4701 i,(dC(j,i),j=1,3),x(ind),xx
4704 if (itype(i).ne.10) then
4706 xx=vbld(i+nres)**2-vbldsc0(1,itype(i))**2
4707 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4708 i,(dC(j,i+nres),j=1,3),x(ind),xx
4711 write (iout,*) "Velocities and violations"
4715 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4716 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4719 if (itype(i).ne.10) then
4721 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4722 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4723 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4730 10 write (iout,*) "Error - singularity in solving the system",&
4731 " of equations for Lagrange multipliers."
4735 "RATTLE inactive; use -DRATTLE option at compile time"
4738 end subroutine rattle_brown
4739 !-----------------------------------------------------------------------------
4741 !-----------------------------------------------------------------------------
4742 subroutine friction_force
4747 ! implicit real*8 (a-h,o-z)
4748 ! include 'DIMENSIONS'
4749 ! include 'COMMON.VAR'
4750 ! include 'COMMON.CHAIN'
4751 ! include 'COMMON.DERIV'
4752 ! include 'COMMON.GEO'
4753 ! include 'COMMON.LOCAL'
4754 ! include 'COMMON.INTERACT'
4755 ! include 'COMMON.MD'
4757 ! include 'COMMON.LANGEVIN'
4759 ! include 'COMMON.LANGEVIN.lang0'
4761 ! include 'COMMON.IOUNITS'
4762 !el real(kind=8),dimension(6*nres) :: gamvec !(MAXRES6) maxres6=6*maxres
4763 !el common /syfek/ gamvec
4764 real(kind=8) :: vv(3),vvtot(3,nres),v_work(6*nres) !,&
4765 !el ginvfric(2*nres,2*nres) !maxres2=2*maxres
4766 !el common /przechowalnia/ ginvfric
4768 logical :: lprn = .false., checkmode = .false.
4769 integer :: i,j,ind,k,nres2,nres6
4773 if(.not.allocated(gamvec)) allocate(gamvec(nres6)) !(MAXRES6)
4774 if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
4782 d_t_work(j)=d_t(j,0)
4787 d_t_work(ind+j)=d_t(j,i)
4792 if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
4794 d_t_work(ind+j)=d_t(j,i+nres)
4800 call fricmat_mult(d_t_work,fric_work)
4802 if (.not.checkmode) return
4805 write (iout,*) "d_t_work and fric_work"
4807 write (iout,'(i3,2e15.5)') i,d_t_work(i),fric_work(i)
4811 friction(j,0)=fric_work(j)
4816 friction(j,i)=fric_work(ind+j)
4821 if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
4823 friction(j,i+nres)=fric_work(ind+j)
4829 write(iout,*) "Friction backbone"
4831 write(iout,'(i5,3e15.5,5x,3e15.5)') &
4832 i,(friction(j,i),j=1,3),(d_t(j,i),j=1,3)
4834 write(iout,*) "Friction side chain"
4836 write(iout,'(i5,3e15.5,5x,3e15.5)') &
4837 i,(friction(j,i+nres),j=1,3),(d_t(j,i+nres),j=1,3)
4846 vvtot(j,i)=vv(j)+0.5d0*d_t(j,i)
4847 vvtot(j,i+nres)=vv(j)+d_t(j,i+nres)
4848 vv(j)=vv(j)+d_t(j,i)
4851 write (iout,*) "vvtot backbone and sidechain"
4853 write (iout,'(i5,3e15.5,5x,3e15.5)') i,(vvtot(j,i),j=1,3),&
4854 (vvtot(j,i+nres),j=1,3)
4859 v_work(ind+j)=vvtot(j,i)
4865 v_work(ind+j)=vvtot(j,i+nres)
4869 write (iout,*) "v_work gamvec and site-based friction forces"
4871 write (iout,'(i5,3e15.5)') i,v_work(i),gamvec(i),&
4875 ! fric_work1(i)=0.0d0
4877 ! fric_work1(i)=fric_work1(i)-A(j,i)*gamvec(j)*v_work(j)
4880 ! write (iout,*) "fric_work and fric_work1"
4882 ! write (iout,'(i5,2e15.5)') i,fric_work(i),fric_work1(i)
4888 ginvfric(i,j)=ginvfric(i,j)+ginv(i,k)*fricmat(k,j)
4892 write (iout,*) "ginvfric"
4894 write (iout,'(i5,100f8.3)') i,(ginvfric(i,j),j=1,dimen)
4896 write (iout,*) "symmetry check"
4899 write (iout,*) i,j,ginvfric(i,j)-ginvfric(j,i)
4904 end subroutine friction_force
4905 !-----------------------------------------------------------------------------
4906 subroutine setup_fricmat
4910 use control_data, only:time_Bcast
4911 use control, only:tcpu
4913 ! implicit real*8 (a-h,o-z)
4917 real(kind=8) :: time00
4919 ! include 'DIMENSIONS'
4920 ! include 'COMMON.VAR'
4921 ! include 'COMMON.CHAIN'
4922 ! include 'COMMON.DERIV'
4923 ! include 'COMMON.GEO'
4924 ! include 'COMMON.LOCAL'
4925 ! include 'COMMON.INTERACT'
4926 ! include 'COMMON.MD'
4927 ! include 'COMMON.SETUP'
4928 ! include 'COMMON.TIME1'
4929 ! integer licznik /0/
4932 ! include 'COMMON.LANGEVIN'
4934 ! include 'COMMON.LANGEVIN.lang0'
4936 ! include 'COMMON.IOUNITS'
4938 integer :: i,j,ind,ind1,m
4939 logical :: lprn = .false.
4940 real(kind=8) :: dtdi !el ,gamvec(2*nres)
4941 !el real(kind=8),dimension(2*nres,2*nres) :: ginvfric,fcopy
4942 real(kind=8),dimension(2*nres,2*nres) :: fcopy
4943 !el real(kind=8),dimension(2*nres*(2*nres+1)/2) :: Ghalf !(mmaxres2) (mmaxres2=(maxres2*(maxres2+1)/2))
4944 !el common /syfek/ gamvec
4945 real(kind=8) :: work(8*2*nres)
4946 integer :: iwork(2*nres)
4947 !el common /przechowalnia/ ginvfric,Ghalf,fcopy
4948 integer :: ii,iti,k,l,nzero,nres2,nres6,ierr
4950 if (fg_rank.ne.king) goto 10
4955 if(.not.allocated(gamvec)) allocate(gamvec(nres2)) !(MAXRES2)
4956 if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
4957 !el if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
4958 !el allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
4959 if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !maxres2=2*maxres
4961 !el if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
4962 ! Zeroing out fricmat
4968 ! Load the friction coefficients corresponding to peptide groups
4974 ! Load the friction coefficients corresponding to side chains
4982 gamvec(ii)=gamsc(iabs(iti))
4984 if (surfarea) call sdarea(gamvec)
4986 ! write (iout,*) "Matrix A and vector gamma"
4988 ! write (iout,'(i2,$)') i
4990 ! write (iout,'(f4.1,$)') A(i,j)
4992 ! write (iout,'(f8.3)') gamvec(i)
4996 write (iout,*) "Vector gamvec"
4998 write (iout,'(i5,f10.5)') i, gamvec(i)
5002 ! The friction matrix
5007 dtdi=dtdi+A(j,k)*A(j,i)*gamvec(j)
5014 write (iout,'(//a)') "Matrix fricmat"
5015 call matout2(dimen,dimen,nres2,nres2,fricmat)
5017 if (lang.eq.2 .or. lang.eq.3) then
5018 ! Mass-scale the friction matrix if non-direct integration will be performed
5024 Ginvfric(i,j)=Ginvfric(i,j)+ &
5025 Gsqrm(i,k)*Gsqrm(l,j)*fricmat(k,l)
5030 ! Diagonalize the friction matrix
5035 Ghalf(ind)=Ginvfric(i,j)
5038 call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
5041 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
5042 " mass-scaled friction matrix"
5043 call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
5045 ! Precompute matrices for tinker stochastic integrator
5052 mt1(i,j)=mt1(i,j)+fricvec(k,i)*gsqrm(k,j)
5053 mt2(i,j)=mt2(i,j)+fricvec(k,i)*gsqrp(k,j)
5059 else if (lang.eq.4) then
5060 ! Diagonalize the friction matrix
5065 Ghalf(ind)=fricmat(i,j)
5068 call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
5071 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
5073 call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
5075 ! Determine the number of zero eigenvalues of the friction matrix
5076 nzero=max0(dimen-dimen1,0)
5077 ! do while (fricgam(nzero+1).le.1.0d-5 .and. nzero.lt.dimen)
5080 write (iout,*) "Number of zero eigenvalues:",nzero
5085 fricmat(i,j)=fricmat(i,j) &
5086 +fricvec(i,k)*fricvec(j,k)/fricgam(k)
5091 write (iout,'(//a)') "Generalized inverse of fricmat"
5092 call matout(dimen,dimen,nres6,nres6,fricmat)
5097 if (nfgtasks.gt.1) then
5098 if (fg_rank.eq.0) then
5099 ! The matching BROADCAST for fg processors is called in ERGASTULUM
5105 call MPI_Bcast(10,1,MPI_INTEGER,king,FG_COMM,IERROR)
5107 time_Bcast=time_Bcast+MPI_Wtime()-time00
5109 time_Bcast=time_Bcast+tcpu()-time00
5111 ! print *,"Processor",myrank,
5112 ! & " BROADCAST iorder in SETUP_FRICMAT"
5115 write (iout,*) "setup_fricmat licznik"!,licznik !sp
5121 ! Scatter the friction matrix
5122 call MPI_Scatterv(fricmat(1,1),nginv_counts(0),&
5123 nginv_start(0),MPI_DOUBLE_PRECISION,fcopy(1,1),&
5124 myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
5127 time_scatter=time_scatter+MPI_Wtime()-time00
5128 time_scatter_fmat=time_scatter_fmat+MPI_Wtime()-time00
5130 time_scatter=time_scatter+tcpu()-time00
5131 time_scatter_fmat=time_scatter_fmat+tcpu()-time00
5135 do j=1,2*my_ng_count
5136 fricmat(j,i)=fcopy(i,j)
5139 ! write (iout,*) "My chunk of fricmat"
5140 ! call MATOUT2(my_ng_count,dimen,maxres2,maxres2,fcopy)
5144 end subroutine setup_fricmat
5145 !-----------------------------------------------------------------------------
5146 subroutine stochastic_force(stochforcvec)
5149 use random, only:anorm_distr
5150 ! implicit real*8 (a-h,o-z)
5151 ! include 'DIMENSIONS'
5152 use control, only: tcpu
5157 ! include 'COMMON.VAR'
5158 ! include 'COMMON.CHAIN'
5159 ! include 'COMMON.DERIV'
5160 ! include 'COMMON.GEO'
5161 ! include 'COMMON.LOCAL'
5162 ! include 'COMMON.INTERACT'
5163 ! include 'COMMON.MD'
5164 ! include 'COMMON.TIME1'
5166 ! include 'COMMON.LANGEVIN'
5168 ! include 'COMMON.LANGEVIN.lang0'
5170 ! include 'COMMON.IOUNITS'
5172 real(kind=8) :: x,sig,lowb,highb
5173 real(kind=8) :: ff(3),force(3,0:2*nres),zeta2,lowb2
5174 real(kind=8) :: highb2,sig2,forcvec(6*nres),stochforcvec(6*nres)
5175 real(kind=8) :: time00
5176 logical :: lprn = .false.
5181 stochforc(j,i)=0.0d0
5191 ! Compute the stochastic forces acting on bodies. Store in force.
5197 force(j,i)=anorm_distr(x,sig,lowb,highb)
5205 force(j,i+nres)=anorm_distr(x,sig2,lowb2,highb2)
5209 time_fsample=time_fsample+MPI_Wtime()-time00
5211 time_fsample=time_fsample+tcpu()-time00
5213 ! Compute the stochastic forces acting on virtual-bond vectors.
5219 stochforc(j,i)=ff(j)+0.5d0*force(j,i)
5222 ff(j)=ff(j)+force(j,i)
5224 if (itype(i+1).ne.ntyp1) then
5226 stochforc(j,i)=stochforc(j,i)+force(j,i+nres+1)
5227 ff(j)=ff(j)+force(j,i+nres+1)
5232 stochforc(j,0)=ff(j)+force(j,nnt+nres)
5235 if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
5237 stochforc(j,i+nres)=force(j,i+nres)
5243 stochforcvec(j)=stochforc(j,0)
5248 stochforcvec(ind+j)=stochforc(j,i)
5253 if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
5255 stochforcvec(ind+j)=stochforc(j,i+nres)
5261 write (iout,*) "stochforcvec"
5263 write(iout,'(i5,e15.5)') i,stochforcvec(i)
5265 write(iout,*) "Stochastic forces backbone"
5267 write(iout,'(i5,3e15.5)') i,(stochforc(j,i),j=1,3)
5269 write(iout,*) "Stochastic forces side chain"
5271 write(iout,'(i5,3e15.5)') &
5272 i,(stochforc(j,i+nres),j=1,3)
5280 write (iout,*) i,ind
5282 forcvec(ind+j)=force(j,i)
5287 write (iout,*) i,ind
5289 forcvec(j+ind)=force(j,i+nres)
5294 write (iout,*) "forcvec"
5298 write (iout,'(2i3,2f10.5)') i,j,force(j,i),&
5305 write (iout,'(2i3,2f10.5)') i,j,force(j,i+nres),&
5314 end subroutine stochastic_force
5315 !-----------------------------------------------------------------------------
5316 subroutine sdarea(gamvec)
5318 ! Scale the friction coefficients according to solvent accessible surface areas
5319 ! Code adapted from TINKER
5323 ! implicit real*8 (a-h,o-z)
5324 ! include 'DIMENSIONS'
5325 ! include 'COMMON.CONTROL'
5326 ! include 'COMMON.VAR'
5327 ! include 'COMMON.MD'
5329 ! include 'COMMON.LANGEVIN'
5331 ! include 'COMMON.LANGEVIN.lang0'
5333 ! include 'COMMON.CHAIN'
5334 ! include 'COMMON.DERIV'
5335 ! include 'COMMON.GEO'
5336 ! include 'COMMON.LOCAL'
5337 ! include 'COMMON.INTERACT'
5338 ! include 'COMMON.IOUNITS'
5339 ! include 'COMMON.NAMES'
5340 real(kind=8),dimension(2*nres) :: radius,gamvec !(maxres2)
5341 real(kind=8),parameter :: twosix = 1.122462048309372981d0
5342 logical :: lprn = .false.
5343 real(kind=8) :: probe,area,ratio
5344 integer :: i,j,ind,iti
5346 ! determine new friction coefficients every few SD steps
5348 ! set the atomic radii to estimates of sigma values
5350 ! print *,"Entered sdarea"
5356 ! Load peptide group radii
5360 ! Load side chain radii
5363 radius(i+nres)=restok(iti)
5366 ! write (iout,*) "i",i," radius",radius(i)
5369 radius(i) = radius(i) / twosix
5370 if (radius(i) .ne. 0.0d0) radius(i) = radius(i) + probe
5373 ! scale atomic friction coefficients by accessible area
5375 if (lprn) write (iout,*) &
5376 "Original gammas, surface areas, scaling factors, new gammas, ",&
5377 "std's of stochastic forces"
5380 if (radius(i).gt.0.0d0) then
5381 call surfatom (i,area,radius)
5382 ratio = dmax1(area/(4.0d0*pi*radius(i)**2),1.0d-1)
5383 if (lprn) write (iout,'(i5,3f10.5,$)') &
5384 i,gamvec(ind+1),area,ratio
5387 gamvec(ind) = ratio * gamvec(ind)
5389 stdforcp(i)=stdfp*dsqrt(gamvec(ind))
5390 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcp(i)
5394 if (radius(i+nres).gt.0.0d0) then
5395 call surfatom (i+nres,area,radius)
5396 ratio = dmax1(area/(4.0d0*pi*radius(i+nres)**2),1.0d-1)
5397 if (lprn) write (iout,'(i5,3f10.5,$)') &
5398 i,gamvec(ind+1),area,ratio
5401 gamvec(ind) = ratio * gamvec(ind)
5403 stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamvec(ind))
5404 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcsc(i)
5409 end subroutine sdarea
5410 !-----------------------------------------------------------------------------
5412 !-----------------------------------------------------------------------------
5415 ! ###################################################
5416 ! ## COPYRIGHT (C) 1996 by Jay William Ponder ##
5417 ! ## All Rights Reserved ##
5418 ! ###################################################
5420 ! ################################################################
5422 ! ## subroutine surfatom -- exposed surface area of an atom ##
5424 ! ################################################################
5427 ! "surfatom" performs an analytical computation of the surface
5428 ! area of a specified atom; a simplified version of "surface"
5430 ! literature references:
5432 ! T. J. Richmond, "Solvent Accessible Surface Area and
5433 ! Excluded Volume in Proteins", Journal of Molecular Biology,
5436 ! L. Wesson and D. Eisenberg, "Atomic Solvation Parameters
5437 ! Applied to Molecular Dynamics of Proteins in Solution",
5438 ! Protein Science, 1, 227-235 (1992)
5440 ! variables and parameters:
5442 ! ir number of atom for which area is desired
5443 ! area accessible surface area of the atom
5444 ! radius radii of each of the individual atoms
5447 subroutine surfatom(ir,area,radius)
5449 ! implicit real*8 (a-h,o-z)
5450 ! include 'DIMENSIONS'
5452 ! include 'COMMON.GEO'
5453 ! include 'COMMON.IOUNITS'
5455 integer :: nsup,nstart_sup
5456 ! double precision c,dc,dc_old,d_c_work,xloc,xrot,dc_norm
5457 ! common /chain/ c(3,maxres2+2),dc(3,0:maxres2),dc_old(3,0:maxres2),
5458 ! & xloc(3,maxres),xrot(3,maxres),dc_norm(3,0:maxres2),
5459 ! & dc_work(MAXRES6),nres,nres0
5460 integer,parameter :: maxarc=300
5464 integer :: mi,ni,narc
5465 integer :: key(maxarc)
5466 integer :: intag(maxarc)
5467 integer :: intag1(maxarc)
5468 real(kind=8) :: area,arcsum
5469 real(kind=8) :: arclen,exang
5470 real(kind=8) :: delta,delta2
5471 real(kind=8) :: eps,rmove
5472 real(kind=8) :: xr,yr,zr
5473 real(kind=8) :: rr,rrsq
5474 real(kind=8) :: rplus,rminus
5475 real(kind=8) :: axx,axy,axz
5476 real(kind=8) :: ayx,ayy
5477 real(kind=8) :: azx,azy,azz
5478 real(kind=8) :: uxj,uyj,uzj
5479 real(kind=8) :: tx,ty,tz
5480 real(kind=8) :: txb,tyb,td
5481 real(kind=8) :: tr2,tr,txr,tyr
5482 real(kind=8) :: tk1,tk2
5483 real(kind=8) :: thec,the,t,tb
5484 real(kind=8) :: txk,tyk,tzk
5485 real(kind=8) :: t1,ti,tf,tt
5486 real(kind=8) :: txj,tyj,tzj
5487 real(kind=8) :: ccsq,cc,xysq
5488 real(kind=8) :: bsqk,bk,cosine
5489 real(kind=8) :: dsqj,gi,pix2
5490 real(kind=8) :: therk,dk,gk
5491 real(kind=8) :: risqk,rik
5492 real(kind=8) :: radius(2*nres) !(maxatm) (maxatm=maxres2)
5493 real(kind=8) :: ri(maxarc),risq(maxarc)
5494 real(kind=8) :: ux(maxarc),uy(maxarc),uz(maxarc)
5495 real(kind=8) :: xc(maxarc),yc(maxarc),zc(maxarc)
5496 real(kind=8) :: xc1(maxarc),yc1(maxarc),zc1(maxarc)
5497 real(kind=8) :: dsq(maxarc),bsq(maxarc)
5498 real(kind=8) :: dsq1(maxarc),bsq1(maxarc)
5499 real(kind=8) :: arci(maxarc),arcf(maxarc)
5500 real(kind=8) :: ex(maxarc),lt(maxarc),gr(maxarc)
5501 real(kind=8) :: b(maxarc),b1(maxarc),bg(maxarc)
5502 real(kind=8) :: kent(maxarc),kout(maxarc)
5503 real(kind=8) :: ther(maxarc)
5504 logical :: moved,top
5505 logical :: omit(maxarc)
5508 ! maxatm = 2*nres !maxres2 maxres2=2*maxres
5509 ! maxlight = 8*maxatm
5512 ! maxtors = 4*maxatm
5514 ! zero out the surface area for the sphere of interest
5517 ! write (2,*) "ir",ir," radius",radius(ir)
5518 if (radius(ir) .eq. 0.0d0) return
5520 ! set the overlap significance and connectivity shift
5524 delta2 = delta * delta
5529 ! store coordinates and radius of the sphere of interest
5537 ! initialize values of some counters and summations
5546 ! test each sphere to see if it overlaps the sphere of interest
5549 if (i.eq.ir .or. radius(i).eq.0.0d0) goto 30
5550 rplus = rr + radius(i)
5552 if (abs(tx) .ge. rplus) goto 30
5554 if (abs(ty) .ge. rplus) goto 30
5556 if (abs(tz) .ge. rplus) goto 30
5558 ! check for sphere overlap by testing distance against radii
5560 xysq = tx*tx + ty*ty
5561 if (xysq .lt. delta2) then
5568 if (rplus-cc .le. delta) goto 30
5569 rminus = rr - radius(i)
5571 ! check to see if sphere of interest is completely buried
5573 if (cc-abs(rminus) .le. delta) then
5574 if (rminus .le. 0.0d0) goto 170
5578 ! check for too many overlaps with sphere of interest
5580 if (io .ge. maxarc) then
5582 20 format (/,' SURFATOM -- Increase the Value of MAXARC')
5586 ! get overlap between current sphere and sphere of interest
5595 gr(io) = (ccsq+rplus*rminus) / (2.0d0*rr*b1(io))
5601 ! case where no other spheres overlap the sphere of interest
5604 area = 4.0d0 * pi * rrsq
5608 ! case where only one sphere overlaps the sphere of interest
5611 area = pix2 * (1.0d0 + gr(1))
5612 area = mod(area,4.0d0*pi) * rrsq
5616 ! case where many spheres intersect the sphere of interest;
5617 ! sort the intersecting spheres by their degree of overlap
5619 call sort2 (io,gr,key)
5622 intag(i) = intag1(k)
5631 ! get radius of each overlap circle on surface of the sphere
5636 risq(i) = rrsq - gi*gi
5637 ri(i) = sqrt(risq(i))
5638 ther(i) = 0.5d0*pi - asin(min(1.0d0,max(-1.0d0,gr(i))))
5641 ! find boundary of inaccessible area on sphere of interest
5644 if (.not. omit(k)) then
5651 ! check to see if J circle is intersecting K circle;
5652 ! get distance between circle centers and sum of radii
5655 if (omit(j)) goto 60
5656 cc = (txk*xc(j)+tyk*yc(j)+tzk*zc(j))/(bk*b(j))
5657 cc = acos(min(1.0d0,max(-1.0d0,cc)))
5658 td = therk + ther(j)
5660 ! check to see if circles enclose separate regions
5662 if (cc .ge. td) goto 60
5664 ! check for circle J completely inside circle K
5666 if (cc+ther(j) .lt. therk) goto 40
5668 ! check for circles that are essentially parallel
5670 if (cc .gt. delta) goto 50
5675 ! check to see if sphere of interest is completely buried
5678 if (pix2-cc .le. td) goto 170
5684 ! find T value of circle intersections
5687 if (omit(k)) goto 110
5702 ! rotation matrix elements
5714 if (.not. omit(j)) then
5719 ! rotate spheres so K vector colinear with z-axis
5721 uxj = txj*axx + tyj*axy - tzj*axz
5722 uyj = tyj*ayy - txj*ayx
5723 uzj = txj*azx + tyj*azy + tzj*azz
5724 cosine = min(1.0d0,max(-1.0d0,uzj/b(j)))
5725 if (acos(cosine) .lt. therk+ther(j)) then
5726 dsqj = uxj*uxj + uyj*uyj
5731 tr2 = risqk*dsqj - tb*tb
5737 ! get T values of intersection for K circle
5740 tb = min(1.0d0,max(-1.0d0,tb))
5742 if (tyb-txr .lt. 0.0d0) tk1 = pix2 - tk1
5744 tb = min(1.0d0,max(-1.0d0,tb))
5746 if (tyb+txr .lt. 0.0d0) tk2 = pix2 - tk2
5747 thec = (rrsq*uzj-gk*bg(j)) / (rik*ri(j)*b(j))
5748 if (abs(thec) .lt. 1.0d0) then
5750 else if (thec .ge. 1.0d0) then
5752 else if (thec .le. -1.0d0) then
5756 ! see if "tk1" is entry or exit point; check t=0 point;
5757 ! "ti" is exit point, "tf" is entry point
5759 cosine = min(1.0d0,max(-1.0d0, &
5760 (uzj*gk-uxj*rik)/(b(j)*rr)))
5761 if ((acos(cosine)-ther(j))*(tk2-tk1) .le. 0.0d0) then
5769 if (narc .ge. maxarc) then
5771 70 format (/,' SURFATOM -- Increase the Value',&
5775 if (tf .le. ti) then
5796 ! special case; K circle without intersections
5798 if (narc .le. 0) goto 90
5800 ! general case; sum up arclength and set connectivity code
5802 call sort2 (narc,arci,key)
5807 if (narc .gt. 1) then
5810 if (t .lt. arci(j)) then
5811 arcsum = arcsum + arci(j) - t
5812 exang = exang + ex(ni)
5814 if (jb .ge. maxarc) then
5816 80 format (/,' SURFATOM -- Increase the Value',&
5821 kent(jb) = maxarc*i + k
5823 kout(jb) = maxarc*k + i
5832 arcsum = arcsum + pix2 - t
5834 exang = exang + ex(ni)
5837 kent(jb) = maxarc*i + k
5839 kout(jb) = maxarc*k + i
5846 arclen = arclen + gr(k)*arcsum
5849 if (arclen .eq. 0.0d0) goto 170
5850 if (jb .eq. 0) goto 150
5852 ! find number of independent boundaries and check connectivity
5856 if (kout(k) .ne. 0) then
5863 if (m .eq. kent(ii)) then
5866 if (j .eq. jb) goto 150
5878 ! attempt to fix connectivity error by moving atom slightly
5882 140 format (/,' SURFATOM -- Connectivity Error at Atom',i6)
5891 ! compute the exposed surface area for the sphere of interest
5894 area = ib*pix2 + exang + arclen
5895 area = mod(area,4.0d0*pi) * rrsq
5897 ! attempt to fix negative area by moving atom slightly
5899 if (area .lt. 0.0d0) then
5902 160 format (/,' SURFATOM -- Negative Area at Atom',i6)
5913 end subroutine surfatom
5914 !----------------------------------------------------------------
5915 !----------------------------------------------------------------
5916 subroutine alloc_MD_arrays
5917 !EL Allocation of arrays used by MD module
5919 integer :: nres2,nres6
5922 !----------------------
5926 allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
5927 allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
5928 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
5929 allocate(fricvec(nres2,nres2))
5930 allocate(pfric_mat(nres2,nres2),vfric_mat(nres2,nres2))
5931 allocate(afric_mat(nres2,nres2),prand_mat(nres2,nres2))
5932 allocate(vrand_mat1(nres2,nres2),vrand_mat2(nres2,nres2)) !(MAXRES2,MAXRES2)
5933 allocate(pfric0_mat(nres2,nres2,0:maxflag_stoch))
5934 allocate(afric0_mat(nres2,nres2,0:maxflag_stoch))
5935 allocate(vfric0_mat(nres2,nres2,0:maxflag_stoch))
5936 allocate(prand0_mat(nres2,nres2,0:maxflag_stoch))
5937 allocate(vrand0_mat1(nres2,nres2,0:maxflag_stoch))
5938 allocate(vrand0_mat2(nres2,nres2,0:maxflag_stoch)) !(MAXRES2,MAXRES2,0:maxflag_stoch)
5939 allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
5941 allocate(mt1(nres2,nres2),mt2(nres2,nres2),mt3(nres2,nres2)) !(maxres2,maxres2)
5942 !----------------------
5944 ! commom.langevin.lang0
5946 allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
5947 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
5948 allocate(fricvec(nres2,nres2)) !(MAXRES2,MAXRES2)
5949 allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
5950 allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
5953 !el if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2))
5954 !----------------------
5955 ! commom.hairpin in CSA module
5956 !----------------------
5957 ! common.mce in MCM_MD module
5958 !----------------------
5960 ! common /mdgrad/ in module.energy
5961 ! common /back_constr/ in module.energy
5962 ! common /qmeas/ in module.energy
5965 allocate(potEcomp(0:n_ene+4)) !(0:n_ene+4)
5967 allocate(d_t(3,0:nres2),d_a(3,0:nres2),d_t_old(3,0:nres2)) !(3,0:MAXRES2)
5968 allocate(d_a_work(nres6)) !(6*MAXRES)
5970 allocate(DM(nres2),DU1(nres2),DU2(nres2))
5971 allocate(DMorig(nres2),DU1orig(nres2),DU2orig(nres2))
5973 allocate(Gmat(nres2,nres2),A(nres2,nres2))
5974 if(.not.allocated(Ginv)) allocate(Ginv(nres2,nres2)) !in control: ergastulum
5975 allocate(Gsqrp(nres2,nres2),Gsqrm(nres2,nres2),Gvec(nres2,nres2)) !(maxres2,maxres2)
5977 allocate(Geigen(nres2)) !(maxres2)
5978 if(.not.allocated(vtot)) allocate(vtot(nres2)) !(maxres2)
5979 ! common /inertia/ in io_conf: parmread
5980 ! real(kind=8),dimension(:),allocatable :: ISC,msc !(ntyp+1)
5981 ! common /langevin/in io read_MDpar
5982 ! real(kind=8),dimension(:),allocatable :: gamsc !(ntyp1)
5983 ! real(kind=8),dimension(:),allocatable :: stdfsc !(ntyp)
5984 ! in io_conf: parmread
5985 ! real(kind=8),dimension(:),allocatable :: restok !(ntyp+1)
5986 ! common /mdpmpi/ in control: ergastulum
5987 if(.not.allocated(ng_start)) allocate(ng_start(0:nfgtasks-1))
5988 if(.not.allocated(ng_counts)) allocate(ng_counts(0:nfgtasks-1))
5989 if(.not.allocated(nginv_counts)) allocate(nginv_counts(0:nfgtasks-1)) !(0:MaxProcs-1)
5990 if(.not.allocated(nginv_start)) allocate(nginv_start(0:nfgtasks)) !(0:MaxProcs)
5991 !----------------------
5992 ! common.muca in read_muca
5993 ! common /double_muca/
5994 ! real(kind=8) :: elow,ehigh,factor,hbin,factor_min
5995 ! real(kind=8),dimension(:),allocatable :: emuca,nemuca,&
5996 ! nemuca2,hist !(4*maxres)
5997 ! common /integer_muca/
5998 ! integer :: nmuca,imtime,muca_smooth
6000 ! real(kind=8),dimension(:),allocatable :: elowi,ehighi !(maxprocs)
6001 !----------------------
6003 ! common /mdgrad/ in module.energy
6004 ! common /back_constr/ in module.energy
6005 ! common /qmeas/ in module.energy
6009 allocate(d_t_work(nres6),d_t_work_new(nres6),d_af_work(nres6))
6010 allocate(d_as_work(nres6),kinetic_force(nres6)) !(MAXRES6)
6011 allocate(d_t_new(3,0:nres2),d_a_old(3,0:nres2),d_a_short(3,0:nres2)) !,d_a !(3,0:MAXRES2)
6012 allocate(stdforcp(nres),stdforcsc(nres)) !(MAXRES)
6013 !----------------------
6015 allocate(D_ban(nres6)) !(MAXRES6) maxres6=6*maxres
6016 ! common /stochcalc/ stochforcvec
6017 allocate(stochforcvec(nres6)) !(MAXRES6) maxres6=6*maxres
6018 !----------------------
6020 end subroutine alloc_MD_arrays
6021 !-----------------------------------------------------------------------------
6022 !-----------------------------------------------------------------------------