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,1).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,1).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,1).ne.10) then
311 Td(i)=Td(i)+vbldsc0(1,itype(k,1))*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,1).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,1).ne.10) then
415 xx=vbld(i+nres)-vbldsc0(1,itype(i,1))
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,1).ne.10) then
443 blen2 = scalar(dc(1,i+nres),dc(1,i+nres))
444 ppvec(ind)=2*vbldsc0(1,itype(i,1))**2-blen2
445 diffbond=dabs(vbldsc0(1,itype(i,1))-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,1).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,1).ne.10) then
530 xx=vbld(i+nres)-vbldsc0(1,itype(i,1))
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)
552 ! Calculate the kinetic and total energy and the kinetic temperature
555 t_enegrad=t_enegrad+MPI_Wtime()-tt0
557 t_enegrad=t_enegrad+tcpu()-tt0
560 kinetic_T=2.0d0/(dimen*Rb)*EK
562 end subroutine brown_step
563 !-----------------------------------------------------------------------------
565 !-----------------------------------------------------------------------------
566 subroutine gauss(RO,AP,MT,M,N,*)
568 ! CALCULATES (RO**(-1))*AP BY GAUSS ELIMINATION
569 ! RO IS A SQUARE MATRIX
570 ! THE CALCULATED PRODUCT IS STORED IN AP
571 ! ABNORMAL EXIT IF RO IS SINGULAR
573 integer :: MT, M, N, M1,I,J,IM,&
575 real(kind=8) :: RO(MT,M),AP(MT,N),X,RM,PR,Y
581 if(dabs(X).le.1.0D-13) return 1
593 if(DABS(RO(J,I)).LE.RM) goto 2
607 if(dabs(X).le.1.0E-13) return 1
616 8 AP(J,K)=AP(J,K)-Y*AP(I,K)
618 9 RO(J,K)=RO(J,K)-Y*RO(I,K)
622 if(dabs(X).le.1.0E-13) return 1
632 15 X=X-AP(K,J)*RO(MI,K)
637 !-----------------------------------------------------------------------------
639 !-----------------------------------------------------------------------------
640 subroutine kinetic(KE_total)
641 !----------------------------------------------------------------
642 ! This subroutine calculates the total kinetic energy of the chain
643 !-----------------------------------------------------------------
645 ! implicit real*8 (a-h,o-z)
646 ! include 'DIMENSIONS'
647 ! include 'COMMON.VAR'
648 ! include 'COMMON.CHAIN'
649 ! include 'COMMON.DERIV'
650 ! include 'COMMON.GEO'
651 ! include 'COMMON.LOCAL'
652 ! include 'COMMON.INTERACT'
653 ! include 'COMMON.MD'
654 ! include 'COMMON.IOUNITS'
655 real(kind=8) :: KE_total
658 real(kind=8) :: KEt_p,KEt_sc,KEr_p,KEr_sc,incr(3),&
661 write (iout,*) "Velocities, kietic"
663 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
664 (d_t(j,i+nres),j=1,3)
669 ! write (iout,*) "ISC",(isc(itype(i,1)),i=1,nres)
670 ! The translational part for peptide virtual bonds
675 ! write (iout,*) "Kinetic trp:",i,(incr(j),j=1,3)
677 v(j)=incr(j)+0.5d0*d_t(j,i)
679 vtot(i)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
680 KEt_p=KEt_p+(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
682 incr(j)=incr(j)+d_t(j,i)
685 ! write(iout,*) 'KEt_p', KEt_p
686 ! The translational part for the side chain virtual bond
687 ! Only now we can initialize incr with zeros. It must be equal
688 ! to the velocities of the first Calpha.
694 if (itype(i,1).eq.10) then
700 v(j)=incr(j)+d_t(j,nres+i)
703 ! write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3)
704 ! write (iout,*) "i",i," msc",msc(iti)," v",(v(j),j=1,3)
705 KEt_sc=KEt_sc+msc(iti,1)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
706 vtot(i+nres)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
708 incr(j)=incr(j)+d_t(j,i)
712 ! write(iout,*) 'KEt_sc', KEt_sc
713 ! The part due to stretching and rotation of the peptide groups
716 ! write (iout,*) "i",i
717 ! write (iout,*) "i",i," mag1",mag1," mag2",mag2
721 ! write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3)
722 KEr_p=KEr_p+(incr(1)*incr(1)+incr(2)*incr(2) &
726 ! write(iout,*) 'KEr_p', KEr_p
727 ! The rotational part of the side chain virtual bond
731 if (itype(i,1).ne.10) then
733 incr(j)=d_t(j,nres+i)
735 ! write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3)
736 KEr_sc=KEr_sc+Isc(iti,1)*(incr(1)*incr(1)+incr(2)*incr(2)+ &
740 ! The total kinetic energy
742 ! write(iout,*) 'KEr_sc', KEr_sc
743 KE_total=0.5d0*(mp(1)*KEt_p+KEt_sc+0.25d0*Ip(1)*KEr_p+KEr_sc)
744 ! write (iout,*) "KE_total",KE_total
746 end subroutine kinetic
747 !-----------------------------------------------------------------------------
749 !-----------------------------------------------------------------------------
751 !------------------------------------------------
752 ! The driver for molecular dynamics subroutines
753 !------------------------------------------------
756 use control, only:tcpu,ovrtim
757 ! use io_comm, only:ilen
759 use compare, only:secondary2,hairpin
760 use io, only:cartout,statout
761 ! implicit real*8 (a-h,o-z)
762 ! include 'DIMENSIONS'
765 integer :: IERROR,ERRCODE
767 ! include 'COMMON.SETUP'
768 ! include 'COMMON.CONTROL'
769 ! include 'COMMON.VAR'
770 ! include 'COMMON.MD'
772 ! include 'COMMON.LANGEVIN'
774 ! include 'COMMON.LANGEVIN.lang0'
776 ! include 'COMMON.CHAIN'
777 ! include 'COMMON.DERIV'
778 ! include 'COMMON.GEO'
779 ! include 'COMMON.LOCAL'
780 ! include 'COMMON.INTERACT'
781 ! include 'COMMON.IOUNITS'
782 ! include 'COMMON.NAMES'
783 ! include 'COMMON.TIME1'
784 ! include 'COMMON.HAIRPIN'
785 real(kind=8),dimension(3) :: L,vcm
787 real(kind=8),dimension(6*nres) :: v_work,v_transf !(maxres6) maxres6=6*maxres
789 integer :: rstcount !ilen,
791 character(len=50) :: tytul
792 !el common /gucio/ cm
793 integer :: itime,i,j,nharp
794 integer,dimension(4,nres/3) :: iharp !(4,nres/3)(4,maxres/3)
796 real(kind=8) :: tt0,scalfac
801 if (ilen(tmpdir).gt.0) &
802 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_" &
803 //liczba(:ilen(liczba))//'.rst')
805 if (ilen(tmpdir).gt.0) &
806 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"//'.rst')
813 write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
819 ! Determine the inverse of the inertia matrix.
820 call setup_MD_matrices
824 t_MDsetup = MPI_Wtime()-tt0
826 t_MDsetup = tcpu()-tt0
829 ! Entering the MD loop
835 if (lang.eq.2 .or. lang.eq.3) then
839 call sd_verlet_p_setup
841 call sd_verlet_ciccotti_setup
845 pfric0_mat(i,j,0)=pfric_mat(i,j)
846 afric0_mat(i,j,0)=afric_mat(i,j)
847 vfric0_mat(i,j,0)=vfric_mat(i,j)
848 prand0_mat(i,j,0)=prand_mat(i,j)
849 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
850 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
855 flag_stoch(i)=.false.
859 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
861 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
865 else if (lang.eq.1 .or. lang.eq.4) then
869 t_langsetup=MPI_Wtime()-tt0
872 t_langsetup=tcpu()-tt0
875 do itime=1,n_timestep
877 if (large.and. mod(itime,ntwe).eq.0) &
878 write (iout,*) "itime",itime
880 if (lang.gt.0 .and. surfarea .and. &
881 mod(itime,reset_fricmat).eq.0) then
882 if (lang.eq.2 .or. lang.eq.3) then
886 call sd_verlet_p_setup
888 call sd_verlet_ciccotti_setup
892 pfric0_mat(i,j,0)=pfric_mat(i,j)
893 afric0_mat(i,j,0)=afric_mat(i,j)
894 vfric0_mat(i,j,0)=vfric_mat(i,j)
895 prand0_mat(i,j,0)=prand_mat(i,j)
896 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
897 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
902 flag_stoch(i)=.false.
905 else if (lang.eq.1 .or. lang.eq.4) then
908 write (iout,'(a,i10)') &
909 "Friction matrix reset based on surface area, itime",itime
911 if (reset_vel .and. tbf .and. lang.eq.0 &
912 .and. mod(itime,count_reset_vel).eq.0) then
914 write(iout,'(a,f20.2)') &
915 "Velocities reset to random values, time",totT
918 d_t_old(j,i)=d_t(j,i)
922 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
926 d_t(j,0)=d_t(j,0)-vcm(j)
929 kinetic_T=2.0d0/(dimen3*Rb)*EK
930 scalfac=dsqrt(T_bath/kinetic_T)
931 write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
934 d_t_old(j,i)=scalfac*d_t(j,i)
940 ! Time-reversible RESPA algorithm
941 ! (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
942 call RESPA_step(itime)
944 ! Variable time step algorithm.
945 call velverlet_step(itime)
949 call brown_step(itime)
951 print *,"Brown dynamics not here!"
953 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
959 if (mod(itime,ntwe).eq.0) call statout(itime)
972 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
975 v_work(ind)=d_t(j,i+nres)
980 write (66,'(80f10.5)') &
981 ((d_t(j,i),j=1,3),i=0,nres-1),((d_t(j,i+nres),j=1,3),i=1,nres)
985 v_transf(i)=v_transf(i)+gvec(j,i)*v_work(j)
987 v_transf(i)= v_transf(i)*dsqrt(geigen(i))
989 write (67,'(80f10.5)') (v_transf(i),i=1,ind)
992 if (mod(itime,ntwx).eq.0) then
993 write (tytul,'("time",f8.2)') totT
995 call hairpin(.true.,nharp,iharp)
996 call secondary2(.true.)
997 call pdbout(potE,tytul,ipdb)
1002 if (rstcount.eq.1000.or.itime.eq.n_timestep) then
1003 open(irest2,file=rest2name,status='unknown')
1004 write(irest2,*) totT,EK,potE,totE,t_bath
1006 ! AL 4/17/17: Now writing d_t(0,:) too
1008 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
1010 ! AL 4/17/17: Now writing d_c(0,:) too
1012 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
1020 t_MD=MPI_Wtime()-tt0
1024 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))') &
1026 'MD calculations setup:',t_MDsetup,&
1027 'Energy & gradient evaluation:',t_enegrad,&
1028 'Stochastic MD setup:',t_langsetup,&
1029 'Stochastic MD step setup:',t_sdsetup,&
1031 write (iout,'(/28(1h=),a25,27(1h=))') &
1032 ' End of MD calculation '
1034 write (iout,*) "time for etotal",t_etotal," elong",t_elong,&
1036 write (iout,*) "time_fric",time_fric," time_stoch",time_stoch,&
1037 " time_fricmatmult",time_fricmatmult," time_fsample ",&
1042 !-----------------------------------------------------------------------------
1043 subroutine velverlet_step(itime)
1044 !-------------------------------------------------------------------------------
1045 ! Perform a single velocity Verlet step; the time step can be rescaled if
1046 ! increments in accelerations exceed the threshold
1047 !-------------------------------------------------------------------------------
1048 ! implicit real*8 (a-h,o-z)
1049 ! include 'DIMENSIONS'
1051 use control, only:tcpu
1055 integer :: ierror,ierrcode
1056 real(kind=8) :: errcode
1058 ! include 'COMMON.SETUP'
1059 ! include 'COMMON.VAR'
1060 ! include 'COMMON.MD'
1062 ! include 'COMMON.LANGEVIN'
1064 ! include 'COMMON.LANGEVIN.lang0'
1066 ! include 'COMMON.CHAIN'
1067 ! include 'COMMON.DERIV'
1068 ! include 'COMMON.GEO'
1069 ! include 'COMMON.LOCAL'
1070 ! include 'COMMON.INTERACT'
1071 ! include 'COMMON.IOUNITS'
1072 ! include 'COMMON.NAMES'
1073 ! include 'COMMON.TIME1'
1074 ! include 'COMMON.MUCA'
1075 real(kind=8),dimension(3) :: vcm,incr
1076 real(kind=8),dimension(3) :: L
1077 integer :: count,rstcount !ilen,
1079 character(len=50) :: tytul
1080 integer :: maxcount_scale = 20
1081 !el common /gucio/ cm
1082 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1083 !el common /stochcalc/ stochforcvec
1084 integer :: itime,icount_scale,itime_scal,i,j,ifac_time
1086 real(kind=8) :: epdrift,tt0,fac_time
1088 if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres)) !(MAXRES6) maxres6=6*maxres
1094 else if (lang.eq.2 .or. lang.eq.3) then
1096 call stochastic_force(stochforcvec)
1099 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
1101 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1108 icount_scale=icount_scale+1
1109 if (icount_scale.gt.maxcount_scale) then
1111 "ERROR: too many attempts at scaling down the time step. ",&
1112 "amax=",amax,"epdrift=",epdrift,&
1113 "damax=",damax,"edriftmax=",edriftmax,&
1117 call MPI_Abort(MPI_COMM_WORLD,IERROR,IERRCODE)
1121 ! First step of the velocity Verlet algorithm
1126 else if (lang.eq.3) then
1128 call sd_verlet1_ciccotti
1130 else if (lang.eq.1) then
1135 ! Build the chain from the newly calculated coordinates
1136 call chainbuild_cart
1137 if (rattle) call rattle1
1139 if (large.and. mod(itime,ntwe).eq.0) then
1140 write (iout,*) "Cartesian and internal coordinates: step 1"
1145 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),&
1146 (dc(j,i+nres),j=1,3)
1148 write (iout,*) "Accelerations"
1150 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1151 (d_a(j,i+nres),j=1,3)
1153 write (iout,*) "Velocities, step 1"
1155 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1156 (d_t(j,i+nres),j=1,3)
1165 ! Calculate energy and forces
1167 call etotal(potEcomp)
1168 ! AL 4/17/17: Reduce the steps if NaNs occurred.
1169 if (potEcomp(0).gt.0.99e20 .or. isnan(potEcomp(0)).gt.0) then
1174 if (large.and. mod(itime,ntwe).eq.0) &
1175 call enerprint(potEcomp)
1178 t_etotal=t_etotal+MPI_Wtime()-tt0
1180 t_etotal=t_etotal+tcpu()-tt0
1183 potE=potEcomp(0)-potEcomp(20)
1185 ! Get the new accelerations
1188 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1190 t_enegrad=t_enegrad+tcpu()-tt0
1192 ! Determine maximum acceleration and scale down the timestep if needed
1194 amax=amax/(itime_scal+1)**2
1195 call predict_edrift(epdrift)
1196 if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then
1197 ! Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
1199 ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax)) &
1201 itime_scal=itime_scal+ifac_time
1202 ! fac_time=dmin1(damax/amax,0.5d0)
1203 fac_time=0.5d0**ifac_time
1204 d_time=d_time*fac_time
1205 if (lang.eq.2 .or. lang.eq.3) then
1207 ! write (iout,*) "Calling sd_verlet_setup: 1"
1208 ! Rescale the stochastic forces and recalculate or restore
1209 ! the matrices of tinker integrator
1210 if (itime_scal.gt.maxflag_stoch) then
1211 if (large) write (iout,'(a,i5,a)') &
1212 "Calculate matrices for stochastic step;",&
1213 " itime_scal ",itime_scal
1215 call sd_verlet_p_setup
1217 call sd_verlet_ciccotti_setup
1219 write (iout,'(2a,i3,a,i3,1h.)') &
1220 "Warning: cannot store matrices for stochastic",&
1221 " integration because the index",itime_scal,&
1222 " is greater than",maxflag_stoch
1223 write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",&
1224 " integration Langevin algorithm for better efficiency."
1225 else if (flag_stoch(itime_scal)) then
1226 if (large) write (iout,'(a,i5,a,l1)') &
1227 "Restore matrices for stochastic step; itime_scal ",&
1228 itime_scal," flag ",flag_stoch(itime_scal)
1231 pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
1232 afric_mat(i,j)=afric0_mat(i,j,itime_scal)
1233 vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
1234 prand_mat(i,j)=prand0_mat(i,j,itime_scal)
1235 vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
1236 vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
1240 if (large) write (iout,'(2a,i5,a,l1)') &
1241 "Calculate & store matrices for stochastic step;",&
1242 " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
1244 call sd_verlet_p_setup
1246 call sd_verlet_ciccotti_setup
1248 flag_stoch(ifac_time)=.true.
1251 pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
1252 afric0_mat(i,j,itime_scal)=afric_mat(i,j)
1253 vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
1254 prand0_mat(i,j,itime_scal)=prand_mat(i,j)
1255 vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
1256 vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
1260 fac_time=1.0d0/dsqrt(fac_time)
1262 stochforcvec(i)=fac_time*stochforcvec(i)
1265 else if (lang.eq.1) then
1266 ! Rescale the accelerations due to stochastic forces
1267 fac_time=1.0d0/dsqrt(fac_time)
1269 d_as_work(i)=d_as_work(i)*fac_time
1272 if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)') &
1273 "itime",itime," Timestep scaled down to ",&
1274 d_time," ifac_time",ifac_time," itime_scal",itime_scal
1276 ! Second step of the velocity Verlet algorithm
1281 else if (lang.eq.3) then
1283 call sd_verlet2_ciccotti
1285 else if (lang.eq.1) then
1290 if (rattle) call rattle2
1293 if (d_time.ne.d_time0) then
1296 if (lang.eq.2 .or. lang.eq.3) then
1297 if (large) write (iout,'(a)') &
1298 "Restore original matrices for stochastic step"
1299 ! write (iout,*) "Calling sd_verlet_setup: 2"
1300 ! Restore the matrices of tinker integrator if the time step has been restored
1303 pfric_mat(i,j)=pfric0_mat(i,j,0)
1304 afric_mat(i,j)=afric0_mat(i,j,0)
1305 vfric_mat(i,j)=vfric0_mat(i,j,0)
1306 prand_mat(i,j)=prand0_mat(i,j,0)
1307 vrand_mat1(i,j)=vrand0_mat1(i,j,0)
1308 vrand_mat2(i,j)=vrand0_mat2(i,j,0)
1317 ! Calculate the kinetic and the total energy and the kinetic temperature
1321 ! call kinetic1(EK1)
1322 ! write (iout,*) "step",itime," EK",EK," EK1",EK1
1324 ! Couple the system to Berendsen bath if needed
1325 if (tbf .and. lang.eq.0) then
1328 kinetic_T=2.0d0/(dimen3*Rb)*EK
1329 ! Backup the coordinates, velocities, and accelerations
1333 d_t_old(j,i)=d_t(j,i)
1334 d_a_old(j,i)=d_a(j,i)
1338 if (mod(itime,ntwe).eq.0 .and. large) then
1339 write (iout,*) "Velocities, step 2"
1341 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1342 (d_t(j,i+nres),j=1,3)
1347 end subroutine velverlet_step
1348 !-----------------------------------------------------------------------------
1349 subroutine RESPA_step(itime)
1350 !-------------------------------------------------------------------------------
1351 ! Perform a single RESPA step.
1352 !-------------------------------------------------------------------------------
1353 ! implicit real*8 (a-h,o-z)
1354 ! include 'DIMENSIONS'
1358 use control, only:tcpu
1360 ! use io_conf, only:cartprint
1363 integer :: IERROR,ERRCODE
1365 ! include 'COMMON.SETUP'
1366 ! include 'COMMON.CONTROL'
1367 ! include 'COMMON.VAR'
1368 ! include 'COMMON.MD'
1370 ! include 'COMMON.LANGEVIN'
1372 ! include 'COMMON.LANGEVIN.lang0'
1374 ! include 'COMMON.CHAIN'
1375 ! include 'COMMON.DERIV'
1376 ! include 'COMMON.GEO'
1377 ! include 'COMMON.LOCAL'
1378 ! include 'COMMON.INTERACT'
1379 ! include 'COMMON.IOUNITS'
1380 ! include 'COMMON.NAMES'
1381 ! include 'COMMON.TIME1'
1382 real(kind=8),dimension(0:n_ene) :: energia_short,energia_long
1383 real(kind=8),dimension(3) :: L,vcm,incr
1384 real(kind=8),dimension(3,0:2*nres) :: dc_old0,d_t_old0,d_a_old0 !(3,0:maxres2) maxres2=2*maxres
1385 logical :: PRINT_AMTS_MSG = .false.
1386 integer :: count,rstcount !ilen,
1388 character(len=50) :: tytul
1389 integer :: maxcount_scale = 10
1390 !el common /gucio/ cm
1391 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1392 !el common /stochcalc/ stochforcvec
1393 integer :: itime,itt,i,j,itsplit
1395 !el common /cipiszcze/ itt
1397 real(kind=8) :: epdrift,tt0,epdriftmax
1400 if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres)) !(MAXRES6) maxres6=6*maxres
1404 if (large.and. mod(itime,ntwe).eq.0) then
1405 write (iout,*) "***************** RESPA itime",itime
1406 write (iout,*) "Cartesian and internal coordinates: step 0"
1408 call pdbout(0.0d0,"cipiszcze",iout)
1410 write (iout,*) "Accelerations from long-range forces"
1412 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1413 (d_a(j,i+nres),j=1,3)
1415 write (iout,*) "Velocities, step 0"
1417 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1418 (d_t(j,i+nres),j=1,3)
1423 ! Perform the initial RESPA step (increment velocities)
1424 ! write (iout,*) "*********************** RESPA ini"
1427 if (mod(itime,ntwe).eq.0 .and. large) then
1428 write (iout,*) "Velocities, end"
1430 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1431 (d_t(j,i+nres),j=1,3)
1435 ! Compute the short-range forces
1441 ! 7/2/2009 commented out
1443 ! call etotal_short(energia_short)
1446 ! 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
1449 d_a(j,i)=d_a_short(j,i)
1453 if (large.and. mod(itime,ntwe).eq.0) then
1454 write (iout,*) "energia_short",energia_short(0)
1455 write (iout,*) "Accelerations from short-range forces"
1457 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1458 (d_a(j,i+nres),j=1,3)
1463 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1465 t_enegrad=t_enegrad+tcpu()-tt0
1470 d_t_old(j,i)=d_t(j,i)
1471 d_a_old(j,i)=d_a(j,i)
1474 ! 6/30/08 A-MTS: attempt at increasing the split number
1477 dc_old0(j,i)=dc_old(j,i)
1478 d_t_old0(j,i)=d_t_old(j,i)
1479 d_a_old0(j,i)=d_a_old(j,i)
1482 if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
1483 if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
1490 ! write (iout,*) "itime",itime," ntime_split",ntime_split
1491 ! Split the time step
1492 d_time=d_time0/ntime_split
1493 ! Perform the short-range RESPA steps (velocity Verlet increments of
1494 ! positions and velocities using short-range forces)
1495 ! write (iout,*) "*********************** RESPA split"
1496 do itsplit=1,ntime_split
1499 else if (lang.eq.2 .or. lang.eq.3) then
1501 call stochastic_force(stochforcvec)
1504 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
1506 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1511 ! First step of the velocity Verlet algorithm
1516 else if (lang.eq.3) then
1518 call sd_verlet1_ciccotti
1520 else if (lang.eq.1) then
1525 ! Build the chain from the newly calculated coordinates
1526 call chainbuild_cart
1527 if (rattle) call rattle1
1529 if (large.and. mod(itime,ntwe).eq.0) then
1530 write (iout,*) "***** ITSPLIT",itsplit
1531 write (iout,*) "Cartesian and internal coordinates: step 1"
1532 call pdbout(0.0d0,"cipiszcze",iout)
1535 write (iout,*) "Velocities, step 1"
1537 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1538 (d_t(j,i+nres),j=1,3)
1547 ! Calculate energy and forces
1549 call etotal_short(energia_short)
1550 ! AL 4/17/17: Exit itime_split loop when energy goes infinite
1551 if (energia_short(0).gt.0.99e20 .or. isnan(energia_short(0)) ) then
1552 if (PRINT_AMTS_MSG) &
1553 write (iout,*) "Infinities/NaNs in energia_short",energia_short(0),"; increasing ntime_split to",ntime_split
1554 ntime_split=ntime_split*2
1555 if (ntime_split.gt.maxtime_split) then
1558 "Cannot rescue the run; aborting job. Retry with a smaller time step"
1560 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1563 "Cannot rescue the run; terminating. Retry with a smaller time step"
1569 if (large.and. mod(itime,ntwe).eq.0) &
1570 call enerprint(energia_short)
1573 t_eshort=t_eshort+MPI_Wtime()-tt0
1575 t_eshort=t_eshort+tcpu()-tt0
1579 ! Get the new accelerations
1581 ! 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
1584 d_a_short(j,i)=d_a(j,i)
1588 if (large.and. mod(itime,ntwe).eq.0) then
1589 write (iout,*)"energia_short",energia_short(0)
1590 write (iout,*) "Accelerations from short-range forces"
1592 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1593 (d_a(j,i+nres),j=1,3)
1598 ! Determine maximum acceleration and scale down the timestep if needed
1600 amax=amax/ntime_split**2
1601 call predict_edrift(epdrift)
1602 if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0) &
1603 write (iout,*) "amax",amax," damax",damax,&
1604 " epdrift",epdrift," epdriftmax",epdriftmax
1605 ! Exit loop and try with increased split number if the change of
1606 ! acceleration is too big
1607 if (amax.gt.damax .or. epdrift.gt.edriftmax) then
1608 if (ntime_split.lt.maxtime_split) then
1610 ntime_split=ntime_split*2
1611 ! AL 4/17/17: We should exit the itime_split loop when acceleration change is too big
1615 dc_old(j,i)=dc_old0(j,i)
1616 d_t_old(j,i)=d_t_old0(j,i)
1617 d_a_old(j,i)=d_a_old0(j,i)
1620 if (PRINT_AMTS_MSG) then
1621 write (iout,*) "acceleration/energy drift too large",amax,&
1622 epdrift," split increased to ",ntime_split," itime",itime,&
1628 "Uh-hu. Bumpy landscape. Maximum splitting number",&
1630 " already reached!!! Trying to carry on!"
1634 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1636 t_enegrad=t_enegrad+tcpu()-tt0
1638 ! Second step of the velocity Verlet algorithm
1643 else if (lang.eq.3) then
1645 call sd_verlet2_ciccotti
1647 else if (lang.eq.1) then
1652 if (rattle) call rattle2
1653 ! Backup the coordinates, velocities, and accelerations
1657 d_t_old(j,i)=d_t(j,i)
1658 d_a_old(j,i)=d_a(j,i)
1665 ! Restore the time step
1667 ! Compute long-range forces
1674 call etotal_long(energia_long)
1675 if (energia_long(0).gt.0.99e20 .or. isnan(energia_long(0))) then
1678 "Infinitied/NaNs in energia_long, Aborting MPI job."
1680 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1682 write (iout,*) "Infinitied/NaNs in energia_long, terminating."
1686 if (large.and. mod(itime,ntwe).eq.0) &
1687 call enerprint(energia_long)
1690 t_elong=t_elong+MPI_Wtime()-tt0
1692 t_elong=t_elong+tcpu()-tt0
1698 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1700 t_enegrad=t_enegrad+tcpu()-tt0
1702 ! Compute accelerations from long-range forces
1704 if (large.and. mod(itime,ntwe).eq.0) then
1705 write (iout,*) "energia_long",energia_long(0)
1706 write (iout,*) "Cartesian and internal coordinates: step 2"
1708 call pdbout(0.0d0,"cipiszcze",iout)
1710 write (iout,*) "Accelerations from long-range forces"
1712 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1713 (d_a(j,i+nres),j=1,3)
1715 write (iout,*) "Velocities, step 2"
1717 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1718 (d_t(j,i+nres),j=1,3)
1722 ! Compute the final RESPA step (increment velocities)
1723 ! write (iout,*) "*********************** RESPA fin"
1725 ! Compute the complete potential energy
1727 potEcomp(i)=energia_short(i)+energia_long(i)
1729 potE=potEcomp(0)-potEcomp(20)
1730 ! potE=energia_short(0)+energia_long(0)
1733 ! Calculate the kinetic and the total energy and the kinetic temperature
1736 ! Couple the system to Berendsen bath if needed
1737 if (tbf .and. lang.eq.0) then
1740 kinetic_T=2.0d0/(dimen3*Rb)*EK
1741 ! Backup the coordinates, velocities, and accelerations
1743 if (mod(itime,ntwe).eq.0 .and. large) then
1744 write (iout,*) "Velocities, end"
1746 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1747 (d_t(j,i+nres),j=1,3)
1752 end subroutine RESPA_step
1753 !-----------------------------------------------------------------------------
1754 subroutine RESPA_vel
1755 ! First and last RESPA step (incrementing velocities using long-range
1758 ! implicit real*8 (a-h,o-z)
1759 ! include 'DIMENSIONS'
1760 ! include 'COMMON.CONTROL'
1761 ! include 'COMMON.VAR'
1762 ! include 'COMMON.MD'
1763 ! include 'COMMON.CHAIN'
1764 ! include 'COMMON.DERIV'
1765 ! include 'COMMON.GEO'
1766 ! include 'COMMON.LOCAL'
1767 ! include 'COMMON.INTERACT'
1768 ! include 'COMMON.IOUNITS'
1769 ! include 'COMMON.NAMES'
1770 integer :: i,j,inres
1773 d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
1777 d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
1781 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
1784 d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
1789 end subroutine RESPA_vel
1790 !-----------------------------------------------------------------------------
1792 ! Applying velocity Verlet algorithm - step 1 to coordinates
1794 ! implicit real*8 (a-h,o-z)
1795 ! include 'DIMENSIONS'
1796 ! include 'COMMON.CONTROL'
1797 ! include 'COMMON.VAR'
1798 ! include 'COMMON.MD'
1799 ! include 'COMMON.CHAIN'
1800 ! include 'COMMON.DERIV'
1801 ! include 'COMMON.GEO'
1802 ! include 'COMMON.LOCAL'
1803 ! include 'COMMON.INTERACT'
1804 ! include 'COMMON.IOUNITS'
1805 ! include 'COMMON.NAMES'
1806 real(kind=8) :: adt,adt2
1807 integer :: i,j,inres
1810 write (iout,*) "VELVERLET1 START: DC"
1812 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),&
1813 (dc(j,i+nres),j=1,3)
1817 adt=d_a_old(j,0)*d_time
1819 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1820 d_t_new(j,0)=d_t_old(j,0)+adt2
1821 d_t(j,0)=d_t_old(j,0)+adt
1825 adt=d_a_old(j,i)*d_time
1827 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1828 d_t_new(j,i)=d_t_old(j,i)+adt2
1829 d_t(j,i)=d_t_old(j,i)+adt
1833 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
1836 adt=d_a_old(j,inres)*d_time
1838 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1839 d_t_new(j,inres)=d_t_old(j,inres)+adt2
1840 d_t(j,inres)=d_t_old(j,inres)+adt
1845 write (iout,*) "VELVERLET1 END: DC"
1847 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),&
1848 (dc(j,i+nres),j=1,3)
1852 end subroutine verlet1
1853 !-----------------------------------------------------------------------------
1855 ! Step 2 of the velocity Verlet algorithm: update velocities
1857 ! implicit real*8 (a-h,o-z)
1858 ! include 'DIMENSIONS'
1859 ! include 'COMMON.CONTROL'
1860 ! include 'COMMON.VAR'
1861 ! include 'COMMON.MD'
1862 ! include 'COMMON.CHAIN'
1863 ! include 'COMMON.DERIV'
1864 ! include 'COMMON.GEO'
1865 ! include 'COMMON.LOCAL'
1866 ! include 'COMMON.INTERACT'
1867 ! include 'COMMON.IOUNITS'
1868 ! include 'COMMON.NAMES'
1869 integer :: i,j,inres
1872 d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
1876 d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
1880 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
1883 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
1888 end subroutine verlet2
1889 !-----------------------------------------------------------------------------
1890 subroutine sddir_precalc
1891 ! Applying velocity Verlet algorithm - step 1 to coordinates
1892 ! implicit real*8 (a-h,o-z)
1893 ! include 'DIMENSIONS'
1899 ! include 'COMMON.CONTROL'
1900 ! include 'COMMON.VAR'
1901 ! include 'COMMON.MD'
1903 ! include 'COMMON.LANGEVIN'
1905 ! include 'COMMON.LANGEVIN.lang0'
1907 ! include 'COMMON.CHAIN'
1908 ! include 'COMMON.DERIV'
1909 ! include 'COMMON.GEO'
1910 ! include 'COMMON.LOCAL'
1911 ! include 'COMMON.INTERACT'
1912 ! include 'COMMON.IOUNITS'
1913 ! include 'COMMON.NAMES'
1914 ! include 'COMMON.TIME1'
1915 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1916 !el common /stochcalc/ stochforcvec
1917 real(kind=8) :: time00
1919 ! Compute friction and stochastic forces
1924 time_fric=time_fric+MPI_Wtime()-time00
1926 call stochastic_force(stochforcvec)
1927 time_stoch=time_stoch+MPI_Wtime()-time00
1930 ! Compute the acceleration due to friction forces (d_af_work) and stochastic
1931 ! forces (d_as_work)
1933 call ginv_mult(fric_work, d_af_work)
1934 call ginv_mult(stochforcvec, d_as_work)
1936 end subroutine sddir_precalc
1937 !-----------------------------------------------------------------------------
1938 subroutine sddir_verlet1
1939 ! Applying velocity Verlet algorithm - step 1 to velocities
1942 ! implicit real*8 (a-h,o-z)
1943 ! include 'DIMENSIONS'
1944 ! include 'COMMON.CONTROL'
1945 ! include 'COMMON.VAR'
1946 ! include 'COMMON.MD'
1948 ! include 'COMMON.LANGEVIN'
1950 ! include 'COMMON.LANGEVIN.lang0'
1952 ! include 'COMMON.CHAIN'
1953 ! include 'COMMON.DERIV'
1954 ! include 'COMMON.GEO'
1955 ! include 'COMMON.LOCAL'
1956 ! include 'COMMON.INTERACT'
1957 ! include 'COMMON.IOUNITS'
1958 ! include 'COMMON.NAMES'
1959 ! Revised 3/31/05 AL: correlation between random contributions to
1960 ! position and velocity increments included.
1961 real(kind=8) :: sqrt13 = 0.57735026918962576451d0 ! 1/sqrt(3)
1962 real(kind=8) :: adt,adt2
1963 integer :: i,j,ind,inres
1965 ! Add the contribution from BOTH friction and stochastic force to the
1966 ! coordinates, but ONLY the contribution from the friction forces to velocities
1969 adt=(d_a_old(j,0)+d_af_work(j))*d_time
1970 adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
1971 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1972 d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
1973 d_t(j,0)=d_t_old(j,0)+adt
1978 adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
1979 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1980 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1981 d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
1982 d_t(j,i)=d_t_old(j,i)+adt
1987 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
1990 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
1991 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1992 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1993 d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
1994 d_t(j,inres)=d_t_old(j,inres)+adt
2000 end subroutine sddir_verlet1
2001 !-----------------------------------------------------------------------------
2002 subroutine sddir_verlet2
2003 ! Calculating the adjusted velocities for accelerations
2006 ! implicit real*8 (a-h,o-z)
2007 ! include 'DIMENSIONS'
2008 ! include 'COMMON.CONTROL'
2009 ! include 'COMMON.VAR'
2010 ! include 'COMMON.MD'
2012 ! include 'COMMON.LANGEVIN'
2014 ! include 'COMMON.LANGEVIN.lang0'
2016 ! include 'COMMON.CHAIN'
2017 ! include 'COMMON.DERIV'
2018 ! include 'COMMON.GEO'
2019 ! include 'COMMON.LOCAL'
2020 ! include 'COMMON.INTERACT'
2021 ! include 'COMMON.IOUNITS'
2022 ! include 'COMMON.NAMES'
2023 real(kind=8),dimension(6*nres) :: stochforcvec,d_as_work1 !(MAXRES6) maxres6=6*maxres
2024 real(kind=8) :: cos60 = 0.5d0, sin60 = 0.86602540378443864676d0
2025 integer :: i,j,ind,inres
2026 ! Revised 3/31/05 AL: correlation between random contributions to
2027 ! position and velocity increments included.
2028 ! The correlation coefficients are calculated at low-friction limit.
2029 ! Also, friction forces are now not calculated with new velocities.
2031 ! call friction_force
2032 call stochastic_force(stochforcvec)
2034 ! Compute the acceleration due to friction forces (d_af_work) and stochastic
2035 ! forces (d_as_work)
2037 call ginv_mult(stochforcvec, d_as_work1)
2043 d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j)) &
2044 +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
2049 d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j)) &
2050 +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
2055 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2058 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres) &
2059 +d_af_work(ind+j))+sin60*d_as_work(ind+j) &
2060 +cos60*d_as_work1(ind+j))*d_time
2066 end subroutine sddir_verlet2
2067 !-----------------------------------------------------------------------------
2068 subroutine max_accel
2070 ! Find the maximum difference in the accelerations of the the sites
2071 ! at the beginning and the end of the time step.
2074 ! implicit real*8 (a-h,o-z)
2075 ! include 'DIMENSIONS'
2076 ! include 'COMMON.CONTROL'
2077 ! include 'COMMON.VAR'
2078 ! include 'COMMON.MD'
2079 ! include 'COMMON.CHAIN'
2080 ! include 'COMMON.DERIV'
2081 ! include 'COMMON.GEO'
2082 ! include 'COMMON.LOCAL'
2083 ! include 'COMMON.INTERACT'
2084 ! include 'COMMON.IOUNITS'
2085 real(kind=8),dimension(3) :: aux,accel,accel_old
2086 real(kind=8) :: dacc
2090 ! aux(j)=d_a(j,0)-d_a_old(j,0)
2091 accel_old(j)=d_a_old(j,0)
2098 ! 7/3/08 changed to asymmetric difference
2100 ! accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
2101 accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
2102 accel(j)=accel(j)+0.5d0*d_a(j,i)
2103 ! if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
2104 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
2105 dacc=dabs(accel(j)-accel_old(j))
2106 ! write (iout,*) i,dacc
2107 if (dacc.gt.amax) amax=dacc
2115 accel_old(j)=d_a_old(j,0)
2120 accel_old(j)=accel_old(j)+d_a_old(j,1)
2121 accel(j)=accel(j)+d_a(j,1)
2125 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2127 ! accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
2128 accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
2129 accel(j)=accel(j)+d_a(j,i+nres)
2133 ! if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
2134 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
2135 dacc=dabs(accel(j)-accel_old(j))
2136 ! write (iout,*) "side-chain",i,dacc
2137 if (dacc.gt.amax) amax=dacc
2141 accel_old(j)=accel_old(j)+d_a_old(j,i)
2142 accel(j)=accel(j)+d_a(j,i)
2143 ! aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
2147 end subroutine max_accel
2148 !-----------------------------------------------------------------------------
2149 subroutine predict_edrift(epdrift)
2151 ! Predict the drift of the potential energy
2154 use control_data, only: lmuca
2155 ! implicit real*8 (a-h,o-z)
2156 ! include 'DIMENSIONS'
2157 ! include 'COMMON.CONTROL'
2158 ! include 'COMMON.VAR'
2159 ! include 'COMMON.MD'
2160 ! include 'COMMON.CHAIN'
2161 ! include 'COMMON.DERIV'
2162 ! include 'COMMON.GEO'
2163 ! include 'COMMON.LOCAL'
2164 ! include 'COMMON.INTERACT'
2165 ! include 'COMMON.IOUNITS'
2166 ! include 'COMMON.MUCA'
2167 real(kind=8) :: epdrift,epdriftij
2169 ! Drift of the potential energy
2175 epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
2176 if (lmuca) epdriftij=epdriftij*factor
2177 ! write (iout,*) "back",i,j,epdriftij
2178 if (epdriftij.gt.epdrift) epdrift=epdriftij
2182 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2185 dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
2186 if (lmuca) epdriftij=epdriftij*factor
2187 ! write (iout,*) "side",i,j,epdriftij
2188 if (epdriftij.gt.epdrift) epdrift=epdriftij
2192 epdrift=0.5d0*epdrift*d_time*d_time
2193 ! write (iout,*) "epdrift",epdrift
2195 end subroutine predict_edrift
2196 !-----------------------------------------------------------------------------
2197 subroutine verlet_bath
2199 ! Coupling to the thermostat by using the Berendsen algorithm
2202 ! implicit real*8 (a-h,o-z)
2203 ! include 'DIMENSIONS'
2204 ! include 'COMMON.CONTROL'
2205 ! include 'COMMON.VAR'
2206 ! include 'COMMON.MD'
2207 ! include 'COMMON.CHAIN'
2208 ! include 'COMMON.DERIV'
2209 ! include 'COMMON.GEO'
2210 ! include 'COMMON.LOCAL'
2211 ! include 'COMMON.INTERACT'
2212 ! include 'COMMON.IOUNITS'
2213 ! include 'COMMON.NAMES'
2214 real(kind=8) :: T_half,fact
2215 integer :: i,j,inres
2217 T_half=2.0d0/(dimen3*Rb)*EK
2218 fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
2219 ! write(iout,*) "T_half", T_half
2220 ! write(iout,*) "EK", EK
2221 ! write(iout,*) "fact", fact
2223 d_t(j,0)=fact*d_t(j,0)
2227 d_t(j,i)=fact*d_t(j,i)
2231 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2234 d_t(j,inres)=fact*d_t(j,inres)
2239 end subroutine verlet_bath
2240 !-----------------------------------------------------------------------------
2242 ! Set up the initial conditions of a MD simulation
2245 use control, only:tcpu
2246 !el use io_basic, only:ilen
2249 use minimm, only:minim_dc,minimize,sc_move
2250 use io_config, only:readrst
2251 use io, only:statout
2252 ! implicit real*8 (a-h,o-z)
2253 ! include 'DIMENSIONS'
2256 character(len=16) :: form
2257 integer :: IERROR,ERRCODE
2259 ! include 'COMMON.SETUP'
2260 ! include 'COMMON.CONTROL'
2261 ! include 'COMMON.VAR'
2262 ! include 'COMMON.MD'
2264 ! include 'COMMON.LANGEVIN'
2266 ! include 'COMMON.LANGEVIN.lang0'
2268 ! include 'COMMON.CHAIN'
2269 ! include 'COMMON.DERIV'
2270 ! include 'COMMON.GEO'
2271 ! include 'COMMON.LOCAL'
2272 ! include 'COMMON.INTERACT'
2273 ! include 'COMMON.IOUNITS'
2274 ! include 'COMMON.NAMES'
2275 ! include 'COMMON.REMD'
2276 real(kind=8),dimension(0:n_ene) :: energia_long,energia_short
2277 real(kind=8),dimension(3) :: vcm,incr,L
2278 real(kind=8) :: xv,sigv,lowb,highb
2279 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
2280 character(len=256) :: qstr
2283 character(len=50) :: tytul
2284 logical :: file_exist
2285 !el common /gucio/ cm
2286 integer :: i,j,ipos,iq,iw,nft_sc,iretcode,nfun,itime,ierr
2287 real(kind=8) :: etot,tt0
2291 ! write(iout,*) "d_time", d_time
2292 ! Compute the standard deviations of stochastic forces for Langevin dynamics
2293 ! if the friction coefficients do not depend on surface area
2294 if (lang.gt.0 .and. .not.surfarea) then
2296 stdforcp(i)=stdfp*dsqrt(gamp)
2299 stdforcsc(i)=stdfsc(iabs(itype(i,1))) &
2300 *dsqrt(gamsc(iabs(itype(i,1))))
2303 ! Open the pdb file for snapshotshots
2306 if (ilen(tmpdir).gt.0) &
2307 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2308 liczba(:ilen(liczba))//".pdb")
2310 file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2314 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) &
2315 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2316 liczba(:ilen(liczba))//".x")
2317 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2320 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) &
2321 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2322 liczba(:ilen(liczba))//".cx")
2323 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2329 if (ilen(tmpdir).gt.0) &
2330 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
2331 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
2333 if (ilen(tmpdir).gt.0) &
2334 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
2335 cartname=prefix(:ilen(prefix))//"_MD.cx"
2339 write (qstr,'(256(1h ))')
2342 iq = qinfrag(i,iset)*10
2343 iw = wfrag(i,iset)/100
2345 if(me.eq.king.or..not.out1file) &
2346 write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
2347 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
2352 iq = qinpair(i,iset)*10
2353 iw = wpair(i,iset)/100
2355 if(me.eq.king.or..not.out1file) &
2356 write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
2357 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
2361 ! pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
2363 ! cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
2365 ! cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
2367 ! statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
2371 if (restart1file) then
2373 inquire(file=mremd_rst_name,exist=file_exist)
2374 write (*,*) me," Before broadcast: file_exist",file_exist
2376 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,&
2379 write (*,*) me," After broadcast: file_exist",file_exist
2380 ! inquire(file=mremd_rst_name,exist=file_exist)
2381 if(me.eq.king.or..not.out1file) &
2382 write(iout,*) "Initial state read by master and distributed"
2384 if (ilen(tmpdir).gt.0) &
2385 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_' &
2386 //liczba(:ilen(liczba))//'.rst')
2387 inquire(file=rest2name,exist=file_exist)
2390 if(.not.restart1file) then
2391 if(me.eq.king.or..not.out1file) &
2392 write(iout,*) "Initial state will be read from file ",&
2393 rest2name(:ilen(rest2name))
2396 call rescale_weights(t_bath)
2398 if(me.eq.king.or..not.out1file)then
2399 if (restart1file) then
2400 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),&
2403 write(iout,*) "File ",rest2name(:ilen(rest2name)),&
2406 write(iout,*) "Initial velocities randomly generated"
2413 ! Generate initial velocities
2414 if(me.eq.king.or..not.out1file) &
2415 write(iout,*) "Initial velocities randomly generated"
2420 ! rest2name = prefix(:ilen(prefix))//'.rst'
2421 if(me.eq.king.or..not.out1file)then
2422 write (iout,*) "Initial velocities"
2424 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2425 (d_t(j,i+nres),j=1,3)
2427 ! Zeroing the total angular momentum of the system
2428 write(iout,*) "Calling the zero-angular momentum subroutine"
2431 ! Getting the potential energy and forces and velocities and accelerations
2433 ! write (iout,*) "velocity of the center of the mass:"
2434 ! write (iout,*) (vcm(j),j=1,3)
2436 d_t(j,0)=d_t(j,0)-vcm(j)
2438 ! Removing the velocity of the center of mass
2440 if(me.eq.king.or..not.out1file)then
2441 write (iout,*) "vcm right after adjustment:"
2442 write (iout,*) (vcm(j),j=1,3)
2444 if ((.not.rest).and.(indpdb.eq.0)) then
2446 if(iranconf.ne.0) then
2448 print *, 'Calling OVERLAP_SC'
2449 call overlap_sc(fail)
2452 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
2453 print *,'SC_move',nft_sc,etot
2454 if(me.eq.king.or..not.out1file) &
2455 write(iout,*) 'SC_move',nft_sc,etot
2459 print *, 'Calling MINIM_DC'
2460 call minim_dc(etot,iretcode,nfun)
2462 call geom_to_var(nvar,varia)
2463 print *,'Calling MINIMIZE.'
2464 call minimize(etot,varia,iretcode,nfun)
2465 call var_to_geom(nvar,varia)
2467 if(me.eq.king.or..not.out1file) &
2468 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2471 call chainbuild_cart
2476 kinetic_T=2.0d0/(dimen3*Rb)*EK
2477 if(me.eq.king.or..not.out1file)then
2487 call etotal(potEcomp)
2488 if (large) call enerprint(potEcomp)
2491 t_etotal=t_etotal+MPI_Wtime()-tt0
2493 t_etotal=t_etotal+tcpu()-tt0
2500 if (amax*d_time .gt. dvmax) then
2501 d_time=d_time*dvmax/amax
2502 if(me.eq.king.or..not.out1file) write (iout,*) &
2503 "Time step reduced to",d_time,&
2504 " because of too large initial acceleration."
2506 if(me.eq.king.or..not.out1file)then
2507 write(iout,*) "Potential energy and its components"
2508 call enerprint(potEcomp)
2509 ! write(iout,*) (potEcomp(i),i=0,n_ene)
2511 potE=potEcomp(0)-potEcomp(20)
2514 if (ntwe.ne.0) call statout(itime)
2515 if(me.eq.king.or..not.out1file) &
2516 write (iout,'(/a/3(a25,1pe14.5/))') "Initial:", &
2517 " Kinetic energy",EK," Potential energy",potE, &
2518 " Total energy",totE," Maximum acceleration ", &
2521 write (iout,*) "Initial coordinates"
2523 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),&
2526 write (iout,*) "Initial dC"
2528 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),&
2529 (dc(j,i+nres),j=1,3)
2531 write (iout,*) "Initial velocities"
2532 write (iout,"(13x,' backbone ',23x,' side chain')")
2534 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2535 (d_t(j,i+nres),j=1,3)
2537 write (iout,*) "Initial accelerations"
2539 ! write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2540 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),&
2541 (d_a(j,i+nres),j=1,3)
2547 d_t_old(j,i)=d_t(j,i)
2548 d_a_old(j,i)=d_a(j,i)
2550 ! write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
2559 call etotal_short(energia_short)
2560 if (large) call enerprint(potEcomp)
2563 t_eshort=t_eshort+MPI_Wtime()-tt0
2565 t_eshort=t_eshort+tcpu()-tt0
2570 if(.not.out1file .and. large) then
2571 write (iout,*) "energia_long",energia_long(0),&
2572 " energia_short",energia_short(0),&
2573 " total",energia_long(0)+energia_short(0)
2574 write (iout,*) "Initial fast-force accelerations"
2576 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
2577 (d_a(j,i+nres),j=1,3)
2580 ! 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
2583 d_a_short(j,i)=d_a(j,i)
2592 call etotal_long(energia_long)
2593 if (large) call enerprint(potEcomp)
2596 t_elong=t_elong+MPI_Wtime()-tt0
2598 t_elong=t_elong+tcpu()-tt0
2603 if(.not.out1file .and. large) then
2604 write (iout,*) "energia_long",energia_long(0)
2605 write (iout,*) "Initial slow-force accelerations"
2607 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
2608 (d_a(j,i+nres),j=1,3)
2612 t_enegrad=t_enegrad+MPI_Wtime()-tt0
2614 t_enegrad=t_enegrad+tcpu()-tt0
2618 end subroutine init_MD
2619 !-----------------------------------------------------------------------------
2620 subroutine random_vel
2622 ! implicit real*8 (a-h,o-z)
2624 use random, only:anorm_distr
2626 ! include 'DIMENSIONS'
2627 ! include 'COMMON.CONTROL'
2628 ! include 'COMMON.VAR'
2629 ! include 'COMMON.MD'
2631 ! include 'COMMON.LANGEVIN'
2633 ! include 'COMMON.LANGEVIN.lang0'
2635 ! include 'COMMON.CHAIN'
2636 ! include 'COMMON.DERIV'
2637 ! include 'COMMON.GEO'
2638 ! include 'COMMON.LOCAL'
2639 ! include 'COMMON.INTERACT'
2640 ! include 'COMMON.IOUNITS'
2641 ! include 'COMMON.NAMES'
2642 ! include 'COMMON.TIME1'
2643 real(kind=8) :: xv,sigv,lowb,highb ,Ek1
2645 real(kind=8) ,allocatable, dimension(:) :: DDU1,DDU2,DL2,DL1,xsolv,DML,rs
2646 real(kind=8) :: sumx
2648 real(kind=8) ,allocatable, dimension(:) :: rsold
2649 real (kind=8),allocatable,dimension(:,:) :: matold
2652 integer :: i,j,ii,k,ind,mark,imark
2653 ! Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
2654 ! First generate velocities in the eigenspace of the G matrix
2655 ! write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
2662 sigv=dsqrt((Rb*t_bath)/geigen(i))
2665 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
2667 write (iout,*) "i",i," ii",ii," geigen",geigen(i),&
2668 " d_t_work_new",d_t_work_new(ii)
2679 Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
2682 write (iout,*) "Ek from eigenvectors",Ek1
2683 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
2687 ! Transform velocities to UNRES coordinate space
2688 allocate (DL1(2*nres))
2689 allocate (DDU1(2*nres))
2690 allocate (DL2(2*nres))
2691 allocate (DDU2(2*nres))
2692 allocate (xsolv(2*nres))
2693 allocate (DML(2*nres))
2694 allocate (rs(2*nres))
2696 allocate (rsold(2*nres))
2697 allocate (matold(2*nres,2*nres))
2699 matold(1,1)=DMorig(1)
2700 matold(1,2)=DU1orig(1)
2701 matold(1,3)=DU2orig(1)
2702 write (*,*) DMorig(1),DU1orig(1),DU2orig(1)
2707 matold(i,j)=DMorig(i)
2708 matold(i,j-1)=DU1orig(i-1)
2710 matold(i,j-2)=DU2orig(i-2)
2718 matold(i,j+1)=DU1orig(i)
2724 matold(i,j+2)=DU2orig(i)
2728 matold(dimen,dimen)=DMorig(dimen)
2729 matold(dimen,dimen-1)=DU1orig(dimen-1)
2730 matold(dimen,dimen-2)=DU2orig(dimen-2)
2731 write(iout,*) "old gmatrix"
2732 call matout(dimen,dimen,2*nres,2*nres,matold)
2736 ! Find the ith eigenvector of the pentadiagonal inertiq matrix
2740 DML(j)=DMorig(j)-geigen(i)
2743 DML(j-1)=DMorig(j)-geigen(i)
2748 DDU1(imark-1)=DU2orig(imark-1)
2749 do j=imark+1,dimen-1
2750 DDU1(j-1)=DU1orig(j)
2758 DDU2(j)=DU2orig(j+1)
2767 write (iout,*) "DL2,DL1,DML,DDU1,DDU2"
2768 write(iout,'(10f10.5)') (DL2(k),k=3,dimen-1)
2769 write(iout,'(10f10.5)') (DL1(k),k=2,dimen-1)
2770 write(iout,'(10f10.5)') (DML(k),k=1,dimen-1)
2771 write(iout,'(10f10.5)') (DDU1(k),k=1,dimen-2)
2772 write(iout,'(10f10.5)') (DDU2(k),k=1,dimen-3)
2775 if (imark.gt.2) rs(imark-2)=-DU2orig(imark-2)
2776 if (imark.gt.1) rs(imark-1)=-DU1orig(imark-1)
2777 if (imark.lt.dimen) rs(imark)=-DU1orig(imark)
2778 if (imark.lt.dimen-1) rs(imark+1)=-DU2orig(imark)
2782 ! write (iout,*) "Vector rs"
2784 ! write (iout,*) j,rs(j)
2787 call FDIAG(dimen-1,DL2,DL1,DML,DDU1,DDU2,rs,xsolv,mark)
2794 sumx=-geigen(i)*xsolv(j)
2796 sumx=sumx+matold(j,k)*xsolv(k)
2799 sumx=sumx+matold(j,k)*xsolv(k-1)
2801 write(iout,'(i5,3f10.5)') j,sumx,rsold(j),sumx-rsold(j)
2804 sumx=-geigen(i)*xsolv(j-1)
2806 sumx=sumx+matold(j,k)*xsolv(k)
2809 sumx=sumx+matold(j,k)*xsolv(k-1)
2811 write(iout,'(i5,3f10.5)') j-1,sumx,rsold(j-1),sumx-rsold(j-1)
2815 "Solution of equations system",i," for eigenvalue",geigen(i)
2817 write(iout,'(i5,f10.5)') j,xsolv(j)
2820 do j=dimen-1,imark,-1
2825 write (iout,*) "Un-normalized eigenvector",i," for eigenvalue",geigen(i)
2827 write(iout,'(i5,f10.5)') j,xsolv(j)
2830 ! Normalize ith eigenvector
2833 sumx=sumx+xsolv(j)**2
2837 xsolv(j)=xsolv(j)/sumx
2840 write (iout,*) "Eigenvector",i," for eigenvalue",geigen(i)
2842 write(iout,'(i5,3f10.5)') j,xsolv(j)
2845 ! All done at this point for eigenvector i, exit loop
2853 write (iout,*) "Unable to find eigenvector",i
2856 ! write (iout,*) "i=",i
2858 ! write (iout,*) "k=",k
2861 ! write(iout,*) "ind",ind," ind1",3*(i-1)+k
2862 d_t_work(ind)=d_t_work(ind) &
2863 +xsolv(j)*d_t_work_new(3*(i-1)+k)
2866 enddo ! i (loop over the eigenvectors)
2869 write (iout,*) "d_t_work"
2871 write (iout,"(i5,f10.5)") i,d_t_work(i)
2876 if (itype(i,1).eq.10) then
2883 ! write (iout,*) "k",k," ii+k",ii+k," ii+j+k",ii+j+k,"EK1",Ek1
2884 Ek1=Ek1+0.5d0*mp(1)*((d_t_work(ii+j+k)+d_t_work_new(ii+k))/2)**2+&
2885 0.5d0*0.25d0*IP(1)*(d_t_work(ii+j+k)-d_t_work_new(ii+k))**2
2888 if (itype(i,1).ne.10) ii=ii+3
2889 write (iout,*) "i",i," itype",itype(i,1)," mass",msc(itype(i,1),1)
2890 write (iout,*) "ii",ii
2893 write (iout,*) "k",k," ii",ii,"EK1",EK1
2894 if (iabs(itype(i,1)).ne.10) Ek1=Ek1+0.5d0*Isc(iabs(itype(i,1),1))*(d_t_work(ii)-d_t_work(ii-3))**2
2895 Ek1=Ek1+0.5d0*msc(iabs(itype(i,1)),1)*d_t_work(ii)**2
2897 write (iout,*) "i",i," ii",ii
2899 write (iout,*) "Ek from d_t_work",Ek1
2900 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
2916 d_t(k,j)=d_t_work(ind)
2919 if (itype(j).ne.10 .and. itype(j).ne.ntyp1) then
2921 d_t(k,j+nres)=d_t_work(ind)
2927 write (iout,*) "Random velocities in the Calpha,SC space"
2929 write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3)
2932 write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3)
2939 if (itype(i,1).eq.10) then
2941 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
2945 d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
2946 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
2951 write (iout,*)"Random velocities in the virtual-bond-vector space"
2953 write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3)
2956 write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3)
2959 write (iout,*) "Ek from d_t_work",Ek1
2960 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
2968 d_t_work(ind)=d_t_work(ind) &
2969 +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
2971 ! write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
2975 ! Transfer to the d_t vector
2977 d_t(j,0)=d_t_work(j)
2983 d_t(j,i)=d_t_work(ind)
2987 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2990 d_t(j,i+nres)=d_t_work(ind)
2996 ! write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",&
2997 ! 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
3001 end subroutine random_vel
3002 !-----------------------------------------------------------------------------
3004 subroutine sd_verlet_p_setup
3005 ! Sets up the parameters of stochastic Verlet algorithm
3006 ! implicit real*8 (a-h,o-z)
3007 ! include 'DIMENSIONS'
3008 use control, only: tcpu
3013 ! include 'COMMON.CONTROL'
3014 ! include 'COMMON.VAR'
3015 ! include 'COMMON.MD'
3017 ! include 'COMMON.LANGEVIN'
3019 ! include 'COMMON.LANGEVIN.lang0'
3021 ! include 'COMMON.CHAIN'
3022 ! include 'COMMON.DERIV'
3023 ! include 'COMMON.GEO'
3024 ! include 'COMMON.LOCAL'
3025 ! include 'COMMON.INTERACT'
3026 ! include 'COMMON.IOUNITS'
3027 ! include 'COMMON.NAMES'
3028 ! include 'COMMON.TIME1'
3029 real(kind=8),dimension(6*nres) :: emgdt !(MAXRES6) maxres6=6*maxres
3030 real(kind=8) :: pterm,vterm,rho,rhoc,vsig
3031 real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
3032 prand_vec,vrand_vec1,vrand_vec2 !(MAXRES6) maxres6=6*maxres
3033 logical :: lprn = .false.
3034 real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
3035 real(kind=8) :: ktm,gdt,egdt,gdt2,gdt3,gdt4,gdt5,gdt6,gdt7,gdt8,&
3037 integer :: i,maxres2
3044 ! AL 8/17/04 Code adapted from tinker
3046 ! Get the frictional and random terms for stochastic dynamics in the
3047 ! eigenspace of mass-scaled UNRES friction matrix
3051 gdt = fricgam(i) * d_time
3053 ! Stochastic dynamics reduces to simple MD for zero friction
3055 if (gdt .le. zero) then
3056 pfric_vec(i) = 1.0d0
3057 vfric_vec(i) = d_time
3058 afric_vec(i) = 0.5d0 * d_time * d_time
3059 prand_vec(i) = 0.0d0
3060 vrand_vec1(i) = 0.0d0
3061 vrand_vec2(i) = 0.0d0
3063 ! Analytical expressions when friction coefficient is large
3066 if (gdt .ge. gdt_radius) then
3069 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
3070 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
3071 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
3072 vterm = 1.0d0 - egdt**2
3073 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
3075 ! Use series expansions when friction coefficient is small
3086 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0 &
3087 - gdt5/120.0d0 + gdt6/720.0d0 &
3088 - gdt7/5040.0d0 + gdt8/40320.0d0 &
3089 - gdt9/362880.0d0) / fricgam(i)**2
3090 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
3091 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
3092 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0 &
3093 + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0 &
3094 + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0 &
3095 + 127.0d0*gdt9/90720.0d0
3096 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0 &
3097 - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0 &
3098 - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0 &
3099 - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
3100 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0 &
3101 - 17.0d0*gdt2/1280.0d0 &
3102 + 17.0d0*gdt3/6144.0d0 &
3103 + 40967.0d0*gdt4/34406400.0d0 &
3104 - 57203.0d0*gdt5/275251200.0d0 &
3105 - 1429487.0d0*gdt6/13212057600.0d0)
3108 ! Compute the scaling factors of random terms for the nonzero friction case
3110 ktm = 0.5d0*d_time/fricgam(i)
3111 psig = dsqrt(ktm*pterm) / fricgam(i)
3112 vsig = dsqrt(ktm*vterm)
3113 rhoc = dsqrt(1.0d0 - rho*rho)
3115 vrand_vec1(i) = vsig * rho
3116 vrand_vec2(i) = vsig * rhoc
3121 "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",&
3124 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),&
3125 afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
3129 ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
3132 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
3133 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
3134 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
3135 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
3136 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
3137 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
3140 t_sdsetup=t_sdsetup+MPI_Wtime()
3142 t_sdsetup=t_sdsetup+tcpu()-tt0
3145 end subroutine sd_verlet_p_setup
3146 !-----------------------------------------------------------------------------
3147 subroutine eigtransf1(n,ndim,ab,d,c)
3151 real(kind=8) :: ab(ndim,ndim,n),c(ndim,n),d(ndim)
3157 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
3162 end subroutine eigtransf1
3163 !-----------------------------------------------------------------------------
3164 subroutine eigtransf(n,ndim,a,b,d,c)
3168 real(kind=8) :: a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
3174 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
3179 end subroutine eigtransf
3180 !-----------------------------------------------------------------------------
3181 subroutine sd_verlet1
3183 ! Applying stochastic velocity Verlet algorithm - step 1 to velocities
3185 ! implicit real*8 (a-h,o-z)
3186 ! include 'DIMENSIONS'
3187 ! include 'COMMON.CONTROL'
3188 ! include 'COMMON.VAR'
3189 ! include 'COMMON.MD'
3191 ! include 'COMMON.LANGEVIN'
3193 ! include 'COMMON.LANGEVIN.lang0'
3195 ! include 'COMMON.CHAIN'
3196 ! include 'COMMON.DERIV'
3197 ! include 'COMMON.GEO'
3198 ! include 'COMMON.LOCAL'
3199 ! include 'COMMON.INTERACT'
3200 ! include 'COMMON.IOUNITS'
3201 ! include 'COMMON.NAMES'
3202 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
3203 !el common /stochcalc/ stochforcvec
3204 logical :: lprn = .false.
3205 real(kind=8) :: ddt1,ddt2
3206 integer :: i,j,ind,inres
3208 ! write (iout,*) "dc_old"
3210 ! write (iout,'(i5,3f10.5,5x,3f10.5)')
3211 ! & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3214 dc_work(j)=dc_old(j,0)
3215 d_t_work(j)=d_t_old(j,0)
3216 d_a_work(j)=d_a_old(j,0)
3221 dc_work(ind+j)=dc_old(j,i)
3222 d_t_work(ind+j)=d_t_old(j,i)
3223 d_a_work(ind+j)=d_a_old(j,i)
3228 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3230 dc_work(ind+j)=dc_old(j,i+nres)
3231 d_t_work(ind+j)=d_t_old(j,i+nres)
3232 d_a_work(ind+j)=d_a_old(j,i+nres)
3240 "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",&
3244 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
3245 vfric_mat(i,j),afric_mat(i,j),&
3246 prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3254 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) &
3255 +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3256 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3257 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3259 d_t_work_new(i)=ddt1+0.5d0*ddt2
3260 d_t_work(i)=ddt1+ddt2
3265 d_t(j,0)=d_t_work(j)
3270 dc(j,i)=dc_work(ind+j)
3271 d_t(j,i)=d_t_work(ind+j)
3276 if (itype(i,1).ne.10) then
3279 dc(j,inres)=dc_work(ind+j)
3280 d_t(j,inres)=d_t_work(ind+j)
3286 end subroutine sd_verlet1
3287 !-----------------------------------------------------------------------------
3288 subroutine sd_verlet2
3290 ! Calculating the adjusted velocities for accelerations
3292 ! implicit real*8 (a-h,o-z)
3293 ! include 'DIMENSIONS'
3294 ! include 'COMMON.CONTROL'
3295 ! include 'COMMON.VAR'
3296 ! include 'COMMON.MD'
3298 ! include 'COMMON.LANGEVIN'
3300 ! include 'COMMON.LANGEVIN.lang0'
3302 ! include 'COMMON.CHAIN'
3303 ! include 'COMMON.DERIV'
3304 ! include 'COMMON.GEO'
3305 ! include 'COMMON.LOCAL'
3306 ! include 'COMMON.INTERACT'
3307 ! include 'COMMON.IOUNITS'
3308 ! include 'COMMON.NAMES'
3309 !el real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV !(MAXRES6) maxres6=6*maxres
3310 real(kind=8),dimension(6*nres) :: stochforcvecV !(MAXRES6) maxres6=6*maxres
3311 !el common /stochcalc/ stochforcvec
3313 real(kind=8) :: ddt1,ddt2
3314 integer :: i,j,ind,inres
3315 ! Compute the stochastic forces which contribute to velocity change
3317 call stochastic_force(stochforcvecV)
3324 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3325 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+ &
3326 vrand_mat2(i,j)*stochforcvecV(j)
3328 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
3332 d_t(j,0)=d_t_work(j)
3337 d_t(j,i)=d_t_work(ind+j)
3342 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3345 d_t(j,inres)=d_t_work(ind+j)
3351 end subroutine sd_verlet2
3352 !-----------------------------------------------------------------------------
3353 subroutine sd_verlet_ciccotti_setup
3355 ! Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
3357 ! implicit real*8 (a-h,o-z)
3358 ! include 'DIMENSIONS'
3359 use control, only: tcpu
3364 ! include 'COMMON.CONTROL'
3365 ! include 'COMMON.VAR'
3366 ! include 'COMMON.MD'
3368 ! include 'COMMON.LANGEVIN'
3370 ! include 'COMMON.LANGEVIN.lang0'
3372 ! include 'COMMON.CHAIN'
3373 ! include 'COMMON.DERIV'
3374 ! include 'COMMON.GEO'
3375 ! include 'COMMON.LOCAL'
3376 ! include 'COMMON.INTERACT'
3377 ! include 'COMMON.IOUNITS'
3378 ! include 'COMMON.NAMES'
3379 ! include 'COMMON.TIME1'
3380 real(kind=8),dimension(6*nres) :: emgdt !(MAXRES6) maxres6=6*maxres
3381 real(kind=8) :: pterm,vterm,rho,rhoc,vsig
3382 real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
3383 prand_vec,vrand_vec1,vrand_vec2 !(MAXRES6) maxres6=6*maxres
3384 logical :: lprn = .false.
3385 real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
3386 real(kind=8) :: ktm,gdt,egdt,tt0
3387 integer :: i,maxres2
3394 ! AL 8/17/04 Code adapted from tinker
3396 ! Get the frictional and random terms for stochastic dynamics in the
3397 ! eigenspace of mass-scaled UNRES friction matrix
3401 write (iout,*) "i",i," fricgam",fricgam(i)
3402 gdt = fricgam(i) * d_time
3404 ! Stochastic dynamics reduces to simple MD for zero friction
3406 if (gdt .le. zero) then
3407 pfric_vec(i) = 1.0d0
3408 vfric_vec(i) = d_time
3409 afric_vec(i) = 0.5d0*d_time*d_time
3410 prand_vec(i) = afric_vec(i)
3411 vrand_vec2(i) = vfric_vec(i)
3413 ! Analytical expressions when friction coefficient is large
3418 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
3419 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
3420 prand_vec(i) = afric_vec(i)
3421 vrand_vec2(i) = vfric_vec(i)
3423 ! Compute the scaling factors of random terms for the nonzero friction case
3425 ! ktm = 0.5d0*d_time/fricgam(i)
3426 ! psig = dsqrt(ktm*pterm) / fricgam(i)
3427 ! vsig = dsqrt(ktm*vterm)
3428 ! prand_vec(i) = psig*afric_vec(i)
3429 ! vrand_vec2(i) = vsig*vfric_vec(i)
3434 "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",&
3437 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),&
3438 afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
3442 ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
3444 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
3445 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
3446 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
3447 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
3448 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
3450 t_sdsetup=t_sdsetup+MPI_Wtime()
3452 t_sdsetup=t_sdsetup+tcpu()-tt0
3455 end subroutine sd_verlet_ciccotti_setup
3456 !-----------------------------------------------------------------------------
3457 subroutine sd_verlet1_ciccotti
3459 ! Applying stochastic velocity Verlet algorithm - step 1 to velocities
3460 ! implicit real*8 (a-h,o-z)
3462 ! include 'DIMENSIONS'
3466 ! include 'COMMON.CONTROL'
3467 ! include 'COMMON.VAR'
3468 ! include 'COMMON.MD'
3470 ! include 'COMMON.LANGEVIN'
3472 ! include 'COMMON.LANGEVIN.lang0'
3474 ! include 'COMMON.CHAIN'
3475 ! include 'COMMON.DERIV'
3476 ! include 'COMMON.GEO'
3477 ! include 'COMMON.LOCAL'
3478 ! include 'COMMON.INTERACT'
3479 ! include 'COMMON.IOUNITS'
3480 ! include 'COMMON.NAMES'
3481 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
3482 !el common /stochcalc/ stochforcvec
3483 logical :: lprn = .false.
3484 real(kind=8) :: ddt1,ddt2
3485 integer :: i,j,ind,inres
3486 ! write (iout,*) "dc_old"
3488 ! write (iout,'(i5,3f10.5,5x,3f10.5)')
3489 ! & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3492 dc_work(j)=dc_old(j,0)
3493 d_t_work(j)=d_t_old(j,0)
3494 d_a_work(j)=d_a_old(j,0)
3499 dc_work(ind+j)=dc_old(j,i)
3500 d_t_work(ind+j)=d_t_old(j,i)
3501 d_a_work(ind+j)=d_a_old(j,i)
3506 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3508 dc_work(ind+j)=dc_old(j,i+nres)
3509 d_t_work(ind+j)=d_t_old(j,i+nres)
3510 d_a_work(ind+j)=d_a_old(j,i+nres)
3519 "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",&
3523 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
3524 vfric_mat(i,j),afric_mat(i,j),&
3525 prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3533 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) &
3534 +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3535 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3536 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3538 d_t_work_new(i)=ddt1+0.5d0*ddt2
3539 d_t_work(i)=ddt1+ddt2
3544 d_t(j,0)=d_t_work(j)
3549 dc(j,i)=dc_work(ind+j)
3550 d_t(j,i)=d_t_work(ind+j)
3555 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3558 dc(j,inres)=dc_work(ind+j)
3559 d_t(j,inres)=d_t_work(ind+j)
3565 end subroutine sd_verlet1_ciccotti
3566 !-----------------------------------------------------------------------------
3567 subroutine sd_verlet2_ciccotti
3569 ! Calculating the adjusted velocities for accelerations
3571 ! implicit real*8 (a-h,o-z)
3572 ! include 'DIMENSIONS'
3573 ! include 'COMMON.CONTROL'
3574 ! include 'COMMON.VAR'
3575 ! include 'COMMON.MD'
3577 ! include 'COMMON.LANGEVIN'
3579 ! include 'COMMON.LANGEVIN.lang0'
3581 ! include 'COMMON.CHAIN'
3582 ! include 'COMMON.DERIV'
3583 ! include 'COMMON.GEO'
3584 ! include 'COMMON.LOCAL'
3585 ! include 'COMMON.INTERACT'
3586 ! include 'COMMON.IOUNITS'
3587 ! include 'COMMON.NAMES'
3588 !el real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV !(MAXRES6) maxres6=6*maxres
3589 real(kind=8),dimension(6*nres) :: stochforcvecV !(MAXRES6) maxres6=6*maxres
3590 !el common /stochcalc/ stochforcvec
3591 real(kind=8) :: ddt1,ddt2
3592 integer :: i,j,ind,inres
3594 ! Compute the stochastic forces which contribute to velocity change
3596 call stochastic_force(stochforcvecV)
3603 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3604 ! ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
3605 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
3607 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
3611 d_t(j,0)=d_t_work(j)
3616 d_t(j,i)=d_t_work(ind+j)
3621 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3624 d_t(j,inres)=d_t_work(ind+j)
3630 end subroutine sd_verlet2_ciccotti
3632 !-----------------------------------------------------------------------------
3634 !-----------------------------------------------------------------------------
3635 subroutine inertia_tensor
3637 ! Calculating the intertia tensor for the entire protein in order to
3638 ! remove the perpendicular components of velocity matrix which cause
3639 ! the molecule to rotate.
3642 ! implicit real*8 (a-h,o-z)
3643 ! include 'DIMENSIONS'
3644 ! include 'COMMON.CONTROL'
3645 ! include 'COMMON.VAR'
3646 ! include 'COMMON.MD'
3647 ! include 'COMMON.CHAIN'
3648 ! include 'COMMON.DERIV'
3649 ! include 'COMMON.GEO'
3650 ! include 'COMMON.LOCAL'
3651 ! include 'COMMON.INTERACT'
3652 ! include 'COMMON.IOUNITS'
3653 ! include 'COMMON.NAMES'
3655 real(kind=8),dimension(3,3) :: Im,Imcp,eigvec,Id
3656 real(kind=8),dimension(3) :: pr,eigval,L,vp,vrot
3657 real(kind=8) :: M_SC,mag,mag2
3658 real(kind=8),dimension(3,0:nres) :: vpp !(3,0:MAXRES)
3659 real(kind=8),dimension(3) :: vs_p,pp,incr,v
3660 real(kind=8),dimension(3,3) :: pr1,pr2
3662 !el common /gucio/ cm
3663 integer :: iti,inres,i,j,k
3674 ! calculating the center of the mass of the protein
3677 cm(j)=cm(j)+c(j,i)+0.5d0*dc(j,i)
3685 iti=iabs(itype(i,1))
3686 M_SC=M_SC+msc(iabs(iti),1)
3689 cm(j)=cm(j)+msc(iabs(iti),1)*c(j,inres)
3693 cm(j)=cm(j)/(M_SC+(nct-nnt)*mp(1))
3698 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
3700 Im(1,1)=Im(1,1)+mp(1)*(pr(2)*pr(2)+pr(3)*pr(3))
3701 Im(1,2)=Im(1,2)-mp(1)*pr(1)*pr(2)
3702 Im(1,3)=Im(1,3)-mp(1)*pr(1)*pr(3)
3703 Im(2,3)=Im(2,3)-mp(1)*pr(2)*pr(3)
3704 Im(2,2)=Im(2,2)+mp(1)*(pr(3)*pr(3)+pr(1)*pr(1))
3705 Im(3,3)=Im(3,3)+mp(1)*(pr(1)*pr(1)+pr(2)*pr(2))
3709 iti=iabs(itype(i,1))
3712 pr(j)=c(j,inres)-cm(j)
3714 Im(1,1)=Im(1,1)+msc(iabs(iti),1)*(pr(2)*pr(2)+pr(3)*pr(3))
3715 Im(1,2)=Im(1,2)-msc(iabs(iti),1)*pr(1)*pr(2)
3716 Im(1,3)=Im(1,3)-msc(iabs(iti),1)*pr(1)*pr(3)
3717 Im(2,3)=Im(2,3)-msc(iabs(iti),1)*pr(2)*pr(3)
3718 Im(2,2)=Im(2,2)+msc(iabs(iti),1)*(pr(3)*pr(3)+pr(1)*pr(1))
3719 Im(3,3)=Im(3,3)+msc(iabs(iti),1)*(pr(1)*pr(1)+pr(2)*pr(2))
3723 Im(1,1)=Im(1,1)+Ip(1)*(1-dc_norm(1,i)*dc_norm(1,i))* &
3724 vbld(i+1)*vbld(i+1)*0.25d0
3725 Im(1,2)=Im(1,2)+Ip(1)*(-dc_norm(1,i)*dc_norm(2,i))* &
3726 vbld(i+1)*vbld(i+1)*0.25d0
3727 Im(1,3)=Im(1,3)+Ip(1)*(-dc_norm(1,i)*dc_norm(3,i))* &
3728 vbld(i+1)*vbld(i+1)*0.25d0
3729 Im(2,3)=Im(2,3)+Ip(1)*(-dc_norm(2,i)*dc_norm(3,i))* &
3730 vbld(i+1)*vbld(i+1)*0.25d0
3731 Im(2,2)=Im(2,2)+Ip(1)*(1-dc_norm(2,i)*dc_norm(2,i))* &
3732 vbld(i+1)*vbld(i+1)*0.25d0
3733 Im(3,3)=Im(3,3)+Ip(1)*(1-dc_norm(3,i)*dc_norm(3,i))* &
3734 vbld(i+1)*vbld(i+1)*0.25d0
3739 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3740 iti=iabs(itype(i,1))
3742 Im(1,1)=Im(1,1)+Isc(iti,1)*(1-dc_norm(1,inres)* &
3743 dc_norm(1,inres))*vbld(inres)*vbld(inres)
3744 Im(1,2)=Im(1,2)-Isc(iti,1)*(dc_norm(1,inres)* &
3745 dc_norm(2,inres))*vbld(inres)*vbld(inres)
3746 Im(1,3)=Im(1,3)-Isc(iti,1)*(dc_norm(1,inres)* &
3747 dc_norm(3,inres))*vbld(inres)*vbld(inres)
3748 Im(2,3)=Im(2,3)-Isc(iti,1)*(dc_norm(2,inres)* &
3749 dc_norm(3,inres))*vbld(inres)*vbld(inres)
3750 Im(2,2)=Im(2,2)+Isc(iti,1)*(1-dc_norm(2,inres)* &
3751 dc_norm(2,inres))*vbld(inres)*vbld(inres)
3752 Im(3,3)=Im(3,3)+Isc(iti,1)*(1-dc_norm(3,inres)* &
3753 dc_norm(3,inres))*vbld(inres)*vbld(inres)
3758 ! write(iout,*) "The angular momentum before adjustment:"
3759 ! write(iout,*) (L(j),j=1,3)
3765 ! Copying the Im matrix for the djacob subroutine
3773 ! Finding the eigenvectors and eignvalues of the inertia tensor
3774 call djacob(3,3,10000,1.0d-10,Imcp,eigvec,eigval)
3775 ! write (iout,*) "Eigenvalues & Eigenvectors"
3776 ! write (iout,'(5x,3f10.5)') (eigval(i),i=1,3)
3779 ! write (iout,'(i5,3f10.5)') i,(eigvec(i,j),j=1,3)
3781 ! Constructing the diagonalized matrix
3783 if (dabs(eigval(i)).gt.1.0d-15) then
3784 Id(i,i)=1.0d0/eigval(i)
3791 Imcp(i,j)=eigvec(j,i)
3797 pr1(i,j)=pr1(i,j)+Id(i,k)*Imcp(k,j)
3804 pr2(i,j)=pr2(i,j)+eigvec(i,k)*pr1(k,j)
3808 ! Calculating the total rotational velocity of the molecule
3811 vrot(i)=vrot(i)+pr2(i,j)*L(j)
3814 ! Resetting the velocities
3816 call vecpr(vrot(1),dc(1,i),vp)
3818 d_t(j,i)=d_t(j,i)-vp(j)
3822 if(itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3824 call vecpr(vrot(1),dc(1,inres),vp)
3826 d_t(j,inres)=d_t(j,inres)-vp(j)
3831 ! write(iout,*) "The angular momentum after adjustment:"
3832 ! write(iout,*) (L(j),j=1,3)
3835 end subroutine inertia_tensor
3836 !-----------------------------------------------------------------------------
3837 subroutine angmom(cm,L)
3840 ! implicit real*8 (a-h,o-z)
3841 ! include 'DIMENSIONS'
3842 ! include 'COMMON.CONTROL'
3843 ! include 'COMMON.VAR'
3844 ! include 'COMMON.MD'
3845 ! include 'COMMON.CHAIN'
3846 ! include 'COMMON.DERIV'
3847 ! include 'COMMON.GEO'
3848 ! include 'COMMON.LOCAL'
3849 ! include 'COMMON.INTERACT'
3850 ! include 'COMMON.IOUNITS'
3851 ! include 'COMMON.NAMES'
3853 real(kind=8),dimension(3) :: L,cm,pr,vp,vrot,incr,v,pp
3854 integer :: iti,inres,i,j
3855 ! Calculate the angular momentum
3864 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
3867 v(j)=incr(j)+0.5d0*d_t(j,i)
3870 incr(j)=incr(j)+d_t(j,i)
3872 call vecpr(pr(1),v(1),vp)
3874 L(j)=L(j)+mp(1)*vp(j)
3878 pp(j)=0.5d0*d_t(j,i)
3880 call vecpr(pr(1),pp(1),vp)
3882 L(j)=L(j)+Ip(1)*vp(j)
3889 iti=iabs(itype(i,1))
3892 pr(j)=c(j,inres)-cm(j)
3894 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3896 v(j)=incr(j)+d_t(j,inres)
3903 call vecpr(pr(1),v(1),vp)
3904 ! write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),
3905 ! & " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
3907 L(j)=L(j)+msc(iabs(iti),1)*vp(j)
3909 ! write (iout,*) "L",(l(j),j=1,3)
3910 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3912 v(j)=incr(j)+d_t(j,inres)
3914 call vecpr(dc(1,inres),d_t(1,inres),vp)
3916 L(j)=L(j)+Isc(iti,1)*vp(j)
3920 incr(j)=incr(j)+d_t(j,i)
3924 end subroutine angmom
3925 !-----------------------------------------------------------------------------
3926 subroutine vcm_vel(vcm)
3929 ! implicit real*8 (a-h,o-z)
3930 ! include 'DIMENSIONS'
3931 ! include 'COMMON.VAR'
3932 ! include 'COMMON.MD'
3933 ! include 'COMMON.CHAIN'
3934 ! include 'COMMON.DERIV'
3935 ! include 'COMMON.GEO'
3936 ! include 'COMMON.LOCAL'
3937 ! include 'COMMON.INTERACT'
3938 ! include 'COMMON.IOUNITS'
3939 real(kind=8),dimension(3) :: vcm,vv
3940 real(kind=8) :: summas,amas
3952 vcm(j)=vcm(j)+mp(1)*(vv(j)+0.5d0*d_t(j,i))
3955 amas=msc(iabs(itype(i,1)),1)
3957 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3959 vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres))
3963 vcm(j)=vcm(j)+amas*vv(j)
3967 vv(j)=vv(j)+d_t(j,i)
3970 ! write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas
3972 vcm(j)=vcm(j)/summas
3975 end subroutine vcm_vel
3976 !-----------------------------------------------------------------------------
3978 !-----------------------------------------------------------------------------
3980 ! RATTLE algorithm for velocity Verlet - step 1, UNRES
3984 ! implicit real*8 (a-h,o-z)
3985 ! include 'DIMENSIONS'
3987 ! include 'COMMON.CONTROL'
3988 ! include 'COMMON.VAR'
3989 ! include 'COMMON.MD'
3991 ! include 'COMMON.LANGEVIN'
3993 ! include 'COMMON.LANGEVIN.lang0'
3995 ! include 'COMMON.CHAIN'
3996 ! include 'COMMON.DERIV'
3997 ! include 'COMMON.GEO'
3998 ! include 'COMMON.LOCAL'
3999 ! include 'COMMON.INTERACT'
4000 ! include 'COMMON.IOUNITS'
4001 ! include 'COMMON.NAMES'
4002 ! include 'COMMON.TIME1'
4003 !el real(kind=8) :: gginv(2*nres,2*nres),&
4004 !el gdc(3,2*nres,2*nres)
4005 real(kind=8) :: dC_uncor(3,2*nres) !,&
4006 !el real(kind=8) :: Cmat(2*nres,2*nres)
4007 real(kind=8) :: x(2*nres),xcorr(3,2*nres) !maxres2=2*maxres
4008 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
4009 !el common /przechowalnia/ nbond
4010 integer :: max_rattle = 5
4011 logical :: lprn = .false., lprn1 = .false., not_done
4012 real(kind=8) :: tol_rattle = 1.0d-5
4014 integer :: ii,i,j,jj,l,ind,ind1,nres2
4017 !el /common/ przechowalnia
4019 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4020 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4021 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4023 if (lprn) write (iout,*) "RATTLE1"
4026 if (itype(i,1).ne.10) nbond=nbond+1
4028 ! Make a folded form of the Ginv-matrix
4041 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4045 if (itype(k,1).ne.10) then
4049 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4056 if (itype(i,1).ne.10) then
4066 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4070 if (itype(k,1).ne.10) then
4074 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4082 write (iout,*) "Matrix GGinv"
4083 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
4089 if (iter.gt.max_rattle) then
4090 write (iout,*) "Error - too many iterations in RATTLE."
4093 ! Calculate the matrix C = GG**(-1) dC_old o dC
4098 dC_uncor(j,ind1)=dC(j,i)
4102 if (itype(i,1).ne.10) then
4105 dC_uncor(j,ind1)=dC(j,i+nres)
4114 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
4118 if (itype(k,1).ne.10) then
4121 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
4126 ! Calculate deviations from standard virtual-bond lengths
4130 x(ind)=vbld(i+1)**2-vbl**2
4133 if (itype(i,1).ne.10) then
4135 x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4139 write (iout,*) "Coordinates and violations"
4141 write(iout,'(i5,3f10.5,5x,e15.5)') &
4142 i,(dC_uncor(j,i),j=1,3),x(i)
4144 write (iout,*) "Velocities and violations"
4148 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4149 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4152 if (itype(i,1).ne.10) then
4154 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4155 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4156 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4159 ! write (iout,*) "gdc"
4161 ! write (iout,*) "i",i
4163 ! write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4169 if (dabs(x(i)).gt.xmax) then
4173 if (xmax.lt.tol_rattle) then
4177 ! Calculate the matrix of the system of equations
4182 Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
4187 write (iout,*) "Matrix Cmat"
4188 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4190 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4191 ! Add constraint term to positions
4198 xx = xx+x(ii)*gdc(j,ind,ii)
4202 d_t_new(j,i)=d_t_new(j,i)-xx/d_time
4206 if (itype(i,1).ne.10) then
4211 xx = xx+x(ii)*gdc(j,ind,ii)
4214 dC(j,i+nres)=dC(j,i+nres)-xx
4215 d_t_new(j,i+nres)=d_t_new(j,i+nres)-xx/d_time
4219 ! Rebuild the chain using the new coordinates
4220 call chainbuild_cart
4222 write (iout,*) "New coordinates, Lagrange multipliers,",&
4223 " and differences between actual and standard bond lengths"
4227 xx=vbld(i+1)**2-vbl**2
4228 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4229 i,(dC(j,i),j=1,3),x(ind),xx
4232 if (itype(i,1).ne.10) then
4234 xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4235 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4236 i,(dC(j,i+nres),j=1,3),x(ind),xx
4239 write (iout,*) "Velocities and violations"
4243 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4244 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4247 if (itype(i,1).ne.10) then
4249 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4250 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4251 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4258 10 write (iout,*) "Error - singularity in solving the system",&
4259 " of equations for Lagrange multipliers."
4263 "RATTLE inactive; use -DRATTLE switch at compile time."
4266 end subroutine rattle1
4267 !-----------------------------------------------------------------------------
4269 ! RATTLE algorithm for velocity Verlet - step 2, UNRES
4273 ! implicit real*8 (a-h,o-z)
4274 ! include 'DIMENSIONS'
4276 ! include 'COMMON.CONTROL'
4277 ! include 'COMMON.VAR'
4278 ! include 'COMMON.MD'
4280 ! include 'COMMON.LANGEVIN'
4282 ! include 'COMMON.LANGEVIN.lang0'
4284 ! include 'COMMON.CHAIN'
4285 ! include 'COMMON.DERIV'
4286 ! include 'COMMON.GEO'
4287 ! include 'COMMON.LOCAL'
4288 ! include 'COMMON.INTERACT'
4289 ! include 'COMMON.IOUNITS'
4290 ! include 'COMMON.NAMES'
4291 ! include 'COMMON.TIME1'
4292 !el real(kind=8) :: gginv(2*nres,2*nres),&
4293 !el gdc(3,2*nres,2*nres)
4294 real(kind=8) :: dC_uncor(3,2*nres) !,&
4295 !el Cmat(2*nres,2*nres)
4296 real(kind=8) :: x(2*nres) !maxres2=2*maxres
4297 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
4298 !el common /przechowalnia/ nbond
4299 integer :: max_rattle = 5
4300 logical :: lprn = .false., lprn1 = .false., not_done
4301 real(kind=8) :: tol_rattle = 1.0d-5
4305 !el /common/ przechowalnia
4306 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4307 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4308 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4310 if (lprn) write (iout,*) "RATTLE2"
4311 if (lprn) write (iout,*) "Velocity correction"
4312 ! Calculate the matrix G dC
4318 gdc(j,i,ind)=GGinv(i,ind)*dC(j,k)
4322 if (itype(k,1).ne.10) then
4325 gdc(j,i,ind)=GGinv(i,ind)*dC(j,k+nres)
4331 ! write (iout,*) "gdc"
4333 ! write (iout,*) "i",i
4335 ! write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4339 ! Calculate the matrix of the system of equations
4346 Cmat(ind,j)=Cmat(ind,j)+dC(k,i)*gdc(k,ind,j)
4351 if (itype(i,1).ne.10) then
4356 Cmat(ind,j)=Cmat(ind,j)+dC(k,i+nres)*gdc(k,ind,j)
4361 ! Calculate the scalar product dC o d_t_new
4365 x(ind)=scalar(d_t(1,i),dC(1,i))
4368 if (itype(i,1).ne.10) then
4370 x(ind)=scalar(d_t(1,i+nres),dC(1,i+nres))
4374 write (iout,*) "Velocities and violations"
4378 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4379 i,ind,(d_t(j,i),j=1,3),x(ind)
4382 if (itype(i,1).ne.10) then
4384 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4385 i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind)
4391 if (dabs(x(i)).gt.xmax) then
4395 if (xmax.lt.tol_rattle) then
4400 write (iout,*) "Matrix Cmat"
4401 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4403 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4404 ! Add constraint term to velocities
4411 xx = xx+x(ii)*gdc(j,ind,ii)
4413 d_t(j,i)=d_t(j,i)-xx
4417 if (itype(i,1).ne.10) then
4422 xx = xx+x(ii)*gdc(j,ind,ii)
4424 d_t(j,i+nres)=d_t(j,i+nres)-xx
4430 "New velocities, Lagrange multipliers violations"
4434 if (lprn) write (iout,'(2i5,3f10.5,5x,2e15.5)') &
4435 i,ind,(d_t(j,i),j=1,3),x(ind),scalar(d_t(1,i),dC(1,i))
4438 if (itype(i,1).ne.10) then
4440 write (iout,'(2i5,3f10.5,5x,2e15.5)') &
4441 i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind),&
4442 scalar(d_t(1,i+nres),dC(1,i+nres))
4448 10 write (iout,*) "Error - singularity in solving the system",&
4449 " of equations for Lagrange multipliers."
4453 "RATTLE inactive; use -DRATTLE option at compile time."
4456 end subroutine rattle2
4457 !-----------------------------------------------------------------------------
4458 subroutine rattle_brown
4459 ! RATTLE/LINCS algorithm for Brownian dynamics, UNRES
4463 ! implicit real*8 (a-h,o-z)
4464 ! include 'DIMENSIONS'
4466 ! include 'COMMON.CONTROL'
4467 ! include 'COMMON.VAR'
4468 ! include 'COMMON.MD'
4470 ! include 'COMMON.LANGEVIN'
4472 ! include 'COMMON.LANGEVIN.lang0'
4474 ! include 'COMMON.CHAIN'
4475 ! include 'COMMON.DERIV'
4476 ! include 'COMMON.GEO'
4477 ! include 'COMMON.LOCAL'
4478 ! include 'COMMON.INTERACT'
4479 ! include 'COMMON.IOUNITS'
4480 ! include 'COMMON.NAMES'
4481 ! include 'COMMON.TIME1'
4482 !el real(kind=8) :: gginv(2*nres,2*nres),&
4483 !el gdc(3,2*nres,2*nres)
4484 real(kind=8) :: dC_uncor(3,2*nres) !,&
4485 !el real(kind=8) :: Cmat(2*nres,2*nres)
4486 real(kind=8) :: x(2*nres) !maxres2=2*maxres
4487 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
4488 !el common /przechowalnia/ nbond
4489 integer :: max_rattle = 5
4490 logical :: lprn = .false., lprn1 = .false., not_done
4491 real(kind=8) :: tol_rattle = 1.0d-5
4495 !el /common/ przechowalnia
4496 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4497 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4498 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4501 if (lprn) write (iout,*) "RATTLE_BROWN"
4504 if (itype(i,1).ne.10) nbond=nbond+1
4506 ! Make a folded form of the Ginv-matrix
4519 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4523 if (itype(k,1).ne.10) then
4527 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4534 if (itype(i,1).ne.10) then
4544 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4548 if (itype(k,1).ne.10) then
4552 if (j.eq.1 .and. l.eq.1)GGinv(ii,jj)=fricmat(ind,ind1)
4560 write (iout,*) "Matrix GGinv"
4561 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
4567 if (iter.gt.max_rattle) then
4568 write (iout,*) "Error - too many iterations in RATTLE."
4571 ! Calculate the matrix C = GG**(-1) dC_old o dC
4576 dC_uncor(j,ind1)=dC(j,i)
4580 if (itype(i,1).ne.10) then
4583 dC_uncor(j,ind1)=dC(j,i+nres)
4592 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
4596 if (itype(k,1).ne.10) then
4599 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
4604 ! Calculate deviations from standard virtual-bond lengths
4608 x(ind)=vbld(i+1)**2-vbl**2
4611 if (itype(i,1).ne.10) then
4613 x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4617 write (iout,*) "Coordinates and violations"
4619 write(iout,'(i5,3f10.5,5x,e15.5)') &
4620 i,(dC_uncor(j,i),j=1,3),x(i)
4622 write (iout,*) "Velocities and violations"
4626 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4627 i,ind,(d_t(j,i),j=1,3),scalar(d_t(1,i),dC_old(1,i))
4630 if (itype(i,1).ne.10) then
4632 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4633 i+nres,ind,(d_t(j,i+nres),j=1,3),&
4634 scalar(d_t(1,i+nres),dC_old(1,i+nres))
4637 write (iout,*) "gdc"
4639 write (iout,*) "i",i
4641 write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4647 if (dabs(x(i)).gt.xmax) then
4651 if (xmax.lt.tol_rattle) then
4655 ! Calculate the matrix of the system of equations
4660 Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
4665 write (iout,*) "Matrix Cmat"
4666 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4668 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4669 ! Add constraint term to positions
4676 xx = xx+x(ii)*gdc(j,ind,ii)
4679 d_t(j,i)=d_t(j,i)+xx/d_time
4684 if (itype(i,1).ne.10) then
4689 xx = xx+x(ii)*gdc(j,ind,ii)
4692 d_t(j,i+nres)=d_t(j,i+nres)+xx/d_time
4693 dC(j,i+nres)=dC(j,i+nres)+xx
4697 ! Rebuild the chain using the new coordinates
4698 call chainbuild_cart
4700 write (iout,*) "New coordinates, Lagrange multipliers,",&
4701 " and differences between actual and standard bond lengths"
4705 xx=vbld(i+1)**2-vbl**2
4706 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4707 i,(dC(j,i),j=1,3),x(ind),xx
4710 if (itype(i,1).ne.10) then
4712 xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4713 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4714 i,(dC(j,i+nres),j=1,3),x(ind),xx
4717 write (iout,*) "Velocities and violations"
4721 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4722 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4725 if (itype(i,1).ne.10) then
4727 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4728 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4729 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4736 10 write (iout,*) "Error - singularity in solving the system",&
4737 " of equations for Lagrange multipliers."
4741 "RATTLE inactive; use -DRATTLE option at compile time"
4744 end subroutine rattle_brown
4745 !-----------------------------------------------------------------------------
4747 !-----------------------------------------------------------------------------
4748 subroutine friction_force
4753 ! implicit real*8 (a-h,o-z)
4754 ! include 'DIMENSIONS'
4755 ! include 'COMMON.VAR'
4756 ! include 'COMMON.CHAIN'
4757 ! include 'COMMON.DERIV'
4758 ! include 'COMMON.GEO'
4759 ! include 'COMMON.LOCAL'
4760 ! include 'COMMON.INTERACT'
4761 ! include 'COMMON.MD'
4763 ! include 'COMMON.LANGEVIN'
4765 ! include 'COMMON.LANGEVIN.lang0'
4767 ! include 'COMMON.IOUNITS'
4768 !el real(kind=8),dimension(6*nres) :: gamvec !(MAXRES6) maxres6=6*maxres
4769 !el common /syfek/ gamvec
4770 real(kind=8) :: vv(3),vvtot(3,nres),v_work(6*nres) !,&
4771 !el ginvfric(2*nres,2*nres) !maxres2=2*maxres
4772 !el common /przechowalnia/ ginvfric
4774 logical :: lprn = .false., checkmode = .false.
4775 integer :: i,j,ind,k,nres2,nres6
4779 if(.not.allocated(gamvec)) allocate(gamvec(nres6)) !(MAXRES6)
4780 if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
4788 d_t_work(j)=d_t(j,0)
4793 d_t_work(ind+j)=d_t(j,i)
4798 if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
4800 d_t_work(ind+j)=d_t(j,i+nres)
4806 call fricmat_mult(d_t_work,fric_work)
4808 if (.not.checkmode) return
4811 write (iout,*) "d_t_work and fric_work"
4813 write (iout,'(i3,2e15.5)') i,d_t_work(i),fric_work(i)
4817 friction(j,0)=fric_work(j)
4822 friction(j,i)=fric_work(ind+j)
4827 if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
4829 friction(j,i+nres)=fric_work(ind+j)
4835 write(iout,*) "Friction backbone"
4837 write(iout,'(i5,3e15.5,5x,3e15.5)') &
4838 i,(friction(j,i),j=1,3),(d_t(j,i),j=1,3)
4840 write(iout,*) "Friction side chain"
4842 write(iout,'(i5,3e15.5,5x,3e15.5)') &
4843 i,(friction(j,i+nres),j=1,3),(d_t(j,i+nres),j=1,3)
4852 vvtot(j,i)=vv(j)+0.5d0*d_t(j,i)
4853 vvtot(j,i+nres)=vv(j)+d_t(j,i+nres)
4854 vv(j)=vv(j)+d_t(j,i)
4857 write (iout,*) "vvtot backbone and sidechain"
4859 write (iout,'(i5,3e15.5,5x,3e15.5)') i,(vvtot(j,i),j=1,3),&
4860 (vvtot(j,i+nres),j=1,3)
4865 v_work(ind+j)=vvtot(j,i)
4871 v_work(ind+j)=vvtot(j,i+nres)
4875 write (iout,*) "v_work gamvec and site-based friction forces"
4877 write (iout,'(i5,3e15.5)') i,v_work(i),gamvec(i),&
4881 ! fric_work1(i)=0.0d0
4883 ! fric_work1(i)=fric_work1(i)-A(j,i)*gamvec(j)*v_work(j)
4886 ! write (iout,*) "fric_work and fric_work1"
4888 ! write (iout,'(i5,2e15.5)') i,fric_work(i),fric_work1(i)
4894 ginvfric(i,j)=ginvfric(i,j)+ginv(i,k)*fricmat(k,j)
4898 write (iout,*) "ginvfric"
4900 write (iout,'(i5,100f8.3)') i,(ginvfric(i,j),j=1,dimen)
4902 write (iout,*) "symmetry check"
4905 write (iout,*) i,j,ginvfric(i,j)-ginvfric(j,i)
4910 end subroutine friction_force
4911 !-----------------------------------------------------------------------------
4912 subroutine setup_fricmat
4916 use control_data, only:time_Bcast
4917 use control, only:tcpu
4919 ! implicit real*8 (a-h,o-z)
4923 real(kind=8) :: time00
4925 ! include 'DIMENSIONS'
4926 ! include 'COMMON.VAR'
4927 ! include 'COMMON.CHAIN'
4928 ! include 'COMMON.DERIV'
4929 ! include 'COMMON.GEO'
4930 ! include 'COMMON.LOCAL'
4931 ! include 'COMMON.INTERACT'
4932 ! include 'COMMON.MD'
4933 ! include 'COMMON.SETUP'
4934 ! include 'COMMON.TIME1'
4935 ! integer licznik /0/
4938 ! include 'COMMON.LANGEVIN'
4940 ! include 'COMMON.LANGEVIN.lang0'
4942 ! include 'COMMON.IOUNITS'
4944 integer :: i,j,ind,ind1,m
4945 logical :: lprn = .false.
4946 real(kind=8) :: dtdi !el ,gamvec(2*nres)
4947 !el real(kind=8),dimension(2*nres,2*nres) :: ginvfric,fcopy
4948 real(kind=8),dimension(2*nres,2*nres) :: fcopy
4949 !el real(kind=8),dimension(2*nres*(2*nres+1)/2) :: Ghalf !(mmaxres2) (mmaxres2=(maxres2*(maxres2+1)/2))
4950 !el common /syfek/ gamvec
4951 real(kind=8) :: work(8*2*nres)
4952 integer :: iwork(2*nres)
4953 !el common /przechowalnia/ ginvfric,Ghalf,fcopy
4954 integer :: ii,iti,k,l,nzero,nres2,nres6,ierr
4956 if (fg_rank.ne.king) goto 10
4961 if(.not.allocated(gamvec)) allocate(gamvec(nres2)) !(MAXRES2)
4962 if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
4963 !el if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
4964 !el allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
4965 if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !maxres2=2*maxres
4967 !el if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
4968 ! Zeroing out fricmat
4974 ! Load the friction coefficients corresponding to peptide groups
4980 ! Load the friction coefficients corresponding to side chains
4988 gamvec(ii)=gamsc(iabs(iti))
4990 if (surfarea) call sdarea(gamvec)
4992 ! write (iout,*) "Matrix A and vector gamma"
4994 ! write (iout,'(i2,$)') i
4996 ! write (iout,'(f4.1,$)') A(i,j)
4998 ! write (iout,'(f8.3)') gamvec(i)
5002 write (iout,*) "Vector gamvec"
5004 write (iout,'(i5,f10.5)') i, gamvec(i)
5008 ! The friction matrix
5013 dtdi=dtdi+A(j,k)*A(j,i)*gamvec(j)
5020 write (iout,'(//a)') "Matrix fricmat"
5021 call matout2(dimen,dimen,nres2,nres2,fricmat)
5023 if (lang.eq.2 .or. lang.eq.3) then
5024 ! Mass-scale the friction matrix if non-direct integration will be performed
5030 Ginvfric(i,j)=Ginvfric(i,j)+ &
5031 Gsqrm(i,k)*Gsqrm(l,j)*fricmat(k,l)
5036 ! Diagonalize the friction matrix
5041 Ghalf(ind)=Ginvfric(i,j)
5044 call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
5047 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
5048 " mass-scaled friction matrix"
5049 call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
5051 ! Precompute matrices for tinker stochastic integrator
5058 mt1(i,j)=mt1(i,j)+fricvec(k,i)*gsqrm(k,j)
5059 mt2(i,j)=mt2(i,j)+fricvec(k,i)*gsqrp(k,j)
5065 else if (lang.eq.4) then
5066 ! Diagonalize the friction matrix
5071 Ghalf(ind)=fricmat(i,j)
5074 call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
5077 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
5079 call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
5081 ! Determine the number of zero eigenvalues of the friction matrix
5082 nzero=max0(dimen-dimen1,0)
5083 ! do while (fricgam(nzero+1).le.1.0d-5 .and. nzero.lt.dimen)
5086 write (iout,*) "Number of zero eigenvalues:",nzero
5091 fricmat(i,j)=fricmat(i,j) &
5092 +fricvec(i,k)*fricvec(j,k)/fricgam(k)
5097 write (iout,'(//a)') "Generalized inverse of fricmat"
5098 call matout(dimen,dimen,nres6,nres6,fricmat)
5103 if (nfgtasks.gt.1) then
5104 if (fg_rank.eq.0) then
5105 ! The matching BROADCAST for fg processors is called in ERGASTULUM
5111 call MPI_Bcast(10,1,MPI_INTEGER,king,FG_COMM,IERROR)
5113 time_Bcast=time_Bcast+MPI_Wtime()-time00
5115 time_Bcast=time_Bcast+tcpu()-time00
5117 ! print *,"Processor",myrank,
5118 ! & " BROADCAST iorder in SETUP_FRICMAT"
5121 write (iout,*) "setup_fricmat licznik"!,licznik !sp
5127 ! Scatter the friction matrix
5128 call MPI_Scatterv(fricmat(1,1),nginv_counts(0),&
5129 nginv_start(0),MPI_DOUBLE_PRECISION,fcopy(1,1),&
5130 myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
5133 time_scatter=time_scatter+MPI_Wtime()-time00
5134 time_scatter_fmat=time_scatter_fmat+MPI_Wtime()-time00
5136 time_scatter=time_scatter+tcpu()-time00
5137 time_scatter_fmat=time_scatter_fmat+tcpu()-time00
5141 do j=1,2*my_ng_count
5142 fricmat(j,i)=fcopy(i,j)
5145 ! write (iout,*) "My chunk of fricmat"
5146 ! call MATOUT2(my_ng_count,dimen,maxres2,maxres2,fcopy)
5150 end subroutine setup_fricmat
5151 !-----------------------------------------------------------------------------
5152 subroutine stochastic_force(stochforcvec)
5155 use random, only:anorm_distr
5156 ! implicit real*8 (a-h,o-z)
5157 ! include 'DIMENSIONS'
5158 use control, only: tcpu
5163 ! include 'COMMON.VAR'
5164 ! include 'COMMON.CHAIN'
5165 ! include 'COMMON.DERIV'
5166 ! include 'COMMON.GEO'
5167 ! include 'COMMON.LOCAL'
5168 ! include 'COMMON.INTERACT'
5169 ! include 'COMMON.MD'
5170 ! include 'COMMON.TIME1'
5172 ! include 'COMMON.LANGEVIN'
5174 ! include 'COMMON.LANGEVIN.lang0'
5176 ! include 'COMMON.IOUNITS'
5178 real(kind=8) :: x,sig,lowb,highb
5179 real(kind=8) :: ff(3),force(3,0:2*nres),zeta2,lowb2
5180 real(kind=8) :: highb2,sig2,forcvec(6*nres),stochforcvec(6*nres)
5181 real(kind=8) :: time00
5182 logical :: lprn = .false.
5187 stochforc(j,i)=0.0d0
5197 ! Compute the stochastic forces acting on bodies. Store in force.
5203 force(j,i)=anorm_distr(x,sig,lowb,highb)
5211 force(j,i+nres)=anorm_distr(x,sig2,lowb2,highb2)
5215 time_fsample=time_fsample+MPI_Wtime()-time00
5217 time_fsample=time_fsample+tcpu()-time00
5219 ! Compute the stochastic forces acting on virtual-bond vectors.
5225 stochforc(j,i)=ff(j)+0.5d0*force(j,i)
5228 ff(j)=ff(j)+force(j,i)
5230 if (itype(i+1,1).ne.ntyp1) then
5232 stochforc(j,i)=stochforc(j,i)+force(j,i+nres+1)
5233 ff(j)=ff(j)+force(j,i+nres+1)
5238 stochforc(j,0)=ff(j)+force(j,nnt+nres)
5241 if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
5243 stochforc(j,i+nres)=force(j,i+nres)
5249 stochforcvec(j)=stochforc(j,0)
5254 stochforcvec(ind+j)=stochforc(j,i)
5259 if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
5261 stochforcvec(ind+j)=stochforc(j,i+nres)
5267 write (iout,*) "stochforcvec"
5269 write(iout,'(i5,e15.5)') i,stochforcvec(i)
5271 write(iout,*) "Stochastic forces backbone"
5273 write(iout,'(i5,3e15.5)') i,(stochforc(j,i),j=1,3)
5275 write(iout,*) "Stochastic forces side chain"
5277 write(iout,'(i5,3e15.5)') &
5278 i,(stochforc(j,i+nres),j=1,3)
5286 write (iout,*) i,ind
5288 forcvec(ind+j)=force(j,i)
5293 write (iout,*) i,ind
5295 forcvec(j+ind)=force(j,i+nres)
5300 write (iout,*) "forcvec"
5304 write (iout,'(2i3,2f10.5)') i,j,force(j,i),&
5311 write (iout,'(2i3,2f10.5)') i,j,force(j,i+nres),&
5320 end subroutine stochastic_force
5321 !-----------------------------------------------------------------------------
5322 subroutine sdarea(gamvec)
5324 ! Scale the friction coefficients according to solvent accessible surface areas
5325 ! Code adapted from TINKER
5329 ! implicit real*8 (a-h,o-z)
5330 ! include 'DIMENSIONS'
5331 ! include 'COMMON.CONTROL'
5332 ! include 'COMMON.VAR'
5333 ! include 'COMMON.MD'
5335 ! include 'COMMON.LANGEVIN'
5337 ! include 'COMMON.LANGEVIN.lang0'
5339 ! include 'COMMON.CHAIN'
5340 ! include 'COMMON.DERIV'
5341 ! include 'COMMON.GEO'
5342 ! include 'COMMON.LOCAL'
5343 ! include 'COMMON.INTERACT'
5344 ! include 'COMMON.IOUNITS'
5345 ! include 'COMMON.NAMES'
5346 real(kind=8),dimension(2*nres) :: radius,gamvec !(maxres2)
5347 real(kind=8),parameter :: twosix = 1.122462048309372981d0
5348 logical :: lprn = .false.
5349 real(kind=8) :: probe,area,ratio
5350 integer :: i,j,ind,iti
5352 ! determine new friction coefficients every few SD steps
5354 ! set the atomic radii to estimates of sigma values
5356 ! print *,"Entered sdarea"
5362 ! Load peptide group radii
5366 ! Load side chain radii
5369 radius(i+nres)=restok(iti,1)
5372 ! write (iout,*) "i",i," radius",radius(i)
5375 radius(i) = radius(i) / twosix
5376 if (radius(i) .ne. 0.0d0) radius(i) = radius(i) + probe
5379 ! scale atomic friction coefficients by accessible area
5381 if (lprn) write (iout,*) &
5382 "Original gammas, surface areas, scaling factors, new gammas, ",&
5383 "std's of stochastic forces"
5386 if (radius(i).gt.0.0d0) then
5387 call surfatom (i,area,radius)
5388 ratio = dmax1(area/(4.0d0*pi*radius(i)**2),1.0d-1)
5389 if (lprn) write (iout,'(i5,3f10.5,$)') &
5390 i,gamvec(ind+1),area,ratio
5393 gamvec(ind) = ratio * gamvec(ind)
5395 stdforcp(i)=stdfp*dsqrt(gamvec(ind))
5396 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcp(i)
5400 if (radius(i+nres).gt.0.0d0) then
5401 call surfatom (i+nres,area,radius)
5402 ratio = dmax1(area/(4.0d0*pi*radius(i+nres)**2),1.0d-1)
5403 if (lprn) write (iout,'(i5,3f10.5,$)') &
5404 i,gamvec(ind+1),area,ratio
5407 gamvec(ind) = ratio * gamvec(ind)
5409 stdforcsc(i)=stdfsc(itype(i,1))*dsqrt(gamvec(ind))
5410 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcsc(i)
5415 end subroutine sdarea
5416 !-----------------------------------------------------------------------------
5418 !-----------------------------------------------------------------------------
5421 ! ###################################################
5422 ! ## COPYRIGHT (C) 1996 by Jay William Ponder ##
5423 ! ## All Rights Reserved ##
5424 ! ###################################################
5426 ! ################################################################
5428 ! ## subroutine surfatom -- exposed surface area of an atom ##
5430 ! ################################################################
5433 ! "surfatom" performs an analytical computation of the surface
5434 ! area of a specified atom; a simplified version of "surface"
5436 ! literature references:
5438 ! T. J. Richmond, "Solvent Accessible Surface Area and
5439 ! Excluded Volume in Proteins", Journal of Molecular Biology,
5442 ! L. Wesson and D. Eisenberg, "Atomic Solvation Parameters
5443 ! Applied to Molecular Dynamics of Proteins in Solution",
5444 ! Protein Science, 1, 227-235 (1992)
5446 ! variables and parameters:
5448 ! ir number of atom for which area is desired
5449 ! area accessible surface area of the atom
5450 ! radius radii of each of the individual atoms
5453 subroutine surfatom(ir,area,radius)
5455 ! implicit real*8 (a-h,o-z)
5456 ! include 'DIMENSIONS'
5458 ! include 'COMMON.GEO'
5459 ! include 'COMMON.IOUNITS'
5461 integer :: nsup,nstart_sup
5462 ! double precision c,dc,dc_old,d_c_work,xloc,xrot,dc_norm
5463 ! common /chain/ c(3,maxres2+2),dc(3,0:maxres2),dc_old(3,0:maxres2),
5464 ! & xloc(3,maxres),xrot(3,maxres),dc_norm(3,0:maxres2),
5465 ! & dc_work(MAXRES6),nres,nres0
5466 integer,parameter :: maxarc=300
5470 integer :: mi,ni,narc
5471 integer :: key(maxarc)
5472 integer :: intag(maxarc)
5473 integer :: intag1(maxarc)
5474 real(kind=8) :: area,arcsum
5475 real(kind=8) :: arclen,exang
5476 real(kind=8) :: delta,delta2
5477 real(kind=8) :: eps,rmove
5478 real(kind=8) :: xr,yr,zr
5479 real(kind=8) :: rr,rrsq
5480 real(kind=8) :: rplus,rminus
5481 real(kind=8) :: axx,axy,axz
5482 real(kind=8) :: ayx,ayy
5483 real(kind=8) :: azx,azy,azz
5484 real(kind=8) :: uxj,uyj,uzj
5485 real(kind=8) :: tx,ty,tz
5486 real(kind=8) :: txb,tyb,td
5487 real(kind=8) :: tr2,tr,txr,tyr
5488 real(kind=8) :: tk1,tk2
5489 real(kind=8) :: thec,the,t,tb
5490 real(kind=8) :: txk,tyk,tzk
5491 real(kind=8) :: t1,ti,tf,tt
5492 real(kind=8) :: txj,tyj,tzj
5493 real(kind=8) :: ccsq,cc,xysq
5494 real(kind=8) :: bsqk,bk,cosine
5495 real(kind=8) :: dsqj,gi,pix2
5496 real(kind=8) :: therk,dk,gk
5497 real(kind=8) :: risqk,rik
5498 real(kind=8) :: radius(2*nres) !(maxatm) (maxatm=maxres2)
5499 real(kind=8) :: ri(maxarc),risq(maxarc)
5500 real(kind=8) :: ux(maxarc),uy(maxarc),uz(maxarc)
5501 real(kind=8) :: xc(maxarc),yc(maxarc),zc(maxarc)
5502 real(kind=8) :: xc1(maxarc),yc1(maxarc),zc1(maxarc)
5503 real(kind=8) :: dsq(maxarc),bsq(maxarc)
5504 real(kind=8) :: dsq1(maxarc),bsq1(maxarc)
5505 real(kind=8) :: arci(maxarc),arcf(maxarc)
5506 real(kind=8) :: ex(maxarc),lt(maxarc),gr(maxarc)
5507 real(kind=8) :: b(maxarc),b1(maxarc),bg(maxarc)
5508 real(kind=8) :: kent(maxarc),kout(maxarc)
5509 real(kind=8) :: ther(maxarc)
5510 logical :: moved,top
5511 logical :: omit(maxarc)
5514 ! maxatm = 2*nres !maxres2 maxres2=2*maxres
5515 ! maxlight = 8*maxatm
5518 ! maxtors = 4*maxatm
5520 ! zero out the surface area for the sphere of interest
5523 ! write (2,*) "ir",ir," radius",radius(ir)
5524 if (radius(ir) .eq. 0.0d0) return
5526 ! set the overlap significance and connectivity shift
5530 delta2 = delta * delta
5535 ! store coordinates and radius of the sphere of interest
5543 ! initialize values of some counters and summations
5552 ! test each sphere to see if it overlaps the sphere of interest
5555 if (i.eq.ir .or. radius(i).eq.0.0d0) goto 30
5556 rplus = rr + radius(i)
5558 if (abs(tx) .ge. rplus) goto 30
5560 if (abs(ty) .ge. rplus) goto 30
5562 if (abs(tz) .ge. rplus) goto 30
5564 ! check for sphere overlap by testing distance against radii
5566 xysq = tx*tx + ty*ty
5567 if (xysq .lt. delta2) then
5574 if (rplus-cc .le. delta) goto 30
5575 rminus = rr - radius(i)
5577 ! check to see if sphere of interest is completely buried
5579 if (cc-abs(rminus) .le. delta) then
5580 if (rminus .le. 0.0d0) goto 170
5584 ! check for too many overlaps with sphere of interest
5586 if (io .ge. maxarc) then
5588 20 format (/,' SURFATOM -- Increase the Value of MAXARC')
5592 ! get overlap between current sphere and sphere of interest
5601 gr(io) = (ccsq+rplus*rminus) / (2.0d0*rr*b1(io))
5607 ! case where no other spheres overlap the sphere of interest
5610 area = 4.0d0 * pi * rrsq
5614 ! case where only one sphere overlaps the sphere of interest
5617 area = pix2 * (1.0d0 + gr(1))
5618 area = mod(area,4.0d0*pi) * rrsq
5622 ! case where many spheres intersect the sphere of interest;
5623 ! sort the intersecting spheres by their degree of overlap
5625 call sort2 (io,gr,key)
5628 intag(i) = intag1(k)
5637 ! get radius of each overlap circle on surface of the sphere
5642 risq(i) = rrsq - gi*gi
5643 ri(i) = sqrt(risq(i))
5644 ther(i) = 0.5d0*pi - asin(min(1.0d0,max(-1.0d0,gr(i))))
5647 ! find boundary of inaccessible area on sphere of interest
5650 if (.not. omit(k)) then
5657 ! check to see if J circle is intersecting K circle;
5658 ! get distance between circle centers and sum of radii
5661 if (omit(j)) goto 60
5662 cc = (txk*xc(j)+tyk*yc(j)+tzk*zc(j))/(bk*b(j))
5663 cc = acos(min(1.0d0,max(-1.0d0,cc)))
5664 td = therk + ther(j)
5666 ! check to see if circles enclose separate regions
5668 if (cc .ge. td) goto 60
5670 ! check for circle J completely inside circle K
5672 if (cc+ther(j) .lt. therk) goto 40
5674 ! check for circles that are essentially parallel
5676 if (cc .gt. delta) goto 50
5681 ! check to see if sphere of interest is completely buried
5684 if (pix2-cc .le. td) goto 170
5690 ! find T value of circle intersections
5693 if (omit(k)) goto 110
5708 ! rotation matrix elements
5720 if (.not. omit(j)) then
5725 ! rotate spheres so K vector colinear with z-axis
5727 uxj = txj*axx + tyj*axy - tzj*axz
5728 uyj = tyj*ayy - txj*ayx
5729 uzj = txj*azx + tyj*azy + tzj*azz
5730 cosine = min(1.0d0,max(-1.0d0,uzj/b(j)))
5731 if (acos(cosine) .lt. therk+ther(j)) then
5732 dsqj = uxj*uxj + uyj*uyj
5737 tr2 = risqk*dsqj - tb*tb
5743 ! get T values of intersection for K circle
5746 tb = min(1.0d0,max(-1.0d0,tb))
5748 if (tyb-txr .lt. 0.0d0) tk1 = pix2 - tk1
5750 tb = min(1.0d0,max(-1.0d0,tb))
5752 if (tyb+txr .lt. 0.0d0) tk2 = pix2 - tk2
5753 thec = (rrsq*uzj-gk*bg(j)) / (rik*ri(j)*b(j))
5754 if (abs(thec) .lt. 1.0d0) then
5756 else if (thec .ge. 1.0d0) then
5758 else if (thec .le. -1.0d0) then
5762 ! see if "tk1" is entry or exit point; check t=0 point;
5763 ! "ti" is exit point, "tf" is entry point
5765 cosine = min(1.0d0,max(-1.0d0, &
5766 (uzj*gk-uxj*rik)/(b(j)*rr)))
5767 if ((acos(cosine)-ther(j))*(tk2-tk1) .le. 0.0d0) then
5775 if (narc .ge. maxarc) then
5777 70 format (/,' SURFATOM -- Increase the Value',&
5781 if (tf .le. ti) then
5802 ! special case; K circle without intersections
5804 if (narc .le. 0) goto 90
5806 ! general case; sum up arclength and set connectivity code
5808 call sort2 (narc,arci,key)
5813 if (narc .gt. 1) then
5816 if (t .lt. arci(j)) then
5817 arcsum = arcsum + arci(j) - t
5818 exang = exang + ex(ni)
5820 if (jb .ge. maxarc) then
5822 80 format (/,' SURFATOM -- Increase the Value',&
5827 kent(jb) = maxarc*i + k
5829 kout(jb) = maxarc*k + i
5838 arcsum = arcsum + pix2 - t
5840 exang = exang + ex(ni)
5843 kent(jb) = maxarc*i + k
5845 kout(jb) = maxarc*k + i
5852 arclen = arclen + gr(k)*arcsum
5855 if (arclen .eq. 0.0d0) goto 170
5856 if (jb .eq. 0) goto 150
5858 ! find number of independent boundaries and check connectivity
5862 if (kout(k) .ne. 0) then
5869 if (m .eq. kent(ii)) then
5872 if (j .eq. jb) goto 150
5884 ! attempt to fix connectivity error by moving atom slightly
5888 140 format (/,' SURFATOM -- Connectivity Error at Atom',i6)
5897 ! compute the exposed surface area for the sphere of interest
5900 area = ib*pix2 + exang + arclen
5901 area = mod(area,4.0d0*pi) * rrsq
5903 ! attempt to fix negative area by moving atom slightly
5905 if (area .lt. 0.0d0) then
5908 160 format (/,' SURFATOM -- Negative Area at Atom',i6)
5919 end subroutine surfatom
5920 !----------------------------------------------------------------
5921 !----------------------------------------------------------------
5922 subroutine alloc_MD_arrays
5923 !EL Allocation of arrays used by MD module
5925 integer :: nres2,nres6
5928 !----------------------
5932 allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
5933 allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
5934 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
5935 allocate(fricvec(nres2,nres2))
5936 allocate(pfric_mat(nres2,nres2),vfric_mat(nres2,nres2))
5937 allocate(afric_mat(nres2,nres2),prand_mat(nres2,nres2))
5938 allocate(vrand_mat1(nres2,nres2),vrand_mat2(nres2,nres2)) !(MAXRES2,MAXRES2)
5939 allocate(pfric0_mat(nres2,nres2,0:maxflag_stoch))
5940 allocate(afric0_mat(nres2,nres2,0:maxflag_stoch))
5941 allocate(vfric0_mat(nres2,nres2,0:maxflag_stoch))
5942 allocate(prand0_mat(nres2,nres2,0:maxflag_stoch))
5943 allocate(vrand0_mat1(nres2,nres2,0:maxflag_stoch))
5944 allocate(vrand0_mat2(nres2,nres2,0:maxflag_stoch)) !(MAXRES2,MAXRES2,0:maxflag_stoch)
5945 allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
5947 allocate(mt1(nres2,nres2),mt2(nres2,nres2),mt3(nres2,nres2)) !(maxres2,maxres2)
5948 !----------------------
5950 ! commom.langevin.lang0
5952 allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
5953 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
5954 allocate(fricvec(nres2,nres2)) !(MAXRES2,MAXRES2)
5955 allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
5956 allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
5959 !el if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2))
5960 !----------------------
5961 ! commom.hairpin in CSA module
5962 !----------------------
5963 ! common.mce in MCM_MD module
5964 !----------------------
5966 ! common /mdgrad/ in module.energy
5967 ! common /back_constr/ in module.energy
5968 ! common /qmeas/ in module.energy
5971 allocate(potEcomp(0:n_ene+4)) !(0:n_ene+4)
5973 allocate(d_t(3,0:nres2),d_a(3,0:nres2),d_t_old(3,0:nres2)) !(3,0:MAXRES2)
5974 allocate(d_a_work(nres6)) !(6*MAXRES)
5976 allocate(DM(nres2),DU1(nres2),DU2(nres2))
5977 allocate(DMorig(nres2),DU1orig(nres2),DU2orig(nres2))
5979 allocate(Gmat(nres2,nres2),A(nres2,nres2))
5980 if(.not.allocated(Ginv)) allocate(Ginv(nres2,nres2)) !in control: ergastulum
5981 allocate(Gsqrp(nres2,nres2),Gsqrm(nres2,nres2),Gvec(nres2,nres2)) !(maxres2,maxres2)
5983 allocate(Geigen(nres2)) !(maxres2)
5984 if(.not.allocated(vtot)) allocate(vtot(nres2)) !(maxres2)
5985 ! common /inertia/ in io_conf: parmread
5986 ! real(kind=8),dimension(:),allocatable :: ISC,msc !(ntyp+1)
5987 ! common /langevin/in io read_MDpar
5988 ! real(kind=8),dimension(:),allocatable :: gamsc !(ntyp1)
5989 ! real(kind=8),dimension(:),allocatable :: stdfsc !(ntyp)
5990 ! in io_conf: parmread
5991 ! real(kind=8),dimension(:),allocatable :: restok !(ntyp+1)
5992 ! common /mdpmpi/ in control: ergastulum
5993 if(.not.allocated(ng_start)) allocate(ng_start(0:nfgtasks-1))
5994 if(.not.allocated(ng_counts)) allocate(ng_counts(0:nfgtasks-1))
5995 if(.not.allocated(nginv_counts)) allocate(nginv_counts(0:nfgtasks-1)) !(0:MaxProcs-1)
5996 if(.not.allocated(nginv_start)) allocate(nginv_start(0:nfgtasks)) !(0:MaxProcs)
5997 !----------------------
5998 ! common.muca in read_muca
5999 ! common /double_muca/
6000 ! real(kind=8) :: elow,ehigh,factor,hbin,factor_min
6001 ! real(kind=8),dimension(:),allocatable :: emuca,nemuca,&
6002 ! nemuca2,hist !(4*maxres)
6003 ! common /integer_muca/
6004 ! integer :: nmuca,imtime,muca_smooth
6006 ! real(kind=8),dimension(:),allocatable :: elowi,ehighi !(maxprocs)
6007 !----------------------
6009 ! common /mdgrad/ in module.energy
6010 ! common /back_constr/ in module.energy
6011 ! common /qmeas/ in module.energy
6015 allocate(d_t_work(nres6),d_t_work_new(nres6),d_af_work(nres6))
6016 allocate(d_as_work(nres6),kinetic_force(nres6)) !(MAXRES6)
6017 allocate(d_t_new(3,0:nres2),d_a_old(3,0:nres2),d_a_short(3,0:nres2)) !,d_a !(3,0:MAXRES2)
6018 allocate(stdforcp(nres),stdforcsc(nres)) !(MAXRES)
6019 !----------------------
6021 allocate(D_ban(nres6)) !(MAXRES6) maxres6=6*maxres
6022 ! common /stochcalc/ stochforcvec
6023 allocate(stochforcvec(nres6)) !(MAXRES6) maxres6=6*maxres
6024 !----------------------
6026 end subroutine alloc_MD_arrays
6027 !-----------------------------------------------------------------------------
6028 !-----------------------------------------------------------------------------