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,mscab
657 integer :: i,j,k,iti,mnum
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
676 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
677 ! write (iout,*) "Kinetic trp:",i,(incr(j),j=1,3)
679 v(j)=incr(j)+0.5d0*d_t(j,i)
681 vtot(i)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
682 KEt_p=KEt_p+mp(mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
684 incr(j)=incr(j)+d_t(j,i)
687 ! write(iout,*) 'KEt_p', KEt_p
688 ! The translational part for the side chain virtual bond
689 ! Only now we can initialize incr with zeros. It must be equal
690 ! to the velocities of the first Calpha.
696 iti=iabs(itype(i,mnum))
702 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
703 if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
704 .or.(mnum.eq.5)) then
710 v(j)=incr(j)+d_t(j,nres+i)
713 ! write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3)
714 ! write (iout,*) "i",i," msc",msc(iti)," v",(v(j),j=1,3)
715 KEt_sc=KEt_sc+mscab*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
716 vtot(i+nres)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
718 incr(j)=incr(j)+d_t(j,i)
722 ! write(iout,*) 'KEt_sc', KEt_sc
723 ! The part due to stretching and rotation of the peptide groups
727 ! write (iout,*) "i",i
728 ! write (iout,*) "i",i," mag1",mag1," mag2",mag2
732 ! write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3)
733 KEr_p=KEr_p+Ip(mnum)*(incr(1)*incr(1)+incr(2)*incr(2) &
737 ! write(iout,*) 'KEr_p', KEr_p
738 ! The rotational part of the side chain virtual bond
742 iti=iabs(itype(i,mnum))
743 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
744 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
745 .and.(mnum.ne.5)) then
747 incr(j)=d_t(j,nres+i)
749 ! write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3)
750 KEr_sc=KEr_sc+Isc(iti,mnum)*(incr(1)*incr(1)+incr(2)*incr(2)+ &
754 ! The total kinetic energy
756 ! write(iout,*) 'KEr_sc', KEr_sc
757 KE_total=0.5d0*(KEt_p+KEt_sc+0.25d0*KEr_p+KEr_sc)
758 ! write (iout,*) "KE_total",KE_total
760 end subroutine kinetic
761 !-----------------------------------------------------------------------------
763 !-----------------------------------------------------------------------------
765 !------------------------------------------------
766 ! The driver for molecular dynamics subroutines
767 !------------------------------------------------
770 use control, only:tcpu,ovrtim
771 ! use io_comm, only:ilen
773 use compare, only:secondary2,hairpin
774 use io, only:cartout,statout
775 ! implicit real*8 (a-h,o-z)
776 ! include 'DIMENSIONS'
779 integer :: IERROR,ERRCODE
781 ! include 'COMMON.SETUP'
782 ! include 'COMMON.CONTROL'
783 ! include 'COMMON.VAR'
784 ! include 'COMMON.MD'
786 ! include 'COMMON.LANGEVIN'
788 ! include 'COMMON.LANGEVIN.lang0'
790 ! include 'COMMON.CHAIN'
791 ! include 'COMMON.DERIV'
792 ! include 'COMMON.GEO'
793 ! include 'COMMON.LOCAL'
794 ! include 'COMMON.INTERACT'
795 ! include 'COMMON.IOUNITS'
796 ! include 'COMMON.NAMES'
797 ! include 'COMMON.TIME1'
798 ! include 'COMMON.HAIRPIN'
799 real(kind=8),dimension(3) :: L,vcm
801 real(kind=8),dimension(6*nres) :: v_work,v_transf !(maxres6) maxres6=6*maxres
803 integer :: rstcount !ilen,
805 character(len=50) :: tytul
806 !el common /gucio/ cm
807 integer :: itime,i,j,nharp
808 integer,dimension(4,nres/3) :: iharp !(4,nres/3)(4,maxres/3)
810 real(kind=8) :: tt0,scalfac
815 if (ilen(tmpdir).gt.0) &
816 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_" &
817 //liczba(:ilen(liczba))//'.rst')
819 if (ilen(tmpdir).gt.0) &
820 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"//'.rst')
827 write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
833 ! Determine the inverse of the inertia matrix.
834 call setup_MD_matrices
838 t_MDsetup = MPI_Wtime()-tt0
840 t_MDsetup = tcpu()-tt0
843 ! Entering the MD loop
849 if (lang.eq.2 .or. lang.eq.3) then
853 call sd_verlet_p_setup
855 call sd_verlet_ciccotti_setup
859 pfric0_mat(i,j,0)=pfric_mat(i,j)
860 afric0_mat(i,j,0)=afric_mat(i,j)
861 vfric0_mat(i,j,0)=vfric_mat(i,j)
862 prand0_mat(i,j,0)=prand_mat(i,j)
863 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
864 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
869 flag_stoch(i)=.false.
873 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
875 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
879 else if (lang.eq.1 .or. lang.eq.4) then
883 t_langsetup=MPI_Wtime()-tt0
886 t_langsetup=tcpu()-tt0
889 do itime=1,n_timestep
891 if (large.and. mod(itime,ntwe).eq.0) &
892 write (iout,*) "itime",itime
894 if (lang.gt.0 .and. surfarea .and. &
895 mod(itime,reset_fricmat).eq.0) then
896 if (lang.eq.2 .or. lang.eq.3) then
900 call sd_verlet_p_setup
902 call sd_verlet_ciccotti_setup
906 pfric0_mat(i,j,0)=pfric_mat(i,j)
907 afric0_mat(i,j,0)=afric_mat(i,j)
908 vfric0_mat(i,j,0)=vfric_mat(i,j)
909 prand0_mat(i,j,0)=prand_mat(i,j)
910 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
911 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
916 flag_stoch(i)=.false.
919 else if (lang.eq.1 .or. lang.eq.4) then
922 write (iout,'(a,i10)') &
923 "Friction matrix reset based on surface area, itime",itime
925 if (reset_vel .and. tbf .and. lang.eq.0 &
926 .and. mod(itime,count_reset_vel).eq.0) then
928 write(iout,'(a,f20.2)') &
929 "Velocities reset to random values, time",totT
932 d_t_old(j,i)=d_t(j,i)
936 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
940 d_t(j,0)=d_t(j,0)-vcm(j)
943 kinetic_T=2.0d0/(dimen3*Rb)*EK
944 scalfac=dsqrt(T_bath/kinetic_T)
945 write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
948 d_t_old(j,i)=scalfac*d_t(j,i)
954 ! Time-reversible RESPA algorithm
955 ! (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
956 call RESPA_step(itime)
958 ! Variable time step algorithm.
959 call velverlet_step(itime)
963 call brown_step(itime)
965 print *,"Brown dynamics not here!"
967 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
973 if (mod(itime,ntwe).eq.0) call statout(itime)
987 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1.and.mnum.ne.5) then
990 v_work(ind)=d_t(j,i+nres)
995 write (66,'(80f10.5)') &
996 ((d_t(j,i),j=1,3),i=0,nres-1),((d_t(j,i+nres),j=1,3),i=1,nres)
1000 v_transf(i)=v_transf(i)+gvec(j,i)*v_work(j)
1002 v_transf(i)= v_transf(i)*dsqrt(geigen(i))
1004 write (67,'(80f10.5)') (v_transf(i),i=1,ind)
1007 if (mod(itime,ntwx).eq.0) then
1008 write (tytul,'("time",f8.2)') totT
1010 call hairpin(.true.,nharp,iharp)
1011 call secondary2(.true.)
1012 call pdbout(potE,tytul,ipdb)
1017 if (rstcount.eq.1000.or.itime.eq.n_timestep) then
1018 open(irest2,file=rest2name,status='unknown')
1019 write(irest2,*) totT,EK,potE,totE,t_bath
1021 ! AL 4/17/17: Now writing d_t(0,:) too
1023 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
1025 ! AL 4/17/17: Now writing d_c(0,:) too
1027 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
1035 t_MD=MPI_Wtime()-tt0
1039 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))') &
1041 'MD calculations setup:',t_MDsetup,&
1042 'Energy & gradient evaluation:',t_enegrad,&
1043 'Stochastic MD setup:',t_langsetup,&
1044 'Stochastic MD step setup:',t_sdsetup,&
1046 write (iout,'(/28(1h=),a25,27(1h=))') &
1047 ' End of MD calculation '
1049 write (iout,*) "time for etotal",t_etotal," elong",t_elong,&
1051 write (iout,*) "time_fric",time_fric," time_stoch",time_stoch,&
1052 " time_fricmatmult",time_fricmatmult," time_fsample ",&
1057 !-----------------------------------------------------------------------------
1058 subroutine velverlet_step(itime)
1059 !-------------------------------------------------------------------------------
1060 ! Perform a single velocity Verlet step; the time step can be rescaled if
1061 ! increments in accelerations exceed the threshold
1062 !-------------------------------------------------------------------------------
1063 ! implicit real*8 (a-h,o-z)
1064 ! include 'DIMENSIONS'
1066 use control, only:tcpu
1070 integer :: ierror,ierrcode
1071 real(kind=8) :: errcode
1073 ! include 'COMMON.SETUP'
1074 ! include 'COMMON.VAR'
1075 ! include 'COMMON.MD'
1077 ! include 'COMMON.LANGEVIN'
1079 ! include 'COMMON.LANGEVIN.lang0'
1081 ! include 'COMMON.CHAIN'
1082 ! include 'COMMON.DERIV'
1083 ! include 'COMMON.GEO'
1084 ! include 'COMMON.LOCAL'
1085 ! include 'COMMON.INTERACT'
1086 ! include 'COMMON.IOUNITS'
1087 ! include 'COMMON.NAMES'
1088 ! include 'COMMON.TIME1'
1089 ! include 'COMMON.MUCA'
1090 real(kind=8),dimension(3) :: vcm,incr
1091 real(kind=8),dimension(3) :: L
1092 integer :: count,rstcount !ilen,
1094 character(len=50) :: tytul
1095 integer :: maxcount_scale = 20
1096 !el common /gucio/ cm
1097 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1098 !el common /stochcalc/ stochforcvec
1099 integer :: itime,icount_scale,itime_scal,i,j,ifac_time
1101 real(kind=8) :: epdrift,tt0,fac_time
1103 if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres)) !(MAXRES6) maxres6=6*maxres
1109 else if (lang.eq.2 .or. lang.eq.3) then
1111 call stochastic_force(stochforcvec)
1114 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
1116 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1123 icount_scale=icount_scale+1
1124 if (icount_scale.gt.maxcount_scale) then
1126 "ERROR: too many attempts at scaling down the time step. ",&
1127 "amax=",amax,"epdrift=",epdrift,&
1128 "damax=",damax,"edriftmax=",edriftmax,&
1132 call MPI_Abort(MPI_COMM_WORLD,IERROR,IERRCODE)
1136 ! First step of the velocity Verlet algorithm
1141 else if (lang.eq.3) then
1143 call sd_verlet1_ciccotti
1145 else if (lang.eq.1) then
1150 ! Build the chain from the newly calculated coordinates
1151 call chainbuild_cart
1152 if (rattle) call rattle1
1154 if (large.and. mod(itime,ntwe).eq.0) then
1155 write (iout,*) "Cartesian and internal coordinates: step 1"
1160 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),&
1161 (dc(j,i+nres),j=1,3)
1163 write (iout,*) "Accelerations"
1165 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1166 (d_a(j,i+nres),j=1,3)
1168 write (iout,*) "Velocities, step 1"
1170 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1171 (d_t(j,i+nres),j=1,3)
1180 ! Calculate energy and forces
1182 call etotal(potEcomp)
1183 ! AL 4/17/17: Reduce the steps if NaNs occurred.
1184 if (potEcomp(0).gt.0.99e20 .or. isnan(potEcomp(0)).gt.0) then
1189 if (large.and. mod(itime,ntwe).eq.0) &
1190 call enerprint(potEcomp)
1193 t_etotal=t_etotal+MPI_Wtime()-tt0
1195 t_etotal=t_etotal+tcpu()-tt0
1198 potE=potEcomp(0)-potEcomp(20)
1200 ! Get the new accelerations
1203 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1205 t_enegrad=t_enegrad+tcpu()-tt0
1207 ! Determine maximum acceleration and scale down the timestep if needed
1209 amax=amax/(itime_scal+1)**2
1210 call predict_edrift(epdrift)
1211 if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then
1212 ! Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
1214 ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax)) &
1216 itime_scal=itime_scal+ifac_time
1217 ! fac_time=dmin1(damax/amax,0.5d0)
1218 fac_time=0.5d0**ifac_time
1219 d_time=d_time*fac_time
1220 if (lang.eq.2 .or. lang.eq.3) then
1222 ! write (iout,*) "Calling sd_verlet_setup: 1"
1223 ! Rescale the stochastic forces and recalculate or restore
1224 ! the matrices of tinker integrator
1225 if (itime_scal.gt.maxflag_stoch) then
1226 if (large) write (iout,'(a,i5,a)') &
1227 "Calculate matrices for stochastic step;",&
1228 " itime_scal ",itime_scal
1230 call sd_verlet_p_setup
1232 call sd_verlet_ciccotti_setup
1234 write (iout,'(2a,i3,a,i3,1h.)') &
1235 "Warning: cannot store matrices for stochastic",&
1236 " integration because the index",itime_scal,&
1237 " is greater than",maxflag_stoch
1238 write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",&
1239 " integration Langevin algorithm for better efficiency."
1240 else if (flag_stoch(itime_scal)) then
1241 if (large) write (iout,'(a,i5,a,l1)') &
1242 "Restore matrices for stochastic step; itime_scal ",&
1243 itime_scal," flag ",flag_stoch(itime_scal)
1246 pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
1247 afric_mat(i,j)=afric0_mat(i,j,itime_scal)
1248 vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
1249 prand_mat(i,j)=prand0_mat(i,j,itime_scal)
1250 vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
1251 vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
1255 if (large) write (iout,'(2a,i5,a,l1)') &
1256 "Calculate & store matrices for stochastic step;",&
1257 " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
1259 call sd_verlet_p_setup
1261 call sd_verlet_ciccotti_setup
1263 flag_stoch(ifac_time)=.true.
1266 pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
1267 afric0_mat(i,j,itime_scal)=afric_mat(i,j)
1268 vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
1269 prand0_mat(i,j,itime_scal)=prand_mat(i,j)
1270 vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
1271 vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
1275 fac_time=1.0d0/dsqrt(fac_time)
1277 stochforcvec(i)=fac_time*stochforcvec(i)
1280 else if (lang.eq.1) then
1281 ! Rescale the accelerations due to stochastic forces
1282 fac_time=1.0d0/dsqrt(fac_time)
1284 d_as_work(i)=d_as_work(i)*fac_time
1287 if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)') &
1288 "itime",itime," Timestep scaled down to ",&
1289 d_time," ifac_time",ifac_time," itime_scal",itime_scal
1291 ! Second step of the velocity Verlet algorithm
1296 else if (lang.eq.3) then
1298 call sd_verlet2_ciccotti
1300 else if (lang.eq.1) then
1305 if (rattle) call rattle2
1308 if (d_time.ne.d_time0) then
1311 if (lang.eq.2 .or. lang.eq.3) then
1312 if (large) write (iout,'(a)') &
1313 "Restore original matrices for stochastic step"
1314 ! write (iout,*) "Calling sd_verlet_setup: 2"
1315 ! Restore the matrices of tinker integrator if the time step has been restored
1318 pfric_mat(i,j)=pfric0_mat(i,j,0)
1319 afric_mat(i,j)=afric0_mat(i,j,0)
1320 vfric_mat(i,j)=vfric0_mat(i,j,0)
1321 prand_mat(i,j)=prand0_mat(i,j,0)
1322 vrand_mat1(i,j)=vrand0_mat1(i,j,0)
1323 vrand_mat2(i,j)=vrand0_mat2(i,j,0)
1332 ! Calculate the kinetic and the total energy and the kinetic temperature
1336 ! call kinetic1(EK1)
1337 ! write (iout,*) "step",itime," EK",EK," EK1",EK1
1339 ! Couple the system to Berendsen bath if needed
1340 if (tbf .and. lang.eq.0) then
1343 kinetic_T=2.0d0/(dimen3*Rb)*EK
1344 ! Backup the coordinates, velocities, and accelerations
1348 d_t_old(j,i)=d_t(j,i)
1349 d_a_old(j,i)=d_a(j,i)
1353 if (mod(itime,ntwe).eq.0 .and. large) then
1354 write (iout,*) "Velocities, step 2"
1356 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1357 (d_t(j,i+nres),j=1,3)
1362 end subroutine velverlet_step
1363 !-----------------------------------------------------------------------------
1364 subroutine RESPA_step(itime)
1365 !-------------------------------------------------------------------------------
1366 ! Perform a single RESPA step.
1367 !-------------------------------------------------------------------------------
1368 ! implicit real*8 (a-h,o-z)
1369 ! include 'DIMENSIONS'
1373 use control, only:tcpu
1375 ! use io_conf, only:cartprint
1378 integer :: IERROR,ERRCODE
1380 ! include 'COMMON.SETUP'
1381 ! include 'COMMON.CONTROL'
1382 ! include 'COMMON.VAR'
1383 ! include 'COMMON.MD'
1385 ! include 'COMMON.LANGEVIN'
1387 ! include 'COMMON.LANGEVIN.lang0'
1389 ! include 'COMMON.CHAIN'
1390 ! include 'COMMON.DERIV'
1391 ! include 'COMMON.GEO'
1392 ! include 'COMMON.LOCAL'
1393 ! include 'COMMON.INTERACT'
1394 ! include 'COMMON.IOUNITS'
1395 ! include 'COMMON.NAMES'
1396 ! include 'COMMON.TIME1'
1397 real(kind=8),dimension(0:n_ene) :: energia_short,energia_long
1398 real(kind=8),dimension(3) :: L,vcm,incr
1399 real(kind=8),dimension(3,0:2*nres) :: dc_old0,d_t_old0,d_a_old0 !(3,0:maxres2) maxres2=2*maxres
1400 logical :: PRINT_AMTS_MSG = .false.
1401 integer :: count,rstcount !ilen,
1403 character(len=50) :: tytul
1404 integer :: maxcount_scale = 10
1405 !el common /gucio/ cm
1406 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1407 !el common /stochcalc/ stochforcvec
1408 integer :: itime,itt,i,j,itsplit
1410 !el common /cipiszcze/ itt
1412 real(kind=8) :: epdrift,tt0,epdriftmax
1415 if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres)) !(MAXRES6) maxres6=6*maxres
1419 if (large.and. mod(itime,ntwe).eq.0) then
1420 write (iout,*) "***************** RESPA itime",itime
1421 write (iout,*) "Cartesian and internal coordinates: step 0"
1423 call pdbout(0.0d0,"cipiszcze",iout)
1425 write (iout,*) "Accelerations from long-range forces"
1427 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1428 (d_a(j,i+nres),j=1,3)
1430 write (iout,*) "Velocities, step 0"
1432 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1433 (d_t(j,i+nres),j=1,3)
1438 ! Perform the initial RESPA step (increment velocities)
1439 ! write (iout,*) "*********************** RESPA ini"
1442 if (mod(itime,ntwe).eq.0 .and. large) then
1443 write (iout,*) "Velocities, end"
1445 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1446 (d_t(j,i+nres),j=1,3)
1450 ! Compute the short-range forces
1456 ! 7/2/2009 commented out
1458 ! call etotal_short(energia_short)
1461 ! 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
1464 d_a(j,i)=d_a_short(j,i)
1468 if (large.and. mod(itime,ntwe).eq.0) then
1469 write (iout,*) "energia_short",energia_short(0)
1470 write (iout,*) "Accelerations from short-range forces"
1472 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1473 (d_a(j,i+nres),j=1,3)
1478 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1480 t_enegrad=t_enegrad+tcpu()-tt0
1485 d_t_old(j,i)=d_t(j,i)
1486 d_a_old(j,i)=d_a(j,i)
1489 ! 6/30/08 A-MTS: attempt at increasing the split number
1492 dc_old0(j,i)=dc_old(j,i)
1493 d_t_old0(j,i)=d_t_old(j,i)
1494 d_a_old0(j,i)=d_a_old(j,i)
1497 if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
1498 if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
1505 ! write (iout,*) "itime",itime," ntime_split",ntime_split
1506 ! Split the time step
1507 d_time=d_time0/ntime_split
1508 ! Perform the short-range RESPA steps (velocity Verlet increments of
1509 ! positions and velocities using short-range forces)
1510 ! write (iout,*) "*********************** RESPA split"
1511 do itsplit=1,ntime_split
1514 else if (lang.eq.2 .or. lang.eq.3) then
1516 call stochastic_force(stochforcvec)
1519 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
1521 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1526 ! First step of the velocity Verlet algorithm
1531 else if (lang.eq.3) then
1533 call sd_verlet1_ciccotti
1535 else if (lang.eq.1) then
1540 ! Build the chain from the newly calculated coordinates
1541 call chainbuild_cart
1542 if (rattle) call rattle1
1544 if (large.and. mod(itime,ntwe).eq.0) then
1545 write (iout,*) "***** ITSPLIT",itsplit
1546 write (iout,*) "Cartesian and internal coordinates: step 1"
1547 call pdbout(0.0d0,"cipiszcze",iout)
1550 write (iout,*) "Velocities, step 1"
1552 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1553 (d_t(j,i+nres),j=1,3)
1562 ! Calculate energy and forces
1564 call etotal_short(energia_short)
1565 ! AL 4/17/17: Exit itime_split loop when energy goes infinite
1566 if (energia_short(0).gt.0.99e20 .or. isnan(energia_short(0)) ) then
1567 if (PRINT_AMTS_MSG) &
1568 write (iout,*) "Infinities/NaNs in energia_short",energia_short(0),"; increasing ntime_split to",ntime_split
1569 ntime_split=ntime_split*2
1570 if (ntime_split.gt.maxtime_split) then
1573 "Cannot rescue the run; aborting job. Retry with a smaller time step"
1575 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1578 "Cannot rescue the run; terminating. Retry with a smaller time step"
1584 if (large.and. mod(itime,ntwe).eq.0) &
1585 call enerprint(energia_short)
1588 t_eshort=t_eshort+MPI_Wtime()-tt0
1590 t_eshort=t_eshort+tcpu()-tt0
1594 ! Get the new accelerations
1596 ! 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
1599 d_a_short(j,i)=d_a(j,i)
1603 if (large.and. mod(itime,ntwe).eq.0) then
1604 write (iout,*)"energia_short",energia_short(0)
1605 write (iout,*) "Accelerations from short-range forces"
1607 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1608 (d_a(j,i+nres),j=1,3)
1613 ! Determine maximum acceleration and scale down the timestep if needed
1615 amax=amax/ntime_split**2
1616 call predict_edrift(epdrift)
1617 if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0) &
1618 write (iout,*) "amax",amax," damax",damax,&
1619 " epdrift",epdrift," epdriftmax",epdriftmax
1620 ! Exit loop and try with increased split number if the change of
1621 ! acceleration is too big
1622 if (amax.gt.damax .or. epdrift.gt.edriftmax) then
1623 if (ntime_split.lt.maxtime_split) then
1625 ntime_split=ntime_split*2
1626 ! AL 4/17/17: We should exit the itime_split loop when acceleration change is too big
1630 dc_old(j,i)=dc_old0(j,i)
1631 d_t_old(j,i)=d_t_old0(j,i)
1632 d_a_old(j,i)=d_a_old0(j,i)
1635 if (PRINT_AMTS_MSG) then
1636 write (iout,*) "acceleration/energy drift too large",amax,&
1637 epdrift," split increased to ",ntime_split," itime",itime,&
1643 "Uh-hu. Bumpy landscape. Maximum splitting number",&
1645 " already reached!!! Trying to carry on!"
1649 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1651 t_enegrad=t_enegrad+tcpu()-tt0
1653 ! Second step of the velocity Verlet algorithm
1658 else if (lang.eq.3) then
1660 call sd_verlet2_ciccotti
1662 else if (lang.eq.1) then
1667 if (rattle) call rattle2
1668 ! Backup the coordinates, velocities, and accelerations
1672 d_t_old(j,i)=d_t(j,i)
1673 d_a_old(j,i)=d_a(j,i)
1680 ! Restore the time step
1682 ! Compute long-range forces
1689 call etotal_long(energia_long)
1690 if (energia_long(0).gt.0.99e20 .or. isnan(energia_long(0))) then
1693 "Infinitied/NaNs in energia_long, Aborting MPI job."
1695 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1697 write (iout,*) "Infinitied/NaNs in energia_long, terminating."
1701 if (large.and. mod(itime,ntwe).eq.0) &
1702 call enerprint(energia_long)
1705 t_elong=t_elong+MPI_Wtime()-tt0
1707 t_elong=t_elong+tcpu()-tt0
1713 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1715 t_enegrad=t_enegrad+tcpu()-tt0
1717 ! Compute accelerations from long-range forces
1719 if (large.and. mod(itime,ntwe).eq.0) then
1720 write (iout,*) "energia_long",energia_long(0)
1721 write (iout,*) "Cartesian and internal coordinates: step 2"
1723 call pdbout(0.0d0,"cipiszcze",iout)
1725 write (iout,*) "Accelerations from long-range forces"
1727 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1728 (d_a(j,i+nres),j=1,3)
1730 write (iout,*) "Velocities, step 2"
1732 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1733 (d_t(j,i+nres),j=1,3)
1737 ! Compute the final RESPA step (increment velocities)
1738 ! write (iout,*) "*********************** RESPA fin"
1740 ! Compute the complete potential energy
1742 potEcomp(i)=energia_short(i)+energia_long(i)
1744 potE=potEcomp(0)-potEcomp(20)
1745 ! potE=energia_short(0)+energia_long(0)
1748 ! Calculate the kinetic and the total energy and the kinetic temperature
1751 ! Couple the system to Berendsen bath if needed
1752 if (tbf .and. lang.eq.0) then
1755 kinetic_T=2.0d0/(dimen3*Rb)*EK
1756 ! Backup the coordinates, velocities, and accelerations
1758 if (mod(itime,ntwe).eq.0 .and. large) then
1759 write (iout,*) "Velocities, end"
1761 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1762 (d_t(j,i+nres),j=1,3)
1767 end subroutine RESPA_step
1768 !-----------------------------------------------------------------------------
1769 subroutine RESPA_vel
1770 ! First and last RESPA step (incrementing velocities using long-range
1773 ! implicit real*8 (a-h,o-z)
1774 ! include 'DIMENSIONS'
1775 ! include 'COMMON.CONTROL'
1776 ! include 'COMMON.VAR'
1777 ! include 'COMMON.MD'
1778 ! include 'COMMON.CHAIN'
1779 ! include 'COMMON.DERIV'
1780 ! include 'COMMON.GEO'
1781 ! include 'COMMON.LOCAL'
1782 ! include 'COMMON.INTERACT'
1783 ! include 'COMMON.IOUNITS'
1784 ! include 'COMMON.NAMES'
1785 integer :: i,j,inres,mnum
1788 d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
1792 d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
1797 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
1798 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
1799 .and.(mnum.ne.5)) then
1802 d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
1807 end subroutine RESPA_vel
1808 !-----------------------------------------------------------------------------
1810 ! Applying velocity Verlet algorithm - step 1 to coordinates
1812 ! implicit real*8 (a-h,o-z)
1813 ! include 'DIMENSIONS'
1814 ! include 'COMMON.CONTROL'
1815 ! include 'COMMON.VAR'
1816 ! include 'COMMON.MD'
1817 ! include 'COMMON.CHAIN'
1818 ! include 'COMMON.DERIV'
1819 ! include 'COMMON.GEO'
1820 ! include 'COMMON.LOCAL'
1821 ! include 'COMMON.INTERACT'
1822 ! include 'COMMON.IOUNITS'
1823 ! include 'COMMON.NAMES'
1824 real(kind=8) :: adt,adt2
1825 integer :: i,j,inres,mnum
1828 write (iout,*) "VELVERLET1 START: DC"
1830 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),&
1831 (dc(j,i+nres),j=1,3)
1835 adt=d_a_old(j,0)*d_time
1837 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1838 d_t_new(j,0)=d_t_old(j,0)+adt2
1839 d_t(j,0)=d_t_old(j,0)+adt
1843 adt=d_a_old(j,i)*d_time
1845 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1846 d_t_new(j,i)=d_t_old(j,i)+adt2
1847 d_t(j,i)=d_t_old(j,i)+adt
1852 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
1853 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
1854 .and.(mnum.ne.5)) then
1857 adt=d_a_old(j,inres)*d_time
1859 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1860 d_t_new(j,inres)=d_t_old(j,inres)+adt2
1861 d_t(j,inres)=d_t_old(j,inres)+adt
1866 write (iout,*) "VELVERLET1 END: DC"
1868 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),&
1869 (dc(j,i+nres),j=1,3)
1873 end subroutine verlet1
1874 !-----------------------------------------------------------------------------
1876 ! Step 2 of the velocity Verlet algorithm: update velocities
1878 ! implicit real*8 (a-h,o-z)
1879 ! include 'DIMENSIONS'
1880 ! include 'COMMON.CONTROL'
1881 ! include 'COMMON.VAR'
1882 ! include 'COMMON.MD'
1883 ! include 'COMMON.CHAIN'
1884 ! include 'COMMON.DERIV'
1885 ! include 'COMMON.GEO'
1886 ! include 'COMMON.LOCAL'
1887 ! include 'COMMON.INTERACT'
1888 ! include 'COMMON.IOUNITS'
1889 ! include 'COMMON.NAMES'
1890 integer :: i,j,inres,mnum
1893 d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
1897 d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
1902 ! iti=iabs(itype(i,mnum))
1903 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
1904 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
1905 .and.(mnum.ne.5)) then
1908 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
1913 end subroutine verlet2
1914 !-----------------------------------------------------------------------------
1915 subroutine sddir_precalc
1916 ! Applying velocity Verlet algorithm - step 1 to coordinates
1917 ! implicit real*8 (a-h,o-z)
1918 ! include 'DIMENSIONS'
1924 ! include 'COMMON.CONTROL'
1925 ! include 'COMMON.VAR'
1926 ! include 'COMMON.MD'
1928 ! include 'COMMON.LANGEVIN'
1930 ! include 'COMMON.LANGEVIN.lang0'
1932 ! include 'COMMON.CHAIN'
1933 ! include 'COMMON.DERIV'
1934 ! include 'COMMON.GEO'
1935 ! include 'COMMON.LOCAL'
1936 ! include 'COMMON.INTERACT'
1937 ! include 'COMMON.IOUNITS'
1938 ! include 'COMMON.NAMES'
1939 ! include 'COMMON.TIME1'
1940 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1941 !el common /stochcalc/ stochforcvec
1942 real(kind=8) :: time00
1944 ! Compute friction and stochastic forces
1949 time_fric=time_fric+MPI_Wtime()-time00
1951 call stochastic_force(stochforcvec)
1952 time_stoch=time_stoch+MPI_Wtime()-time00
1955 ! Compute the acceleration due to friction forces (d_af_work) and stochastic
1956 ! forces (d_as_work)
1958 call ginv_mult(fric_work, d_af_work)
1959 call ginv_mult(stochforcvec, d_as_work)
1961 end subroutine sddir_precalc
1962 !-----------------------------------------------------------------------------
1963 subroutine sddir_verlet1
1964 ! Applying velocity Verlet algorithm - step 1 to velocities
1967 ! implicit real*8 (a-h,o-z)
1968 ! include 'DIMENSIONS'
1969 ! include 'COMMON.CONTROL'
1970 ! include 'COMMON.VAR'
1971 ! include 'COMMON.MD'
1973 ! include 'COMMON.LANGEVIN'
1975 ! include 'COMMON.LANGEVIN.lang0'
1977 ! include 'COMMON.CHAIN'
1978 ! include 'COMMON.DERIV'
1979 ! include 'COMMON.GEO'
1980 ! include 'COMMON.LOCAL'
1981 ! include 'COMMON.INTERACT'
1982 ! include 'COMMON.IOUNITS'
1983 ! include 'COMMON.NAMES'
1984 ! Revised 3/31/05 AL: correlation between random contributions to
1985 ! position and velocity increments included.
1986 real(kind=8) :: sqrt13 = 0.57735026918962576451d0 ! 1/sqrt(3)
1987 real(kind=8) :: adt,adt2
1988 integer :: i,j,ind,inres,mnum
1990 ! Add the contribution from BOTH friction and stochastic force to the
1991 ! coordinates, but ONLY the contribution from the friction forces to velocities
1994 adt=(d_a_old(j,0)+d_af_work(j))*d_time
1995 adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
1996 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1997 d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
1998 d_t(j,0)=d_t_old(j,0)+adt
2003 adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
2004 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
2005 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
2006 d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
2007 d_t(j,i)=d_t_old(j,i)+adt
2013 ! iti=iabs(itype(i,mnum))
2014 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2015 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2016 .and.(mnum.ne.5)) then
2019 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
2020 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
2021 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
2022 d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
2023 d_t(j,inres)=d_t_old(j,inres)+adt
2029 end subroutine sddir_verlet1
2030 !-----------------------------------------------------------------------------
2031 subroutine sddir_verlet2
2032 ! Calculating the adjusted velocities for accelerations
2035 ! implicit real*8 (a-h,o-z)
2036 ! include 'DIMENSIONS'
2037 ! include 'COMMON.CONTROL'
2038 ! include 'COMMON.VAR'
2039 ! include 'COMMON.MD'
2041 ! include 'COMMON.LANGEVIN'
2043 ! include 'COMMON.LANGEVIN.lang0'
2045 ! include 'COMMON.CHAIN'
2046 ! include 'COMMON.DERIV'
2047 ! include 'COMMON.GEO'
2048 ! include 'COMMON.LOCAL'
2049 ! include 'COMMON.INTERACT'
2050 ! include 'COMMON.IOUNITS'
2051 ! include 'COMMON.NAMES'
2052 real(kind=8),dimension(6*nres) :: stochforcvec,d_as_work1 !(MAXRES6) maxres6=6*maxres
2053 real(kind=8) :: cos60 = 0.5d0, sin60 = 0.86602540378443864676d0
2054 integer :: i,j,ind,inres,mnum
2055 ! Revised 3/31/05 AL: correlation between random contributions to
2056 ! position and velocity increments included.
2057 ! The correlation coefficients are calculated at low-friction limit.
2058 ! Also, friction forces are now not calculated with new velocities.
2060 ! call friction_force
2061 call stochastic_force(stochforcvec)
2063 ! Compute the acceleration due to friction forces (d_af_work) and stochastic
2064 ! forces (d_as_work)
2066 call ginv_mult(stochforcvec, d_as_work1)
2072 d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j)) &
2073 +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
2078 d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j)) &
2079 +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
2085 ! iti=iabs(itype(i,mnum))
2086 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2087 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2088 .and.(mnum.ne.5)) then
2091 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres) &
2092 +d_af_work(ind+j))+sin60*d_as_work(ind+j) &
2093 +cos60*d_as_work1(ind+j))*d_time
2099 end subroutine sddir_verlet2
2100 !-----------------------------------------------------------------------------
2101 subroutine max_accel
2103 ! Find the maximum difference in the accelerations of the the sites
2104 ! at the beginning and the end of the time step.
2107 ! implicit real*8 (a-h,o-z)
2108 ! include 'DIMENSIONS'
2109 ! include 'COMMON.CONTROL'
2110 ! include 'COMMON.VAR'
2111 ! include 'COMMON.MD'
2112 ! include 'COMMON.CHAIN'
2113 ! include 'COMMON.DERIV'
2114 ! include 'COMMON.GEO'
2115 ! include 'COMMON.LOCAL'
2116 ! include 'COMMON.INTERACT'
2117 ! include 'COMMON.IOUNITS'
2118 real(kind=8),dimension(3) :: aux,accel,accel_old
2119 real(kind=8) :: dacc
2123 ! aux(j)=d_a(j,0)-d_a_old(j,0)
2124 accel_old(j)=d_a_old(j,0)
2131 ! 7/3/08 changed to asymmetric difference
2133 ! accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
2134 accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
2135 accel(j)=accel(j)+0.5d0*d_a(j,i)
2136 ! if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
2137 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
2138 dacc=dabs(accel(j)-accel_old(j))
2139 ! write (iout,*) i,dacc
2140 if (dacc.gt.amax) amax=dacc
2148 accel_old(j)=d_a_old(j,0)
2153 accel_old(j)=accel_old(j)+d_a_old(j,1)
2154 accel(j)=accel(j)+d_a(j,1)
2159 ! iti=iabs(itype(i,mnum))
2160 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2161 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2162 .and.(mnum.ne.5)) then
2164 ! accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
2165 accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
2166 accel(j)=accel(j)+d_a(j,i+nres)
2170 ! if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
2171 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
2172 dacc=dabs(accel(j)-accel_old(j))
2173 ! write (iout,*) "side-chain",i,dacc
2174 if (dacc.gt.amax) amax=dacc
2178 accel_old(j)=accel_old(j)+d_a_old(j,i)
2179 accel(j)=accel(j)+d_a(j,i)
2180 ! aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
2184 end subroutine max_accel
2185 !-----------------------------------------------------------------------------
2186 subroutine predict_edrift(epdrift)
2188 ! Predict the drift of the potential energy
2191 use control_data, only: lmuca
2192 ! implicit real*8 (a-h,o-z)
2193 ! include 'DIMENSIONS'
2194 ! include 'COMMON.CONTROL'
2195 ! include 'COMMON.VAR'
2196 ! include 'COMMON.MD'
2197 ! include 'COMMON.CHAIN'
2198 ! include 'COMMON.DERIV'
2199 ! include 'COMMON.GEO'
2200 ! include 'COMMON.LOCAL'
2201 ! include 'COMMON.INTERACT'
2202 ! include 'COMMON.IOUNITS'
2203 ! include 'COMMON.MUCA'
2204 real(kind=8) :: epdrift,epdriftij
2206 ! Drift of the potential energy
2212 epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
2213 if (lmuca) epdriftij=epdriftij*factor
2214 ! write (iout,*) "back",i,j,epdriftij
2215 if (epdriftij.gt.epdrift) epdrift=epdriftij
2219 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1.and.&
2220 molnum(i).ne.5) then
2223 dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
2224 if (lmuca) epdriftij=epdriftij*factor
2225 ! write (iout,*) "side",i,j,epdriftij
2226 if (epdriftij.gt.epdrift) epdrift=epdriftij
2230 epdrift=0.5d0*epdrift*d_time*d_time
2231 ! write (iout,*) "epdrift",epdrift
2233 end subroutine predict_edrift
2234 !-----------------------------------------------------------------------------
2235 subroutine verlet_bath
2237 ! Coupling to the thermostat by using the Berendsen algorithm
2240 ! implicit real*8 (a-h,o-z)
2241 ! include 'DIMENSIONS'
2242 ! include 'COMMON.CONTROL'
2243 ! include 'COMMON.VAR'
2244 ! include 'COMMON.MD'
2245 ! include 'COMMON.CHAIN'
2246 ! include 'COMMON.DERIV'
2247 ! include 'COMMON.GEO'
2248 ! include 'COMMON.LOCAL'
2249 ! include 'COMMON.INTERACT'
2250 ! include 'COMMON.IOUNITS'
2251 ! include 'COMMON.NAMES'
2252 real(kind=8) :: T_half,fact
2253 integer :: i,j,inres,mnum
2255 T_half=2.0d0/(dimen3*Rb)*EK
2256 fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
2257 ! write(iout,*) "T_half", T_half
2258 ! write(iout,*) "EK", EK
2259 ! write(iout,*) "fact", fact
2261 d_t(j,0)=fact*d_t(j,0)
2265 d_t(j,i)=fact*d_t(j,i)
2270 ! iti=iabs(itype(i,mnum))
2271 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2272 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2273 .and.(mnum.ne.5)) then
2276 d_t(j,inres)=fact*d_t(j,inres)
2281 end subroutine verlet_bath
2282 !-----------------------------------------------------------------------------
2284 ! Set up the initial conditions of a MD simulation
2287 use control, only:tcpu
2288 !el use io_basic, only:ilen
2291 use minimm, only:minim_dc,minimize,sc_move
2292 use io_config, only:readrst
2293 use io, only:statout
2294 ! implicit real*8 (a-h,o-z)
2295 ! include 'DIMENSIONS'
2298 character(len=16) :: form
2299 integer :: IERROR,ERRCODE
2301 ! include 'COMMON.SETUP'
2302 ! include 'COMMON.CONTROL'
2303 ! include 'COMMON.VAR'
2304 ! include 'COMMON.MD'
2306 ! include 'COMMON.LANGEVIN'
2308 ! include 'COMMON.LANGEVIN.lang0'
2310 ! include 'COMMON.CHAIN'
2311 ! include 'COMMON.DERIV'
2312 ! include 'COMMON.GEO'
2313 ! include 'COMMON.LOCAL'
2314 ! include 'COMMON.INTERACT'
2315 ! include 'COMMON.IOUNITS'
2316 ! include 'COMMON.NAMES'
2317 ! include 'COMMON.REMD'
2318 real(kind=8),dimension(0:n_ene) :: energia_long,energia_short
2319 real(kind=8),dimension(3) :: vcm,incr,L
2320 real(kind=8) :: xv,sigv,lowb,highb
2321 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
2322 character(len=256) :: qstr
2325 character(len=50) :: tytul
2326 logical :: file_exist
2327 !el common /gucio/ cm
2328 integer :: i,j,ipos,iq,iw,nft_sc,iretcode,nfun,itime,ierr,mnum
2329 real(kind=8) :: etot,tt0
2333 ! write(iout,*) "d_time", d_time
2334 ! Compute the standard deviations of stochastic forces for Langevin dynamics
2335 ! if the friction coefficients do not depend on surface area
2336 if (lang.gt.0 .and. .not.surfarea) then
2339 stdforcp(i)=stdfp(mnum)*dsqrt(gamp(mnum))
2343 stdforcsc(i)=stdfsc(iabs(itype(i,mnum)),mnum) &
2344 *dsqrt(gamsc(iabs(itype(i,mnum)),mnum))
2347 ! Open the pdb file for snapshotshots
2350 if (ilen(tmpdir).gt.0) &
2351 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2352 liczba(:ilen(liczba))//".pdb")
2354 file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2358 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) &
2359 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2360 liczba(:ilen(liczba))//".x")
2361 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2364 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) &
2365 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2366 liczba(:ilen(liczba))//".cx")
2367 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2373 if (ilen(tmpdir).gt.0) &
2374 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
2375 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
2377 if (ilen(tmpdir).gt.0) &
2378 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
2379 cartname=prefix(:ilen(prefix))//"_MD.cx"
2383 write (qstr,'(256(1h ))')
2386 iq = qinfrag(i,iset)*10
2387 iw = wfrag(i,iset)/100
2389 if(me.eq.king.or..not.out1file) &
2390 write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
2391 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
2396 iq = qinpair(i,iset)*10
2397 iw = wpair(i,iset)/100
2399 if(me.eq.king.or..not.out1file) &
2400 write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
2401 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
2405 ! pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
2407 ! cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
2409 ! cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
2411 ! statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
2415 if (restart1file) then
2417 inquire(file=mremd_rst_name,exist=file_exist)
2418 write (*,*) me," Before broadcast: file_exist",file_exist
2420 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,&
2423 write (*,*) me," After broadcast: file_exist",file_exist
2424 ! inquire(file=mremd_rst_name,exist=file_exist)
2425 if(me.eq.king.or..not.out1file) &
2426 write(iout,*) "Initial state read by master and distributed"
2428 if (ilen(tmpdir).gt.0) &
2429 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_' &
2430 //liczba(:ilen(liczba))//'.rst')
2431 inquire(file=rest2name,exist=file_exist)
2434 if(.not.restart1file) then
2435 if(me.eq.king.or..not.out1file) &
2436 write(iout,*) "Initial state will be read from file ",&
2437 rest2name(:ilen(rest2name))
2440 call rescale_weights(t_bath)
2442 if(me.eq.king.or..not.out1file)then
2443 if (restart1file) then
2444 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),&
2447 write(iout,*) "File ",rest2name(:ilen(rest2name)),&
2450 write(iout,*) "Initial velocities randomly generated"
2457 ! Generate initial velocities
2458 if(me.eq.king.or..not.out1file) &
2459 write(iout,*) "Initial velocities randomly generated"
2464 ! rest2name = prefix(:ilen(prefix))//'.rst'
2465 if(me.eq.king.or..not.out1file)then
2466 write (iout,*) "Initial velocities"
2468 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2469 (d_t(j,i+nres),j=1,3)
2471 ! Zeroing the total angular momentum of the system
2472 write(iout,*) "Calling the zero-angular momentum subroutine"
2475 ! Getting the potential energy and forces and velocities and accelerations
2477 ! write (iout,*) "velocity of the center of the mass:"
2478 ! write (iout,*) (vcm(j),j=1,3)
2480 d_t(j,0)=d_t(j,0)-vcm(j)
2482 ! Removing the velocity of the center of mass
2484 if(me.eq.king.or..not.out1file)then
2485 write (iout,*) "vcm right after adjustment:"
2486 write (iout,*) (vcm(j),j=1,3)
2488 if ((.not.rest).and.(indpdb.eq.0)) then
2490 if(iranconf.ne.0) then
2492 print *, 'Calling OVERLAP_SC'
2493 call overlap_sc(fail)
2496 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
2497 print *,'SC_move',nft_sc,etot
2498 if(me.eq.king.or..not.out1file) &
2499 write(iout,*) 'SC_move',nft_sc,etot
2503 print *, 'Calling MINIM_DC'
2504 call minim_dc(etot,iretcode,nfun)
2506 call geom_to_var(nvar,varia)
2507 print *,'Calling MINIMIZE.'
2508 call minimize(etot,varia,iretcode,nfun)
2509 call var_to_geom(nvar,varia)
2511 if(me.eq.king.or..not.out1file) &
2512 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2515 call chainbuild_cart
2520 kinetic_T=2.0d0/(dimen3*Rb)*EK
2521 if(me.eq.king.or..not.out1file)then
2531 call etotal(potEcomp)
2532 if (large) call enerprint(potEcomp)
2535 t_etotal=t_etotal+MPI_Wtime()-tt0
2537 t_etotal=t_etotal+tcpu()-tt0
2544 if (amax*d_time .gt. dvmax) then
2545 d_time=d_time*dvmax/amax
2546 if(me.eq.king.or..not.out1file) write (iout,*) &
2547 "Time step reduced to",d_time,&
2548 " because of too large initial acceleration."
2550 if(me.eq.king.or..not.out1file)then
2551 write(iout,*) "Potential energy and its components"
2552 call enerprint(potEcomp)
2553 ! write(iout,*) (potEcomp(i),i=0,n_ene)
2555 potE=potEcomp(0)-potEcomp(20)
2558 if (ntwe.ne.0) call statout(itime)
2559 if(me.eq.king.or..not.out1file) &
2560 write (iout,'(/a/3(a25,1pe14.5/))') "Initial:", &
2561 " Kinetic energy",EK," Potential energy",potE, &
2562 " Total energy",totE," Maximum acceleration ", &
2565 write (iout,*) "Initial coordinates"
2567 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),&
2570 write (iout,*) "Initial dC"
2572 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),&
2573 (dc(j,i+nres),j=1,3)
2575 write (iout,*) "Initial velocities"
2576 write (iout,"(13x,' backbone ',23x,' side chain')")
2578 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2579 (d_t(j,i+nres),j=1,3)
2581 write (iout,*) "Initial accelerations"
2583 ! write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2584 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),&
2585 (d_a(j,i+nres),j=1,3)
2591 d_t_old(j,i)=d_t(j,i)
2592 d_a_old(j,i)=d_a(j,i)
2594 ! write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
2603 call etotal_short(energia_short)
2604 if (large) call enerprint(potEcomp)
2607 t_eshort=t_eshort+MPI_Wtime()-tt0
2609 t_eshort=t_eshort+tcpu()-tt0
2614 if(.not.out1file .and. large) then
2615 write (iout,*) "energia_long",energia_long(0),&
2616 " energia_short",energia_short(0),&
2617 " total",energia_long(0)+energia_short(0)
2618 write (iout,*) "Initial fast-force accelerations"
2620 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
2621 (d_a(j,i+nres),j=1,3)
2624 ! 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
2627 d_a_short(j,i)=d_a(j,i)
2636 call etotal_long(energia_long)
2637 if (large) call enerprint(potEcomp)
2640 t_elong=t_elong+MPI_Wtime()-tt0
2642 t_elong=t_elong+tcpu()-tt0
2647 if(.not.out1file .and. large) then
2648 write (iout,*) "energia_long",energia_long(0)
2649 write (iout,*) "Initial slow-force accelerations"
2651 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
2652 (d_a(j,i+nres),j=1,3)
2656 t_enegrad=t_enegrad+MPI_Wtime()-tt0
2658 t_enegrad=t_enegrad+tcpu()-tt0
2662 end subroutine init_MD
2663 !-----------------------------------------------------------------------------
2664 subroutine random_vel
2666 ! implicit real*8 (a-h,o-z)
2668 use random, only:anorm_distr
2670 ! include 'DIMENSIONS'
2671 ! include 'COMMON.CONTROL'
2672 ! include 'COMMON.VAR'
2673 ! include 'COMMON.MD'
2675 ! include 'COMMON.LANGEVIN'
2677 ! include 'COMMON.LANGEVIN.lang0'
2679 ! include 'COMMON.CHAIN'
2680 ! include 'COMMON.DERIV'
2681 ! include 'COMMON.GEO'
2682 ! include 'COMMON.LOCAL'
2683 ! include 'COMMON.INTERACT'
2684 ! include 'COMMON.IOUNITS'
2685 ! include 'COMMON.NAMES'
2686 ! include 'COMMON.TIME1'
2687 real(kind=8) :: xv,sigv,lowb,highb ,Ek1
2689 real(kind=8) ,allocatable, dimension(:) :: DDU1,DDU2,DL2,DL1,xsolv,DML,rs
2690 real(kind=8) :: sumx
2692 real(kind=8) ,allocatable, dimension(:) :: rsold
2693 real (kind=8),allocatable,dimension(:,:) :: matold
2696 integer :: i,j,ii,k,ind,mark,imark,mnum
2697 ! Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
2698 ! First generate velocities in the eigenspace of the G matrix
2699 ! write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
2706 sigv=dsqrt((Rb*t_bath)/geigen(i))
2709 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
2711 write (iout,*) "i",i," ii",ii," geigen",geigen(i),&
2712 " d_t_work_new",d_t_work_new(ii)
2723 Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
2726 write (iout,*) "Ek from eigenvectors",Ek1
2727 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
2731 ! Transform velocities to UNRES coordinate space
2732 allocate (DL1(2*nres))
2733 allocate (DDU1(2*nres))
2734 allocate (DL2(2*nres))
2735 allocate (DDU2(2*nres))
2736 allocate (xsolv(2*nres))
2737 allocate (DML(2*nres))
2738 allocate (rs(2*nres))
2740 allocate (rsold(2*nres))
2741 allocate (matold(2*nres,2*nres))
2743 matold(1,1)=DMorig(1)
2744 matold(1,2)=DU1orig(1)
2745 matold(1,3)=DU2orig(1)
2746 write (*,*) DMorig(1),DU1orig(1),DU2orig(1)
2751 matold(i,j)=DMorig(i)
2752 matold(i,j-1)=DU1orig(i-1)
2754 matold(i,j-2)=DU2orig(i-2)
2762 matold(i,j+1)=DU1orig(i)
2768 matold(i,j+2)=DU2orig(i)
2772 matold(dimen,dimen)=DMorig(dimen)
2773 matold(dimen,dimen-1)=DU1orig(dimen-1)
2774 matold(dimen,dimen-2)=DU2orig(dimen-2)
2775 write(iout,*) "old gmatrix"
2776 call matout(dimen,dimen,2*nres,2*nres,matold)
2780 ! Find the ith eigenvector of the pentadiagonal inertiq matrix
2784 DML(j)=DMorig(j)-geigen(i)
2787 DML(j-1)=DMorig(j)-geigen(i)
2792 DDU1(imark-1)=DU2orig(imark-1)
2793 do j=imark+1,dimen-1
2794 DDU1(j-1)=DU1orig(j)
2802 DDU2(j)=DU2orig(j+1)
2811 write (iout,*) "DL2,DL1,DML,DDU1,DDU2"
2812 write(iout,'(10f10.5)') (DL2(k),k=3,dimen-1)
2813 write(iout,'(10f10.5)') (DL1(k),k=2,dimen-1)
2814 write(iout,'(10f10.5)') (DML(k),k=1,dimen-1)
2815 write(iout,'(10f10.5)') (DDU1(k),k=1,dimen-2)
2816 write(iout,'(10f10.5)') (DDU2(k),k=1,dimen-3)
2819 if (imark.gt.2) rs(imark-2)=-DU2orig(imark-2)
2820 if (imark.gt.1) rs(imark-1)=-DU1orig(imark-1)
2821 if (imark.lt.dimen) rs(imark)=-DU1orig(imark)
2822 if (imark.lt.dimen-1) rs(imark+1)=-DU2orig(imark)
2826 ! write (iout,*) "Vector rs"
2828 ! write (iout,*) j,rs(j)
2831 call FDIAG(dimen-1,DL2,DL1,DML,DDU1,DDU2,rs,xsolv,mark)
2838 sumx=-geigen(i)*xsolv(j)
2840 sumx=sumx+matold(j,k)*xsolv(k)
2843 sumx=sumx+matold(j,k)*xsolv(k-1)
2845 write(iout,'(i5,3f10.5)') j,sumx,rsold(j),sumx-rsold(j)
2848 sumx=-geigen(i)*xsolv(j-1)
2850 sumx=sumx+matold(j,k)*xsolv(k)
2853 sumx=sumx+matold(j,k)*xsolv(k-1)
2855 write(iout,'(i5,3f10.5)') j-1,sumx,rsold(j-1),sumx-rsold(j-1)
2859 "Solution of equations system",i," for eigenvalue",geigen(i)
2861 write(iout,'(i5,f10.5)') j,xsolv(j)
2864 do j=dimen-1,imark,-1
2869 write (iout,*) "Un-normalized eigenvector",i," for eigenvalue",geigen(i)
2871 write(iout,'(i5,f10.5)') j,xsolv(j)
2874 ! Normalize ith eigenvector
2877 sumx=sumx+xsolv(j)**2
2881 xsolv(j)=xsolv(j)/sumx
2884 write (iout,*) "Eigenvector",i," for eigenvalue",geigen(i)
2886 write(iout,'(i5,3f10.5)') j,xsolv(j)
2889 ! All done at this point for eigenvector i, exit loop
2897 write (iout,*) "Unable to find eigenvector",i
2900 ! write (iout,*) "i=",i
2902 ! write (iout,*) "k=",k
2905 ! write(iout,*) "ind",ind," ind1",3*(i-1)+k
2906 d_t_work(ind)=d_t_work(ind) &
2907 +xsolv(j)*d_t_work_new(3*(i-1)+k)
2910 enddo ! i (loop over the eigenvectors)
2913 write (iout,*) "d_t_work"
2915 write (iout,"(i5,f10.5)") i,d_t_work(i)
2920 ! if (itype(i,1).eq.10) then
2922 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
2923 iti=iabs(itype(i,mnum))
2924 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2925 if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
2926 .or.(mnum.eq.5)) then
2933 ! write (iout,*) "k",k," ii+k",ii+k," ii+j+k",ii+j+k,"EK1",Ek1
2934 Ek1=Ek1+0.5d0*mp(mnum)*((d_t_work(ii+j+k)+d_t_work_new(ii+k))/2)**2+&
2935 0.5d0*0.25d0*IP(mnum)*(d_t_work(ii+j+k)-d_t_work_new(ii+k))**2
2938 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2939 .and.(mnum.ne.5)) ii=ii+3
2940 write (iout,*) "i",i," itype",itype(i,mnum)," mass",msc(itype(i,mnum),mnum)
2941 write (iout,*) "ii",ii
2944 write (iout,*) "k",k," ii",ii,"EK1",EK1
2945 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2947 Ek1=Ek1+0.5d0*Isc(iabs(itype(i,mnum),mnum))*(d_t_work(ii)-d_t_work(ii-3))**2
2948 Ek1=Ek1+0.5d0*msc(iabs(itype(i,mnum)),mnum)*d_t_work(ii)**2
2950 write (iout,*) "i",i," ii",ii
2952 write (iout,*) "Ek from d_t_work",Ek1
2953 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
2969 d_t(k,j)=d_t_work(ind)
2973 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2976 d_t(k,j+nres)=d_t_work(ind)
2982 write (iout,*) "Random velocities in the Calpha,SC space"
2984 write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3)
2987 write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3)
2994 ! if (itype(i,1).eq.10) then
2996 if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
2999 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
3003 d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
3004 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
3009 write (iout,*)"Random velocities in the virtual-bond-vector space"
3011 write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3)
3014 write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3)
3017 write (iout,*) "Ek from d_t_work",Ek1
3018 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
3026 d_t_work(ind)=d_t_work(ind) &
3027 +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
3029 ! write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
3033 ! Transfer to the d_t vector
3035 d_t(j,0)=d_t_work(j)
3041 d_t(j,i)=d_t_work(ind)
3046 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3047 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3048 .and.(mnum.ne.5)) then
3051 d_t(j,i+nres)=d_t_work(ind)
3057 ! write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",&
3058 ! 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
3062 end subroutine random_vel
3063 !-----------------------------------------------------------------------------
3065 subroutine sd_verlet_p_setup
3066 ! Sets up the parameters of stochastic Verlet algorithm
3067 ! implicit real*8 (a-h,o-z)
3068 ! include 'DIMENSIONS'
3069 use control, only: tcpu
3074 ! include 'COMMON.CONTROL'
3075 ! include 'COMMON.VAR'
3076 ! include 'COMMON.MD'
3078 ! include 'COMMON.LANGEVIN'
3080 ! include 'COMMON.LANGEVIN.lang0'
3082 ! include 'COMMON.CHAIN'
3083 ! include 'COMMON.DERIV'
3084 ! include 'COMMON.GEO'
3085 ! include 'COMMON.LOCAL'
3086 ! include 'COMMON.INTERACT'
3087 ! include 'COMMON.IOUNITS'
3088 ! include 'COMMON.NAMES'
3089 ! include 'COMMON.TIME1'
3090 real(kind=8),dimension(6*nres) :: emgdt !(MAXRES6) maxres6=6*maxres
3091 real(kind=8) :: pterm,vterm,rho,rhoc,vsig
3092 real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
3093 prand_vec,vrand_vec1,vrand_vec2 !(MAXRES6) maxres6=6*maxres
3094 logical :: lprn = .false.
3095 real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
3096 real(kind=8) :: ktm,gdt,egdt,gdt2,gdt3,gdt4,gdt5,gdt6,gdt7,gdt8,&
3098 integer :: i,maxres2
3105 ! AL 8/17/04 Code adapted from tinker
3107 ! Get the frictional and random terms for stochastic dynamics in the
3108 ! eigenspace of mass-scaled UNRES friction matrix
3112 gdt = fricgam(i) * d_time
3114 ! Stochastic dynamics reduces to simple MD for zero friction
3116 if (gdt .le. zero) then
3117 pfric_vec(i) = 1.0d0
3118 vfric_vec(i) = d_time
3119 afric_vec(i) = 0.5d0 * d_time * d_time
3120 prand_vec(i) = 0.0d0
3121 vrand_vec1(i) = 0.0d0
3122 vrand_vec2(i) = 0.0d0
3124 ! Analytical expressions when friction coefficient is large
3127 if (gdt .ge. gdt_radius) then
3130 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
3131 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
3132 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
3133 vterm = 1.0d0 - egdt**2
3134 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
3136 ! Use series expansions when friction coefficient is small
3147 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0 &
3148 - gdt5/120.0d0 + gdt6/720.0d0 &
3149 - gdt7/5040.0d0 + gdt8/40320.0d0 &
3150 - gdt9/362880.0d0) / fricgam(i)**2
3151 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
3152 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
3153 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0 &
3154 + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0 &
3155 + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0 &
3156 + 127.0d0*gdt9/90720.0d0
3157 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0 &
3158 - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0 &
3159 - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0 &
3160 - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
3161 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0 &
3162 - 17.0d0*gdt2/1280.0d0 &
3163 + 17.0d0*gdt3/6144.0d0 &
3164 + 40967.0d0*gdt4/34406400.0d0 &
3165 - 57203.0d0*gdt5/275251200.0d0 &
3166 - 1429487.0d0*gdt6/13212057600.0d0)
3169 ! Compute the scaling factors of random terms for the nonzero friction case
3171 ktm = 0.5d0*d_time/fricgam(i)
3172 psig = dsqrt(ktm*pterm) / fricgam(i)
3173 vsig = dsqrt(ktm*vterm)
3174 rhoc = dsqrt(1.0d0 - rho*rho)
3176 vrand_vec1(i) = vsig * rho
3177 vrand_vec2(i) = vsig * rhoc
3182 "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",&
3185 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),&
3186 afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
3190 ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
3193 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
3194 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
3195 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
3196 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
3197 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
3198 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
3201 t_sdsetup=t_sdsetup+MPI_Wtime()
3203 t_sdsetup=t_sdsetup+tcpu()-tt0
3206 end subroutine sd_verlet_p_setup
3207 !-----------------------------------------------------------------------------
3208 subroutine eigtransf1(n,ndim,ab,d,c)
3212 real(kind=8) :: ab(ndim,ndim,n),c(ndim,n),d(ndim)
3218 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
3223 end subroutine eigtransf1
3224 !-----------------------------------------------------------------------------
3225 subroutine eigtransf(n,ndim,a,b,d,c)
3229 real(kind=8) :: a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
3235 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
3240 end subroutine eigtransf
3241 !-----------------------------------------------------------------------------
3242 subroutine sd_verlet1
3244 ! Applying stochastic velocity Verlet algorithm - step 1 to velocities
3246 ! implicit real*8 (a-h,o-z)
3247 ! include 'DIMENSIONS'
3248 ! include 'COMMON.CONTROL'
3249 ! include 'COMMON.VAR'
3250 ! include 'COMMON.MD'
3252 ! include 'COMMON.LANGEVIN'
3254 ! include 'COMMON.LANGEVIN.lang0'
3256 ! include 'COMMON.CHAIN'
3257 ! include 'COMMON.DERIV'
3258 ! include 'COMMON.GEO'
3259 ! include 'COMMON.LOCAL'
3260 ! include 'COMMON.INTERACT'
3261 ! include 'COMMON.IOUNITS'
3262 ! include 'COMMON.NAMES'
3263 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
3264 !el common /stochcalc/ stochforcvec
3265 logical :: lprn = .false.
3266 real(kind=8) :: ddt1,ddt2
3267 integer :: i,j,ind,inres
3269 ! write (iout,*) "dc_old"
3271 ! write (iout,'(i5,3f10.5,5x,3f10.5)')
3272 ! & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3275 dc_work(j)=dc_old(j,0)
3276 d_t_work(j)=d_t_old(j,0)
3277 d_a_work(j)=d_a_old(j,0)
3282 dc_work(ind+j)=dc_old(j,i)
3283 d_t_work(ind+j)=d_t_old(j,i)
3284 d_a_work(ind+j)=d_a_old(j,i)
3290 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3291 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3292 .and.(mnum.ne.5)) then
3294 dc_work(ind+j)=dc_old(j,i+nres)
3295 d_t_work(ind+j)=d_t_old(j,i+nres)
3296 d_a_work(ind+j)=d_a_old(j,i+nres)
3304 "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",&
3308 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
3309 vfric_mat(i,j),afric_mat(i,j),&
3310 prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3318 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) &
3319 +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3320 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3321 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3323 d_t_work_new(i)=ddt1+0.5d0*ddt2
3324 d_t_work(i)=ddt1+ddt2
3329 d_t(j,0)=d_t_work(j)
3334 dc(j,i)=dc_work(ind+j)
3335 d_t(j,i)=d_t_work(ind+j)
3341 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3342 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3343 .and.(mnum.ne.5)) then
3346 dc(j,inres)=dc_work(ind+j)
3347 d_t(j,inres)=d_t_work(ind+j)
3353 end subroutine sd_verlet1
3354 !-----------------------------------------------------------------------------
3355 subroutine sd_verlet2
3357 ! Calculating the adjusted velocities for accelerations
3359 ! implicit real*8 (a-h,o-z)
3360 ! include 'DIMENSIONS'
3361 ! include 'COMMON.CONTROL'
3362 ! include 'COMMON.VAR'
3363 ! include 'COMMON.MD'
3365 ! include 'COMMON.LANGEVIN'
3367 ! include 'COMMON.LANGEVIN.lang0'
3369 ! include 'COMMON.CHAIN'
3370 ! include 'COMMON.DERIV'
3371 ! include 'COMMON.GEO'
3372 ! include 'COMMON.LOCAL'
3373 ! include 'COMMON.INTERACT'
3374 ! include 'COMMON.IOUNITS'
3375 ! include 'COMMON.NAMES'
3376 !el real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV !(MAXRES6) maxres6=6*maxres
3377 real(kind=8),dimension(6*nres) :: stochforcvecV !(MAXRES6) maxres6=6*maxres
3378 !el common /stochcalc/ stochforcvec
3380 real(kind=8) :: ddt1,ddt2
3381 integer :: i,j,ind,inres
3382 ! Compute the stochastic forces which contribute to velocity change
3384 call stochastic_force(stochforcvecV)
3391 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3392 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+ &
3393 vrand_mat2(i,j)*stochforcvecV(j)
3395 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
3399 d_t(j,0)=d_t_work(j)
3404 d_t(j,i)=d_t_work(ind+j)
3410 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3411 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3412 .and.(mnum.ne.5)) then
3415 d_t(j,inres)=d_t_work(ind+j)
3421 end subroutine sd_verlet2
3422 !-----------------------------------------------------------------------------
3423 subroutine sd_verlet_ciccotti_setup
3425 ! Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
3427 ! implicit real*8 (a-h,o-z)
3428 ! include 'DIMENSIONS'
3429 use control, only: tcpu
3434 ! include 'COMMON.CONTROL'
3435 ! include 'COMMON.VAR'
3436 ! include 'COMMON.MD'
3438 ! include 'COMMON.LANGEVIN'
3440 ! include 'COMMON.LANGEVIN.lang0'
3442 ! include 'COMMON.CHAIN'
3443 ! include 'COMMON.DERIV'
3444 ! include 'COMMON.GEO'
3445 ! include 'COMMON.LOCAL'
3446 ! include 'COMMON.INTERACT'
3447 ! include 'COMMON.IOUNITS'
3448 ! include 'COMMON.NAMES'
3449 ! include 'COMMON.TIME1'
3450 real(kind=8),dimension(6*nres) :: emgdt !(MAXRES6) maxres6=6*maxres
3451 real(kind=8) :: pterm,vterm,rho,rhoc,vsig
3452 real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
3453 prand_vec,vrand_vec1,vrand_vec2 !(MAXRES6) maxres6=6*maxres
3454 logical :: lprn = .false.
3455 real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
3456 real(kind=8) :: ktm,gdt,egdt,tt0
3457 integer :: i,maxres2
3464 ! AL 8/17/04 Code adapted from tinker
3466 ! Get the frictional and random terms for stochastic dynamics in the
3467 ! eigenspace of mass-scaled UNRES friction matrix
3471 write (iout,*) "i",i," fricgam",fricgam(i)
3472 gdt = fricgam(i) * d_time
3474 ! Stochastic dynamics reduces to simple MD for zero friction
3476 if (gdt .le. zero) then
3477 pfric_vec(i) = 1.0d0
3478 vfric_vec(i) = d_time
3479 afric_vec(i) = 0.5d0*d_time*d_time
3480 prand_vec(i) = afric_vec(i)
3481 vrand_vec2(i) = vfric_vec(i)
3483 ! Analytical expressions when friction coefficient is large
3488 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
3489 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
3490 prand_vec(i) = afric_vec(i)
3491 vrand_vec2(i) = vfric_vec(i)
3493 ! Compute the scaling factors of random terms for the nonzero friction case
3495 ! ktm = 0.5d0*d_time/fricgam(i)
3496 ! psig = dsqrt(ktm*pterm) / fricgam(i)
3497 ! vsig = dsqrt(ktm*vterm)
3498 ! prand_vec(i) = psig*afric_vec(i)
3499 ! vrand_vec2(i) = vsig*vfric_vec(i)
3504 "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",&
3507 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),&
3508 afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
3512 ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
3514 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
3515 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
3516 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
3517 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
3518 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
3520 t_sdsetup=t_sdsetup+MPI_Wtime()
3522 t_sdsetup=t_sdsetup+tcpu()-tt0
3525 end subroutine sd_verlet_ciccotti_setup
3526 !-----------------------------------------------------------------------------
3527 subroutine sd_verlet1_ciccotti
3529 ! Applying stochastic velocity Verlet algorithm - step 1 to velocities
3530 ! implicit real*8 (a-h,o-z)
3532 ! include 'DIMENSIONS'
3536 ! include 'COMMON.CONTROL'
3537 ! include 'COMMON.VAR'
3538 ! include 'COMMON.MD'
3540 ! include 'COMMON.LANGEVIN'
3542 ! include 'COMMON.LANGEVIN.lang0'
3544 ! include 'COMMON.CHAIN'
3545 ! include 'COMMON.DERIV'
3546 ! include 'COMMON.GEO'
3547 ! include 'COMMON.LOCAL'
3548 ! include 'COMMON.INTERACT'
3549 ! include 'COMMON.IOUNITS'
3550 ! include 'COMMON.NAMES'
3551 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
3552 !el common /stochcalc/ stochforcvec
3553 logical :: lprn = .false.
3554 real(kind=8) :: ddt1,ddt2
3555 integer :: i,j,ind,inres
3556 ! write (iout,*) "dc_old"
3558 ! write (iout,'(i5,3f10.5,5x,3f10.5)')
3559 ! & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3562 dc_work(j)=dc_old(j,0)
3563 d_t_work(j)=d_t_old(j,0)
3564 d_a_work(j)=d_a_old(j,0)
3569 dc_work(ind+j)=dc_old(j,i)
3570 d_t_work(ind+j)=d_t_old(j,i)
3571 d_a_work(ind+j)=d_a_old(j,i)
3576 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3578 dc_work(ind+j)=dc_old(j,i+nres)
3579 d_t_work(ind+j)=d_t_old(j,i+nres)
3580 d_a_work(ind+j)=d_a_old(j,i+nres)
3589 "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",&
3593 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
3594 vfric_mat(i,j),afric_mat(i,j),&
3595 prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3603 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) &
3604 +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3605 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3606 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3608 d_t_work_new(i)=ddt1+0.5d0*ddt2
3609 d_t_work(i)=ddt1+ddt2
3614 d_t(j,0)=d_t_work(j)
3619 dc(j,i)=dc_work(ind+j)
3620 d_t(j,i)=d_t_work(ind+j)
3626 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3627 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3628 .and.(mnum.ne.5)) then
3631 dc(j,inres)=dc_work(ind+j)
3632 d_t(j,inres)=d_t_work(ind+j)
3638 end subroutine sd_verlet1_ciccotti
3639 !-----------------------------------------------------------------------------
3640 subroutine sd_verlet2_ciccotti
3642 ! Calculating the adjusted velocities for accelerations
3644 ! implicit real*8 (a-h,o-z)
3645 ! include 'DIMENSIONS'
3646 ! include 'COMMON.CONTROL'
3647 ! include 'COMMON.VAR'
3648 ! include 'COMMON.MD'
3650 ! include 'COMMON.LANGEVIN'
3652 ! include 'COMMON.LANGEVIN.lang0'
3654 ! include 'COMMON.CHAIN'
3655 ! include 'COMMON.DERIV'
3656 ! include 'COMMON.GEO'
3657 ! include 'COMMON.LOCAL'
3658 ! include 'COMMON.INTERACT'
3659 ! include 'COMMON.IOUNITS'
3660 ! include 'COMMON.NAMES'
3661 !el real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV !(MAXRES6) maxres6=6*maxres
3662 real(kind=8),dimension(6*nres) :: stochforcvecV !(MAXRES6) maxres6=6*maxres
3663 !el common /stochcalc/ stochforcvec
3664 real(kind=8) :: ddt1,ddt2
3665 integer :: i,j,ind,inres
3667 ! Compute the stochastic forces which contribute to velocity change
3669 call stochastic_force(stochforcvecV)
3676 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3677 ! ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
3678 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
3680 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
3684 d_t(j,0)=d_t_work(j)
3689 d_t(j,i)=d_t_work(ind+j)
3695 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3697 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3700 d_t(j,inres)=d_t_work(ind+j)
3706 end subroutine sd_verlet2_ciccotti
3708 !-----------------------------------------------------------------------------
3710 !-----------------------------------------------------------------------------
3711 subroutine inertia_tensor
3713 ! Calculating the intertia tensor for the entire protein in order to
3714 ! remove the perpendicular components of velocity matrix which cause
3715 ! the molecule to rotate.
3718 ! implicit real*8 (a-h,o-z)
3719 ! include 'DIMENSIONS'
3720 ! include 'COMMON.CONTROL'
3721 ! include 'COMMON.VAR'
3722 ! include 'COMMON.MD'
3723 ! include 'COMMON.CHAIN'
3724 ! include 'COMMON.DERIV'
3725 ! include 'COMMON.GEO'
3726 ! include 'COMMON.LOCAL'
3727 ! include 'COMMON.INTERACT'
3728 ! include 'COMMON.IOUNITS'
3729 ! include 'COMMON.NAMES'
3731 real(kind=8),dimension(3,3) :: Im,Imcp,eigvec,Id
3732 real(kind=8),dimension(3) :: pr,eigval,L,vp,vrot
3733 real(kind=8) :: M_SC,mag,mag2,M_PEP
3734 real(kind=8),dimension(3,0:nres) :: vpp !(3,0:MAXRES)
3735 real(kind=8),dimension(3) :: vs_p,pp,incr,v
3736 real(kind=8),dimension(3,3) :: pr1,pr2
3738 !el common /gucio/ cm
3739 integer :: iti,inres,i,j,k,mnum
3750 ! calculating the center of the mass of the protein
3754 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
3755 if (itype(i,mnum).eq.ntyp1_molec(mnum)) cycle
3756 M_PEP=M_PEP+mp(mnum)
3758 cm(j)=cm(j)+(c(j,i)+0.5d0*dc(j,i))*mp(mnum)
3767 if (mnum.eq.5) cycle
3768 iti=iabs(itype(i,mnum))
3769 M_SC=M_SC+msc(iabs(iti),mnum)
3772 cm(j)=cm(j)+msc(iabs(iti),mnum)*c(j,inres)
3776 cm(j)=cm(j)/(M_SC+M_PEP)
3781 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
3783 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
3785 Im(1,1)=Im(1,1)+mp(mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
3786 Im(1,2)=Im(1,2)-mp(mnum)*pr(1)*pr(2)
3787 Im(1,3)=Im(1,3)-mp(mnum)*pr(1)*pr(3)
3788 Im(2,3)=Im(2,3)-mp(mnum)*pr(2)*pr(3)
3789 Im(2,2)=Im(2,2)+mp(mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
3790 Im(3,3)=Im(3,3)+mp(mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
3795 iti=iabs(itype(i,mnum))
3796 if (mnum.eq.5) cycle
3799 pr(j)=c(j,inres)-cm(j)
3801 Im(1,1)=Im(1,1)+msc(iabs(iti),mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
3802 Im(1,2)=Im(1,2)-msc(iabs(iti),mnum)*pr(1)*pr(2)
3803 Im(1,3)=Im(1,3)-msc(iabs(iti),mnum)*pr(1)*pr(3)
3804 Im(2,3)=Im(2,3)-msc(iabs(iti),mnum)*pr(2)*pr(3)
3805 Im(2,2)=Im(2,2)+msc(iabs(iti),mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
3806 Im(3,3)=Im(3,3)+msc(iabs(iti),mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
3811 Im(1,1)=Im(1,1)+Ip(mnum)*(1-dc_norm(1,i)*dc_norm(1,i))* &
3812 vbld(i+1)*vbld(i+1)*0.25d0
3813 Im(1,2)=Im(1,2)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(2,i))* &
3814 vbld(i+1)*vbld(i+1)*0.25d0
3815 Im(1,3)=Im(1,3)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(3,i))* &
3816 vbld(i+1)*vbld(i+1)*0.25d0
3817 Im(2,3)=Im(2,3)+Ip(mnum)*(-dc_norm(2,i)*dc_norm(3,i))* &
3818 vbld(i+1)*vbld(i+1)*0.25d0
3819 Im(2,2)=Im(2,2)+Ip(mnum)*(1-dc_norm(2,i)*dc_norm(2,i))* &
3820 vbld(i+1)*vbld(i+1)*0.25d0
3821 Im(3,3)=Im(3,3)+Ip(mnum)*(1-dc_norm(3,i)*dc_norm(3,i))* &
3822 vbld(i+1)*vbld(i+1)*0.25d0
3828 ! if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
3829 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3830 .and.(mnum.ne.5)) then
3831 iti=iabs(itype(i,mnum))
3833 Im(1,1)=Im(1,1)+Isc(iti,mnum)*(1-dc_norm(1,inres)* &
3834 dc_norm(1,inres))*vbld(inres)*vbld(inres)
3835 Im(1,2)=Im(1,2)-Isc(iti,mnum)*(dc_norm(1,inres)* &
3836 dc_norm(2,inres))*vbld(inres)*vbld(inres)
3837 Im(1,3)=Im(1,3)-Isc(iti,mnum)*(dc_norm(1,inres)* &
3838 dc_norm(3,inres))*vbld(inres)*vbld(inres)
3839 Im(2,3)=Im(2,3)-Isc(iti,mnum)*(dc_norm(2,inres)* &
3840 dc_norm(3,inres))*vbld(inres)*vbld(inres)
3841 Im(2,2)=Im(2,2)+Isc(iti,mnum)*(1-dc_norm(2,inres)* &
3842 dc_norm(2,inres))*vbld(inres)*vbld(inres)
3843 Im(3,3)=Im(3,3)+Isc(iti,mnum)*(1-dc_norm(3,inres)* &
3844 dc_norm(3,inres))*vbld(inres)*vbld(inres)
3849 ! write(iout,*) "The angular momentum before adjustment:"
3850 ! write(iout,*) (L(j),j=1,3)
3856 ! Copying the Im matrix for the djacob subroutine
3864 ! Finding the eigenvectors and eignvalues of the inertia tensor
3865 call djacob(3,3,10000,1.0d-10,Imcp,eigvec,eigval)
3866 ! write (iout,*) "Eigenvalues & Eigenvectors"
3867 ! write (iout,'(5x,3f10.5)') (eigval(i),i=1,3)
3870 ! write (iout,'(i5,3f10.5)') i,(eigvec(i,j),j=1,3)
3872 ! Constructing the diagonalized matrix
3874 if (dabs(eigval(i)).gt.1.0d-15) then
3875 Id(i,i)=1.0d0/eigval(i)
3882 Imcp(i,j)=eigvec(j,i)
3888 pr1(i,j)=pr1(i,j)+Id(i,k)*Imcp(k,j)
3895 pr2(i,j)=pr2(i,j)+eigvec(i,k)*pr1(k,j)
3899 ! Calculating the total rotational velocity of the molecule
3902 vrot(i)=vrot(i)+pr2(i,j)*L(j)
3905 ! Resetting the velocities
3907 call vecpr(vrot(1),dc(1,i),vp)
3909 d_t(j,i)=d_t(j,i)-vp(j)
3914 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3915 .and.(mnum.ne.5)) then
3916 ! if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
3917 ! if(itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3919 call vecpr(vrot(1),dc(1,inres),vp)
3921 d_t(j,inres)=d_t(j,inres)-vp(j)
3926 ! write(iout,*) "The angular momentum after adjustment:"
3927 ! write(iout,*) (L(j),j=1,3)
3930 end subroutine inertia_tensor
3931 !-----------------------------------------------------------------------------
3932 subroutine angmom(cm,L)
3935 ! implicit real*8 (a-h,o-z)
3936 ! include 'DIMENSIONS'
3937 ! include 'COMMON.CONTROL'
3938 ! include 'COMMON.VAR'
3939 ! include 'COMMON.MD'
3940 ! include 'COMMON.CHAIN'
3941 ! include 'COMMON.DERIV'
3942 ! include 'COMMON.GEO'
3943 ! include 'COMMON.LOCAL'
3944 ! include 'COMMON.INTERACT'
3945 ! include 'COMMON.IOUNITS'
3946 ! include 'COMMON.NAMES'
3947 real(kind=8) :: mscab
3948 real(kind=8),dimension(3) :: L,cm,pr,vp,vrot,incr,v,pp
3949 integer :: iti,inres,i,j,mnum
3950 ! Calculate the angular momentum
3959 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
3961 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
3964 v(j)=incr(j)+0.5d0*d_t(j,i)
3967 incr(j)=incr(j)+d_t(j,i)
3969 call vecpr(pr(1),v(1),vp)
3971 L(j)=L(j)+mp(mnum)*vp(j)
3975 pp(j)=0.5d0*d_t(j,i)
3977 call vecpr(pr(1),pp(1),vp)
3979 L(j)=L(j)+Ip(mnum)*vp(j)
3987 iti=iabs(itype(i,mnum))
3995 pr(j)=c(j,inres)-cm(j)
3997 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3998 .and.(mnum.ne.5)) then
4000 v(j)=incr(j)+d_t(j,inres)
4007 call vecpr(pr(1),v(1),vp)
4008 ! write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),
4009 ! & " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
4011 L(j)=L(j)+mscab*vp(j)
4013 ! write (iout,*) "L",(l(j),j=1,3)
4014 if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4015 .and.(mnum.ne.5)) then
4017 v(j)=incr(j)+d_t(j,inres)
4019 call vecpr(dc(1,inres),d_t(1,inres),vp)
4021 L(j)=L(j)+Isc(iti,mnum)*vp(j)
4025 incr(j)=incr(j)+d_t(j,i)
4029 end subroutine angmom
4030 !-----------------------------------------------------------------------------
4031 subroutine vcm_vel(vcm)
4034 ! implicit real*8 (a-h,o-z)
4035 ! include 'DIMENSIONS'
4036 ! include 'COMMON.VAR'
4037 ! include 'COMMON.MD'
4038 ! include 'COMMON.CHAIN'
4039 ! include 'COMMON.DERIV'
4040 ! include 'COMMON.GEO'
4041 ! include 'COMMON.LOCAL'
4042 ! include 'COMMON.INTERACT'
4043 ! include 'COMMON.IOUNITS'
4044 real(kind=8),dimension(3) :: vcm,vv
4045 real(kind=8) :: summas,amas
4055 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
4057 summas=summas+mp(mnum)
4059 vcm(j)=vcm(j)+mp(mnum)*(vv(j)+0.5d0*d_t(j,i))
4063 amas=msc(iabs(itype(i,mnum)),mnum)
4068 if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4069 .and.(mnum.ne.5)) then
4071 vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres))
4075 vcm(j)=vcm(j)+amas*vv(j)
4079 vv(j)=vv(j)+d_t(j,i)
4082 ! write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas
4084 vcm(j)=vcm(j)/summas
4087 end subroutine vcm_vel
4088 !-----------------------------------------------------------------------------
4090 !-----------------------------------------------------------------------------
4092 ! RATTLE algorithm for velocity Verlet - step 1, UNRES
4096 ! implicit real*8 (a-h,o-z)
4097 ! include 'DIMENSIONS'
4099 ! include 'COMMON.CONTROL'
4100 ! include 'COMMON.VAR'
4101 ! include 'COMMON.MD'
4103 ! include 'COMMON.LANGEVIN'
4105 ! include 'COMMON.LANGEVIN.lang0'
4107 ! include 'COMMON.CHAIN'
4108 ! include 'COMMON.DERIV'
4109 ! include 'COMMON.GEO'
4110 ! include 'COMMON.LOCAL'
4111 ! include 'COMMON.INTERACT'
4112 ! include 'COMMON.IOUNITS'
4113 ! include 'COMMON.NAMES'
4114 ! include 'COMMON.TIME1'
4115 !el real(kind=8) :: gginv(2*nres,2*nres),&
4116 !el gdc(3,2*nres,2*nres)
4117 real(kind=8) :: dC_uncor(3,2*nres) !,&
4118 !el real(kind=8) :: Cmat(2*nres,2*nres)
4119 real(kind=8) :: x(2*nres),xcorr(3,2*nres) !maxres2=2*maxres
4120 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
4121 !el common /przechowalnia/ nbond
4122 integer :: max_rattle = 5
4123 logical :: lprn = .false., lprn1 = .false., not_done
4124 real(kind=8) :: tol_rattle = 1.0d-5
4126 integer :: ii,i,j,jj,l,ind,ind1,nres2
4129 !el /common/ przechowalnia
4131 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4132 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4133 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4135 if (lprn) write (iout,*) "RATTLE1"
4139 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4140 .and.(mnum.ne.5)) nbond=nbond+1
4142 ! Make a folded form of the Ginv-matrix
4155 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4160 if (itype(k,1).ne.10 .and. itype(k,mnum).ne.ntyp1_molec(mnum)&
4161 .and.(mnum.ne.5)) then
4165 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4173 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4184 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4188 if (itype(k,1).ne.10) then
4192 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4200 write (iout,*) "Matrix GGinv"
4201 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
4207 if (iter.gt.max_rattle) then
4208 write (iout,*) "Error - too many iterations in RATTLE."
4211 ! Calculate the matrix C = GG**(-1) dC_old o dC
4216 dC_uncor(j,ind1)=dC(j,i)
4220 if (itype(i,1).ne.10) then
4223 dC_uncor(j,ind1)=dC(j,i+nres)
4232 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
4236 if (itype(k,1).ne.10) then
4239 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
4244 ! Calculate deviations from standard virtual-bond lengths
4248 x(ind)=vbld(i+1)**2-vbl**2
4251 if (itype(i,1).ne.10) then
4253 x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4257 write (iout,*) "Coordinates and violations"
4259 write(iout,'(i5,3f10.5,5x,e15.5)') &
4260 i,(dC_uncor(j,i),j=1,3),x(i)
4262 write (iout,*) "Velocities and violations"
4266 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4267 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4271 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4272 .and.(mnum.ne.5)) then
4275 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4276 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4277 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4280 ! write (iout,*) "gdc"
4282 ! write (iout,*) "i",i
4284 ! write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4290 if (dabs(x(i)).gt.xmax) then
4294 if (xmax.lt.tol_rattle) then
4298 ! Calculate the matrix of the system of equations
4303 Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
4308 write (iout,*) "Matrix Cmat"
4309 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4311 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4312 ! Add constraint term to positions
4319 xx = xx+x(ii)*gdc(j,ind,ii)
4323 d_t_new(j,i)=d_t_new(j,i)-xx/d_time
4327 if (itype(i,1).ne.10) then
4332 xx = xx+x(ii)*gdc(j,ind,ii)
4335 dC(j,i+nres)=dC(j,i+nres)-xx
4336 d_t_new(j,i+nres)=d_t_new(j,i+nres)-xx/d_time
4340 ! Rebuild the chain using the new coordinates
4341 call chainbuild_cart
4343 write (iout,*) "New coordinates, Lagrange multipliers,",&
4344 " and differences between actual and standard bond lengths"
4348 xx=vbld(i+1)**2-vbl**2
4349 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4350 i,(dC(j,i),j=1,3),x(ind),xx
4354 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4357 xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4358 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4359 i,(dC(j,i+nres),j=1,3),x(ind),xx
4362 write (iout,*) "Velocities and violations"
4366 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4367 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4370 if (itype(i,1).ne.10) then
4372 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4373 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4374 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4381 10 write (iout,*) "Error - singularity in solving the system",&
4382 " of equations for Lagrange multipliers."
4386 "RATTLE inactive; use -DRATTLE switch at compile time."
4389 end subroutine rattle1
4390 !-----------------------------------------------------------------------------
4392 ! RATTLE algorithm for velocity Verlet - step 2, UNRES
4396 ! implicit real*8 (a-h,o-z)
4397 ! include 'DIMENSIONS'
4399 ! include 'COMMON.CONTROL'
4400 ! include 'COMMON.VAR'
4401 ! include 'COMMON.MD'
4403 ! include 'COMMON.LANGEVIN'
4405 ! include 'COMMON.LANGEVIN.lang0'
4407 ! include 'COMMON.CHAIN'
4408 ! include 'COMMON.DERIV'
4409 ! include 'COMMON.GEO'
4410 ! include 'COMMON.LOCAL'
4411 ! include 'COMMON.INTERACT'
4412 ! include 'COMMON.IOUNITS'
4413 ! include 'COMMON.NAMES'
4414 ! include 'COMMON.TIME1'
4415 !el real(kind=8) :: gginv(2*nres,2*nres),&
4416 !el gdc(3,2*nres,2*nres)
4417 real(kind=8) :: dC_uncor(3,2*nres) !,&
4418 !el Cmat(2*nres,2*nres)
4419 real(kind=8) :: x(2*nres) !maxres2=2*maxres
4420 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
4421 !el common /przechowalnia/ nbond
4422 integer :: max_rattle = 5
4423 logical :: lprn = .false., lprn1 = .false., not_done
4424 real(kind=8) :: tol_rattle = 1.0d-5
4428 !el /common/ przechowalnia
4429 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4430 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4431 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4433 if (lprn) write (iout,*) "RATTLE2"
4434 if (lprn) write (iout,*) "Velocity correction"
4435 ! Calculate the matrix G dC
4441 gdc(j,i,ind)=GGinv(i,ind)*dC(j,k)
4446 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4447 .and.(mnum.ne.5)) then
4450 gdc(j,i,ind)=GGinv(i,ind)*dC(j,k+nres)
4456 ! write (iout,*) "gdc"
4458 ! write (iout,*) "i",i
4460 ! write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4464 ! Calculate the matrix of the system of equations
4471 Cmat(ind,j)=Cmat(ind,j)+dC(k,i)*gdc(k,ind,j)
4477 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4478 .and.(mnum.ne.5)) then
4483 Cmat(ind,j)=Cmat(ind,j)+dC(k,i+nres)*gdc(k,ind,j)
4488 ! Calculate the scalar product dC o d_t_new
4492 x(ind)=scalar(d_t(1,i),dC(1,i))
4496 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4497 .and.(mnum.ne.5)) then
4499 x(ind)=scalar(d_t(1,i+nres),dC(1,i+nres))
4503 write (iout,*) "Velocities and violations"
4507 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4508 i,ind,(d_t(j,i),j=1,3),x(ind)
4512 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4513 .and.(mnum.ne.5)) then
4515 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4516 i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind)
4522 if (dabs(x(i)).gt.xmax) then
4526 if (xmax.lt.tol_rattle) then
4531 write (iout,*) "Matrix Cmat"
4532 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4534 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4535 ! Add constraint term to velocities
4542 xx = xx+x(ii)*gdc(j,ind,ii)
4544 d_t(j,i)=d_t(j,i)-xx
4549 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4550 .and.(mnum.ne.5)) then
4555 xx = xx+x(ii)*gdc(j,ind,ii)
4557 d_t(j,i+nres)=d_t(j,i+nres)-xx
4563 "New velocities, Lagrange multipliers violations"
4567 if (lprn) write (iout,'(2i5,3f10.5,5x,2e15.5)') &
4568 i,ind,(d_t(j,i),j=1,3),x(ind),scalar(d_t(1,i),dC(1,i))
4572 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4575 write (iout,'(2i5,3f10.5,5x,2e15.5)') &
4576 i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind),&
4577 scalar(d_t(1,i+nres),dC(1,i+nres))
4583 10 write (iout,*) "Error - singularity in solving the system",&
4584 " of equations for Lagrange multipliers."
4588 "RATTLE inactive; use -DRATTLE option at compile time."
4591 end subroutine rattle2
4592 !-----------------------------------------------------------------------------
4593 subroutine rattle_brown
4594 ! RATTLE/LINCS algorithm for Brownian dynamics, UNRES
4598 ! implicit real*8 (a-h,o-z)
4599 ! include 'DIMENSIONS'
4601 ! include 'COMMON.CONTROL'
4602 ! include 'COMMON.VAR'
4603 ! include 'COMMON.MD'
4605 ! include 'COMMON.LANGEVIN'
4607 ! include 'COMMON.LANGEVIN.lang0'
4609 ! include 'COMMON.CHAIN'
4610 ! include 'COMMON.DERIV'
4611 ! include 'COMMON.GEO'
4612 ! include 'COMMON.LOCAL'
4613 ! include 'COMMON.INTERACT'
4614 ! include 'COMMON.IOUNITS'
4615 ! include 'COMMON.NAMES'
4616 ! include 'COMMON.TIME1'
4617 !el real(kind=8) :: gginv(2*nres,2*nres),&
4618 !el gdc(3,2*nres,2*nres)
4619 real(kind=8) :: dC_uncor(3,2*nres) !,&
4620 !el real(kind=8) :: Cmat(2*nres,2*nres)
4621 real(kind=8) :: x(2*nres) !maxres2=2*maxres
4622 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
4623 !el common /przechowalnia/ nbond
4624 integer :: max_rattle = 5
4625 logical :: lprn = .false., lprn1 = .false., not_done
4626 real(kind=8) :: tol_rattle = 1.0d-5
4630 !el /common/ przechowalnia
4631 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4632 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4633 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4636 if (lprn) write (iout,*) "RATTLE_BROWN"
4639 if (itype(i,1).ne.10) nbond=nbond+1
4641 ! Make a folded form of the Ginv-matrix
4654 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4658 if (itype(k,1).ne.10) then
4662 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4669 if (itype(i,1).ne.10) then
4679 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4683 if (itype(k,1).ne.10) then
4687 if (j.eq.1 .and. l.eq.1)GGinv(ii,jj)=fricmat(ind,ind1)
4695 write (iout,*) "Matrix GGinv"
4696 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
4702 if (iter.gt.max_rattle) then
4703 write (iout,*) "Error - too many iterations in RATTLE."
4706 ! Calculate the matrix C = GG**(-1) dC_old o dC
4711 dC_uncor(j,ind1)=dC(j,i)
4715 if (itype(i,1).ne.10) then
4718 dC_uncor(j,ind1)=dC(j,i+nres)
4727 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
4731 if (itype(k,1).ne.10) then
4734 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
4739 ! Calculate deviations from standard virtual-bond lengths
4743 x(ind)=vbld(i+1)**2-vbl**2
4746 if (itype(i,1).ne.10) then
4748 x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4752 write (iout,*) "Coordinates and violations"
4754 write(iout,'(i5,3f10.5,5x,e15.5)') &
4755 i,(dC_uncor(j,i),j=1,3),x(i)
4757 write (iout,*) "Velocities and violations"
4761 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4762 i,ind,(d_t(j,i),j=1,3),scalar(d_t(1,i),dC_old(1,i))
4765 if (itype(i,1).ne.10) then
4767 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4768 i+nres,ind,(d_t(j,i+nres),j=1,3),&
4769 scalar(d_t(1,i+nres),dC_old(1,i+nres))
4772 write (iout,*) "gdc"
4774 write (iout,*) "i",i
4776 write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4782 if (dabs(x(i)).gt.xmax) then
4786 if (xmax.lt.tol_rattle) then
4790 ! Calculate the matrix of the system of equations
4795 Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
4800 write (iout,*) "Matrix Cmat"
4801 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4803 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4804 ! Add constraint term to positions
4811 xx = xx+x(ii)*gdc(j,ind,ii)
4814 d_t(j,i)=d_t(j,i)+xx/d_time
4819 if (itype(i,1).ne.10) then
4824 xx = xx+x(ii)*gdc(j,ind,ii)
4827 d_t(j,i+nres)=d_t(j,i+nres)+xx/d_time
4828 dC(j,i+nres)=dC(j,i+nres)+xx
4832 ! Rebuild the chain using the new coordinates
4833 call chainbuild_cart
4835 write (iout,*) "New coordinates, Lagrange multipliers,",&
4836 " and differences between actual and standard bond lengths"
4840 xx=vbld(i+1)**2-vbl**2
4841 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4842 i,(dC(j,i),j=1,3),x(ind),xx
4845 if (itype(i,1).ne.10) then
4847 xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4848 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4849 i,(dC(j,i+nres),j=1,3),x(ind),xx
4852 write (iout,*) "Velocities and violations"
4856 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4857 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4860 if (itype(i,1).ne.10) then
4862 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4863 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4864 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4871 10 write (iout,*) "Error - singularity in solving the system",&
4872 " of equations for Lagrange multipliers."
4876 "RATTLE inactive; use -DRATTLE option at compile time"
4879 end subroutine rattle_brown
4880 !-----------------------------------------------------------------------------
4882 !-----------------------------------------------------------------------------
4883 subroutine friction_force
4888 ! implicit real*8 (a-h,o-z)
4889 ! include 'DIMENSIONS'
4890 ! include 'COMMON.VAR'
4891 ! include 'COMMON.CHAIN'
4892 ! include 'COMMON.DERIV'
4893 ! include 'COMMON.GEO'
4894 ! include 'COMMON.LOCAL'
4895 ! include 'COMMON.INTERACT'
4896 ! include 'COMMON.MD'
4898 ! include 'COMMON.LANGEVIN'
4900 ! include 'COMMON.LANGEVIN.lang0'
4902 ! include 'COMMON.IOUNITS'
4903 !el real(kind=8),dimension(6*nres) :: gamvec !(MAXRES6) maxres6=6*maxres
4904 !el common /syfek/ gamvec
4905 real(kind=8) :: vv(3),vvtot(3,nres),v_work(6*nres) !,&
4906 !el ginvfric(2*nres,2*nres) !maxres2=2*maxres
4907 !el common /przechowalnia/ ginvfric
4909 logical :: lprn = .false., checkmode = .false.
4910 integer :: i,j,ind,k,nres2,nres6,mnum
4914 if(.not.allocated(gamvec)) allocate(gamvec(nres6)) !(MAXRES6)
4915 if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
4923 d_t_work(j)=d_t(j,0)
4928 d_t_work(ind+j)=d_t(j,i)
4934 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
4935 .and.(mnum.ne.5)) then
4937 d_t_work(ind+j)=d_t(j,i+nres)
4943 call fricmat_mult(d_t_work,fric_work)
4945 if (.not.checkmode) return
4948 write (iout,*) "d_t_work and fric_work"
4950 write (iout,'(i3,2e15.5)') i,d_t_work(i),fric_work(i)
4954 friction(j,0)=fric_work(j)
4959 friction(j,i)=fric_work(ind+j)
4965 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
4966 .and.(mnum.ne.5)) then
4967 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
4969 friction(j,i+nres)=fric_work(ind+j)
4975 write(iout,*) "Friction backbone"
4977 write(iout,'(i5,3e15.5,5x,3e15.5)') &
4978 i,(friction(j,i),j=1,3),(d_t(j,i),j=1,3)
4980 write(iout,*) "Friction side chain"
4982 write(iout,'(i5,3e15.5,5x,3e15.5)') &
4983 i,(friction(j,i+nres),j=1,3),(d_t(j,i+nres),j=1,3)
4992 vvtot(j,i)=vv(j)+0.5d0*d_t(j,i)
4993 vvtot(j,i+nres)=vv(j)+d_t(j,i+nres)
4994 vv(j)=vv(j)+d_t(j,i)
4997 write (iout,*) "vvtot backbone and sidechain"
4999 write (iout,'(i5,3e15.5,5x,3e15.5)') i,(vvtot(j,i),j=1,3),&
5000 (vvtot(j,i+nres),j=1,3)
5005 v_work(ind+j)=vvtot(j,i)
5011 v_work(ind+j)=vvtot(j,i+nres)
5015 write (iout,*) "v_work gamvec and site-based friction forces"
5017 write (iout,'(i5,3e15.5)') i,v_work(i),gamvec(i),&
5021 ! fric_work1(i)=0.0d0
5023 ! fric_work1(i)=fric_work1(i)-A(j,i)*gamvec(j)*v_work(j)
5026 ! write (iout,*) "fric_work and fric_work1"
5028 ! write (iout,'(i5,2e15.5)') i,fric_work(i),fric_work1(i)
5034 ginvfric(i,j)=ginvfric(i,j)+ginv(i,k)*fricmat(k,j)
5038 write (iout,*) "ginvfric"
5040 write (iout,'(i5,100f8.3)') i,(ginvfric(i,j),j=1,dimen)
5042 write (iout,*) "symmetry check"
5045 write (iout,*) i,j,ginvfric(i,j)-ginvfric(j,i)
5050 end subroutine friction_force
5051 !-----------------------------------------------------------------------------
5052 subroutine setup_fricmat
5056 use control_data, only:time_Bcast
5057 use control, only:tcpu
5059 ! implicit real*8 (a-h,o-z)
5063 real(kind=8) :: time00
5065 ! include 'DIMENSIONS'
5066 ! include 'COMMON.VAR'
5067 ! include 'COMMON.CHAIN'
5068 ! include 'COMMON.DERIV'
5069 ! include 'COMMON.GEO'
5070 ! include 'COMMON.LOCAL'
5071 ! include 'COMMON.INTERACT'
5072 ! include 'COMMON.MD'
5073 ! include 'COMMON.SETUP'
5074 ! include 'COMMON.TIME1'
5075 ! integer licznik /0/
5078 ! include 'COMMON.LANGEVIN'
5080 ! include 'COMMON.LANGEVIN.lang0'
5082 ! include 'COMMON.IOUNITS'
5084 integer :: i,j,ind,ind1,m
5085 logical :: lprn = .false.
5086 real(kind=8) :: dtdi !el ,gamvec(2*nres)
5087 !el real(kind=8),dimension(2*nres,2*nres) :: ginvfric,fcopy
5088 real(kind=8),dimension(2*nres,2*nres) :: fcopy
5089 !el real(kind=8),dimension(2*nres*(2*nres+1)/2) :: Ghalf !(mmaxres2) (mmaxres2=(maxres2*(maxres2+1)/2))
5090 !el common /syfek/ gamvec
5091 real(kind=8) :: work(8*2*nres)
5092 integer :: iwork(2*nres)
5093 !el common /przechowalnia/ ginvfric,Ghalf,fcopy
5094 integer :: ii,iti,k,l,nzero,nres2,nres6,ierr,mnum
5096 if (fg_rank.ne.king) goto 10
5101 if(.not.allocated(gamvec)) allocate(gamvec(nres2)) !(MAXRES2)
5102 if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
5103 !el if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
5104 !el allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
5105 if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !maxres2=2*maxres
5107 !el if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
5108 ! Zeroing out fricmat
5114 ! Load the friction coefficients corresponding to peptide groups
5119 gamvec(ind1)=gamp(mnum)
5121 ! Load the friction coefficients corresponding to side chains
5125 gamsc(ntyp1_molec(j),j)=1.0d0
5132 gamvec(ii)=gamsc(iabs(iti),mnum)
5134 if (surfarea) call sdarea(gamvec)
5136 ! write (iout,*) "Matrix A and vector gamma"
5138 ! write (iout,'(i2,$)') i
5140 ! write (iout,'(f4.1,$)') A(i,j)
5142 ! write (iout,'(f8.3)') gamvec(i)
5146 write (iout,*) "Vector gamvec"
5148 write (iout,'(i5,f10.5)') i, gamvec(i)
5152 ! The friction matrix
5157 dtdi=dtdi+A(j,k)*A(j,i)*gamvec(j)
5164 write (iout,'(//a)') "Matrix fricmat"
5165 call matout2(dimen,dimen,nres2,nres2,fricmat)
5167 if (lang.eq.2 .or. lang.eq.3) then
5168 ! Mass-scale the friction matrix if non-direct integration will be performed
5174 Ginvfric(i,j)=Ginvfric(i,j)+ &
5175 Gsqrm(i,k)*Gsqrm(l,j)*fricmat(k,l)
5180 ! Diagonalize the friction matrix
5185 Ghalf(ind)=Ginvfric(i,j)
5188 call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
5191 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
5192 " mass-scaled friction matrix"
5193 call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
5195 ! Precompute matrices for tinker stochastic integrator
5202 mt1(i,j)=mt1(i,j)+fricvec(k,i)*gsqrm(k,j)
5203 mt2(i,j)=mt2(i,j)+fricvec(k,i)*gsqrp(k,j)
5209 else if (lang.eq.4) then
5210 ! Diagonalize the friction matrix
5215 Ghalf(ind)=fricmat(i,j)
5218 call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
5221 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
5223 call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
5225 ! Determine the number of zero eigenvalues of the friction matrix
5226 nzero=max0(dimen-dimen1,0)
5227 ! do while (fricgam(nzero+1).le.1.0d-5 .and. nzero.lt.dimen)
5230 write (iout,*) "Number of zero eigenvalues:",nzero
5235 fricmat(i,j)=fricmat(i,j) &
5236 +fricvec(i,k)*fricvec(j,k)/fricgam(k)
5241 write (iout,'(//a)') "Generalized inverse of fricmat"
5242 call matout(dimen,dimen,nres6,nres6,fricmat)
5247 if (nfgtasks.gt.1) then
5248 if (fg_rank.eq.0) then
5249 ! The matching BROADCAST for fg processors is called in ERGASTULUM
5255 call MPI_Bcast(10,1,MPI_INTEGER,king,FG_COMM,IERROR)
5257 time_Bcast=time_Bcast+MPI_Wtime()-time00
5259 time_Bcast=time_Bcast+tcpu()-time00
5261 ! print *,"Processor",myrank,
5262 ! & " BROADCAST iorder in SETUP_FRICMAT"
5265 write (iout,*) "setup_fricmat licznik"!,licznik !sp
5271 ! Scatter the friction matrix
5272 call MPI_Scatterv(fricmat(1,1),nginv_counts(0),&
5273 nginv_start(0),MPI_DOUBLE_PRECISION,fcopy(1,1),&
5274 myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
5277 time_scatter=time_scatter+MPI_Wtime()-time00
5278 time_scatter_fmat=time_scatter_fmat+MPI_Wtime()-time00
5280 time_scatter=time_scatter+tcpu()-time00
5281 time_scatter_fmat=time_scatter_fmat+tcpu()-time00
5285 do j=1,2*my_ng_count
5286 fricmat(j,i)=fcopy(i,j)
5289 ! write (iout,*) "My chunk of fricmat"
5290 ! call MATOUT2(my_ng_count,dimen,maxres2,maxres2,fcopy)
5294 end subroutine setup_fricmat
5295 !-----------------------------------------------------------------------------
5296 subroutine stochastic_force(stochforcvec)
5299 use random, only:anorm_distr
5300 ! implicit real*8 (a-h,o-z)
5301 ! include 'DIMENSIONS'
5302 use control, only: tcpu
5307 ! include 'COMMON.VAR'
5308 ! include 'COMMON.CHAIN'
5309 ! include 'COMMON.DERIV'
5310 ! include 'COMMON.GEO'
5311 ! include 'COMMON.LOCAL'
5312 ! include 'COMMON.INTERACT'
5313 ! include 'COMMON.MD'
5314 ! include 'COMMON.TIME1'
5316 ! include 'COMMON.LANGEVIN'
5318 ! include 'COMMON.LANGEVIN.lang0'
5320 ! include 'COMMON.IOUNITS'
5322 real(kind=8) :: x,sig,lowb,highb
5323 real(kind=8) :: ff(3),force(3,0:2*nres),zeta2,lowb2
5324 real(kind=8) :: highb2,sig2,forcvec(6*nres),stochforcvec(6*nres)
5325 real(kind=8) :: time00
5326 logical :: lprn = .false.
5327 integer :: i,j,ind,mnum
5331 stochforc(j,i)=0.0d0
5341 ! Compute the stochastic forces acting on bodies. Store in force.
5347 force(j,i)=anorm_distr(x,sig,lowb,highb)
5355 force(j,i+nres)=anorm_distr(x,sig2,lowb2,highb2)
5359 time_fsample=time_fsample+MPI_Wtime()-time00
5361 time_fsample=time_fsample+tcpu()-time00
5363 ! Compute the stochastic forces acting on virtual-bond vectors.
5369 stochforc(j,i)=ff(j)+0.5d0*force(j,i)
5372 ff(j)=ff(j)+force(j,i)
5374 ! if (itype(i+1,1).ne.ntyp1) then
5376 if (itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
5378 stochforc(j,i)=stochforc(j,i)+force(j,i+nres+1)
5379 ff(j)=ff(j)+force(j,i+nres+1)
5384 stochforc(j,0)=ff(j)+force(j,nnt+nres)
5388 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5389 .and.(mnum.ne.5)) then
5390 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
5392 stochforc(j,i+nres)=force(j,i+nres)
5398 stochforcvec(j)=stochforc(j,0)
5403 stochforcvec(ind+j)=stochforc(j,i)
5409 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5410 .and.(mnum.ne.5)) then
5411 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
5413 stochforcvec(ind+j)=stochforc(j,i+nres)
5419 write (iout,*) "stochforcvec"
5421 write(iout,'(i5,e15.5)') i,stochforcvec(i)
5423 write(iout,*) "Stochastic forces backbone"
5425 write(iout,'(i5,3e15.5)') i,(stochforc(j,i),j=1,3)
5427 write(iout,*) "Stochastic forces side chain"
5429 write(iout,'(i5,3e15.5)') &
5430 i,(stochforc(j,i+nres),j=1,3)
5438 write (iout,*) i,ind
5440 forcvec(ind+j)=force(j,i)
5445 write (iout,*) i,ind
5447 forcvec(j+ind)=force(j,i+nres)
5452 write (iout,*) "forcvec"
5456 write (iout,'(2i3,2f10.5)') i,j,force(j,i),&
5463 write (iout,'(2i3,2f10.5)') i,j,force(j,i+nres),&
5472 end subroutine stochastic_force
5473 !-----------------------------------------------------------------------------
5474 subroutine sdarea(gamvec)
5476 ! Scale the friction coefficients according to solvent accessible surface areas
5477 ! Code adapted from TINKER
5481 ! implicit real*8 (a-h,o-z)
5482 ! include 'DIMENSIONS'
5483 ! include 'COMMON.CONTROL'
5484 ! include 'COMMON.VAR'
5485 ! include 'COMMON.MD'
5487 ! include 'COMMON.LANGEVIN'
5489 ! include 'COMMON.LANGEVIN.lang0'
5491 ! include 'COMMON.CHAIN'
5492 ! include 'COMMON.DERIV'
5493 ! include 'COMMON.GEO'
5494 ! include 'COMMON.LOCAL'
5495 ! include 'COMMON.INTERACT'
5496 ! include 'COMMON.IOUNITS'
5497 ! include 'COMMON.NAMES'
5498 real(kind=8),dimension(2*nres) :: radius,gamvec !(maxres2)
5499 real(kind=8),parameter :: twosix = 1.122462048309372981d0
5500 logical :: lprn = .false.
5501 real(kind=8) :: probe,area,ratio
5502 integer :: i,j,ind,iti,mnum
5504 ! determine new friction coefficients every few SD steps
5506 ! set the atomic radii to estimates of sigma values
5508 ! print *,"Entered sdarea"
5514 ! Load peptide group radii
5517 radius(i)=pstok(mnum)
5519 ! Load side chain radii
5523 radius(i+nres)=restok(iti,mnum)
5526 ! write (iout,*) "i",i," radius",radius(i)
5529 radius(i) = radius(i) / twosix
5530 if (radius(i) .ne. 0.0d0) radius(i) = radius(i) + probe
5533 ! scale atomic friction coefficients by accessible area
5535 if (lprn) write (iout,*) &
5536 "Original gammas, surface areas, scaling factors, new gammas, ",&
5537 "std's of stochastic forces"
5540 if (radius(i).gt.0.0d0) then
5541 call surfatom (i,area,radius)
5542 ratio = dmax1(area/(4.0d0*pi*radius(i)**2),1.0d-1)
5543 if (lprn) write (iout,'(i5,3f10.5,$)') &
5544 i,gamvec(ind+1),area,ratio
5547 gamvec(ind) = ratio * gamvec(ind)
5550 stdforcp(i)=stdfp(mnum)*dsqrt(gamvec(ind))
5551 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcp(i)
5555 if (radius(i+nres).gt.0.0d0) then
5556 call surfatom (i+nres,area,radius)
5557 ratio = dmax1(area/(4.0d0*pi*radius(i+nres)**2),1.0d-1)
5558 if (lprn) write (iout,'(i5,3f10.5,$)') &
5559 i,gamvec(ind+1),area,ratio
5562 gamvec(ind) = ratio * gamvec(ind)
5565 stdforcsc(i)=stdfsc(itype(i,mnum),mnum)*dsqrt(gamvec(ind))
5566 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcsc(i)
5571 end subroutine sdarea
5572 !-----------------------------------------------------------------------------
5574 !-----------------------------------------------------------------------------
5577 ! ###################################################
5578 ! ## COPYRIGHT (C) 1996 by Jay William Ponder ##
5579 ! ## All Rights Reserved ##
5580 ! ###################################################
5582 ! ################################################################
5584 ! ## subroutine surfatom -- exposed surface area of an atom ##
5586 ! ################################################################
5589 ! "surfatom" performs an analytical computation of the surface
5590 ! area of a specified atom; a simplified version of "surface"
5592 ! literature references:
5594 ! T. J. Richmond, "Solvent Accessible Surface Area and
5595 ! Excluded Volume in Proteins", Journal of Molecular Biology,
5598 ! L. Wesson and D. Eisenberg, "Atomic Solvation Parameters
5599 ! Applied to Molecular Dynamics of Proteins in Solution",
5600 ! Protein Science, 1, 227-235 (1992)
5602 ! variables and parameters:
5604 ! ir number of atom for which area is desired
5605 ! area accessible surface area of the atom
5606 ! radius radii of each of the individual atoms
5609 subroutine surfatom(ir,area,radius)
5611 ! implicit real*8 (a-h,o-z)
5612 ! include 'DIMENSIONS'
5614 ! include 'COMMON.GEO'
5615 ! include 'COMMON.IOUNITS'
5617 integer :: nsup,nstart_sup
5618 ! double precision c,dc,dc_old,d_c_work,xloc,xrot,dc_norm
5619 ! common /chain/ c(3,maxres2+2),dc(3,0:maxres2),dc_old(3,0:maxres2),
5620 ! & xloc(3,maxres),xrot(3,maxres),dc_norm(3,0:maxres2),
5621 ! & dc_work(MAXRES6),nres,nres0
5622 integer,parameter :: maxarc=300
5626 integer :: mi,ni,narc
5627 integer :: key(maxarc)
5628 integer :: intag(maxarc)
5629 integer :: intag1(maxarc)
5630 real(kind=8) :: area,arcsum
5631 real(kind=8) :: arclen,exang
5632 real(kind=8) :: delta,delta2
5633 real(kind=8) :: eps,rmove
5634 real(kind=8) :: xr,yr,zr
5635 real(kind=8) :: rr,rrsq
5636 real(kind=8) :: rplus,rminus
5637 real(kind=8) :: axx,axy,axz
5638 real(kind=8) :: ayx,ayy
5639 real(kind=8) :: azx,azy,azz
5640 real(kind=8) :: uxj,uyj,uzj
5641 real(kind=8) :: tx,ty,tz
5642 real(kind=8) :: txb,tyb,td
5643 real(kind=8) :: tr2,tr,txr,tyr
5644 real(kind=8) :: tk1,tk2
5645 real(kind=8) :: thec,the,t,tb
5646 real(kind=8) :: txk,tyk,tzk
5647 real(kind=8) :: t1,ti,tf,tt
5648 real(kind=8) :: txj,tyj,tzj
5649 real(kind=8) :: ccsq,cc,xysq
5650 real(kind=8) :: bsqk,bk,cosine
5651 real(kind=8) :: dsqj,gi,pix2
5652 real(kind=8) :: therk,dk,gk
5653 real(kind=8) :: risqk,rik
5654 real(kind=8) :: radius(2*nres) !(maxatm) (maxatm=maxres2)
5655 real(kind=8) :: ri(maxarc),risq(maxarc)
5656 real(kind=8) :: ux(maxarc),uy(maxarc),uz(maxarc)
5657 real(kind=8) :: xc(maxarc),yc(maxarc),zc(maxarc)
5658 real(kind=8) :: xc1(maxarc),yc1(maxarc),zc1(maxarc)
5659 real(kind=8) :: dsq(maxarc),bsq(maxarc)
5660 real(kind=8) :: dsq1(maxarc),bsq1(maxarc)
5661 real(kind=8) :: arci(maxarc),arcf(maxarc)
5662 real(kind=8) :: ex(maxarc),lt(maxarc),gr(maxarc)
5663 real(kind=8) :: b(maxarc),b1(maxarc),bg(maxarc)
5664 real(kind=8) :: kent(maxarc),kout(maxarc)
5665 real(kind=8) :: ther(maxarc)
5666 logical :: moved,top
5667 logical :: omit(maxarc)
5670 ! maxatm = 2*nres !maxres2 maxres2=2*maxres
5671 ! maxlight = 8*maxatm
5674 ! maxtors = 4*maxatm
5676 ! zero out the surface area for the sphere of interest
5679 ! write (2,*) "ir",ir," radius",radius(ir)
5680 if (radius(ir) .eq. 0.0d0) return
5682 ! set the overlap significance and connectivity shift
5686 delta2 = delta * delta
5691 ! store coordinates and radius of the sphere of interest
5699 ! initialize values of some counters and summations
5708 ! test each sphere to see if it overlaps the sphere of interest
5711 if (i.eq.ir .or. radius(i).eq.0.0d0) goto 30
5712 rplus = rr + radius(i)
5714 if (abs(tx) .ge. rplus) goto 30
5716 if (abs(ty) .ge. rplus) goto 30
5718 if (abs(tz) .ge. rplus) goto 30
5720 ! check for sphere overlap by testing distance against radii
5722 xysq = tx*tx + ty*ty
5723 if (xysq .lt. delta2) then
5730 if (rplus-cc .le. delta) goto 30
5731 rminus = rr - radius(i)
5733 ! check to see if sphere of interest is completely buried
5735 if (cc-abs(rminus) .le. delta) then
5736 if (rminus .le. 0.0d0) goto 170
5740 ! check for too many overlaps with sphere of interest
5742 if (io .ge. maxarc) then
5744 20 format (/,' SURFATOM -- Increase the Value of MAXARC')
5748 ! get overlap between current sphere and sphere of interest
5757 gr(io) = (ccsq+rplus*rminus) / (2.0d0*rr*b1(io))
5763 ! case where no other spheres overlap the sphere of interest
5766 area = 4.0d0 * pi * rrsq
5770 ! case where only one sphere overlaps the sphere of interest
5773 area = pix2 * (1.0d0 + gr(1))
5774 area = mod(area,4.0d0*pi) * rrsq
5778 ! case where many spheres intersect the sphere of interest;
5779 ! sort the intersecting spheres by their degree of overlap
5781 call sort2 (io,gr,key)
5784 intag(i) = intag1(k)
5793 ! get radius of each overlap circle on surface of the sphere
5798 risq(i) = rrsq - gi*gi
5799 ri(i) = sqrt(risq(i))
5800 ther(i) = 0.5d0*pi - asin(min(1.0d0,max(-1.0d0,gr(i))))
5803 ! find boundary of inaccessible area on sphere of interest
5806 if (.not. omit(k)) then
5813 ! check to see if J circle is intersecting K circle;
5814 ! get distance between circle centers and sum of radii
5817 if (omit(j)) goto 60
5818 cc = (txk*xc(j)+tyk*yc(j)+tzk*zc(j))/(bk*b(j))
5819 cc = acos(min(1.0d0,max(-1.0d0,cc)))
5820 td = therk + ther(j)
5822 ! check to see if circles enclose separate regions
5824 if (cc .ge. td) goto 60
5826 ! check for circle J completely inside circle K
5828 if (cc+ther(j) .lt. therk) goto 40
5830 ! check for circles that are essentially parallel
5832 if (cc .gt. delta) goto 50
5837 ! check to see if sphere of interest is completely buried
5840 if (pix2-cc .le. td) goto 170
5846 ! find T value of circle intersections
5849 if (omit(k)) goto 110
5864 ! rotation matrix elements
5876 if (.not. omit(j)) then
5881 ! rotate spheres so K vector colinear with z-axis
5883 uxj = txj*axx + tyj*axy - tzj*axz
5884 uyj = tyj*ayy - txj*ayx
5885 uzj = txj*azx + tyj*azy + tzj*azz
5886 cosine = min(1.0d0,max(-1.0d0,uzj/b(j)))
5887 if (acos(cosine) .lt. therk+ther(j)) then
5888 dsqj = uxj*uxj + uyj*uyj
5893 tr2 = risqk*dsqj - tb*tb
5899 ! get T values of intersection for K circle
5902 tb = min(1.0d0,max(-1.0d0,tb))
5904 if (tyb-txr .lt. 0.0d0) tk1 = pix2 - tk1
5906 tb = min(1.0d0,max(-1.0d0,tb))
5908 if (tyb+txr .lt. 0.0d0) tk2 = pix2 - tk2
5909 thec = (rrsq*uzj-gk*bg(j)) / (rik*ri(j)*b(j))
5910 if (abs(thec) .lt. 1.0d0) then
5912 else if (thec .ge. 1.0d0) then
5914 else if (thec .le. -1.0d0) then
5918 ! see if "tk1" is entry or exit point; check t=0 point;
5919 ! "ti" is exit point, "tf" is entry point
5921 cosine = min(1.0d0,max(-1.0d0, &
5922 (uzj*gk-uxj*rik)/(b(j)*rr)))
5923 if ((acos(cosine)-ther(j))*(tk2-tk1) .le. 0.0d0) then
5931 if (narc .ge. maxarc) then
5933 70 format (/,' SURFATOM -- Increase the Value',&
5937 if (tf .le. ti) then
5958 ! special case; K circle without intersections
5960 if (narc .le. 0) goto 90
5962 ! general case; sum up arclength and set connectivity code
5964 call sort2 (narc,arci,key)
5969 if (narc .gt. 1) then
5972 if (t .lt. arci(j)) then
5973 arcsum = arcsum + arci(j) - t
5974 exang = exang + ex(ni)
5976 if (jb .ge. maxarc) then
5978 80 format (/,' SURFATOM -- Increase the Value',&
5983 kent(jb) = maxarc*i + k
5985 kout(jb) = maxarc*k + i
5994 arcsum = arcsum + pix2 - t
5996 exang = exang + ex(ni)
5999 kent(jb) = maxarc*i + k
6001 kout(jb) = maxarc*k + i
6008 arclen = arclen + gr(k)*arcsum
6011 if (arclen .eq. 0.0d0) goto 170
6012 if (jb .eq. 0) goto 150
6014 ! find number of independent boundaries and check connectivity
6018 if (kout(k) .ne. 0) then
6025 if (m .eq. kent(ii)) then
6028 if (j .eq. jb) goto 150
6040 ! attempt to fix connectivity error by moving atom slightly
6044 140 format (/,' SURFATOM -- Connectivity Error at Atom',i6)
6053 ! compute the exposed surface area for the sphere of interest
6056 area = ib*pix2 + exang + arclen
6057 area = mod(area,4.0d0*pi) * rrsq
6059 ! attempt to fix negative area by moving atom slightly
6061 if (area .lt. 0.0d0) then
6064 160 format (/,' SURFATOM -- Negative Area at Atom',i6)
6075 end subroutine surfatom
6076 !----------------------------------------------------------------
6077 !----------------------------------------------------------------
6078 subroutine alloc_MD_arrays
6079 !EL Allocation of arrays used by MD module
6081 integer :: nres2,nres6
6084 !----------------------
6088 allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
6089 allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
6090 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
6091 allocate(fricvec(nres2,nres2))
6092 allocate(pfric_mat(nres2,nres2),vfric_mat(nres2,nres2))
6093 allocate(afric_mat(nres2,nres2),prand_mat(nres2,nres2))
6094 allocate(vrand_mat1(nres2,nres2),vrand_mat2(nres2,nres2)) !(MAXRES2,MAXRES2)
6095 allocate(pfric0_mat(nres2,nres2,0:maxflag_stoch))
6096 allocate(afric0_mat(nres2,nres2,0:maxflag_stoch))
6097 allocate(vfric0_mat(nres2,nres2,0:maxflag_stoch))
6098 allocate(prand0_mat(nres2,nres2,0:maxflag_stoch))
6099 allocate(vrand0_mat1(nres2,nres2,0:maxflag_stoch))
6100 allocate(vrand0_mat2(nres2,nres2,0:maxflag_stoch)) !(MAXRES2,MAXRES2,0:maxflag_stoch)
6101 allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
6103 allocate(mt1(nres2,nres2),mt2(nres2,nres2),mt3(nres2,nres2)) !(maxres2,maxres2)
6104 !----------------------
6106 ! commom.langevin.lang0
6108 allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
6109 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
6110 allocate(fricvec(nres2,nres2)) !(MAXRES2,MAXRES2)
6111 allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
6112 allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
6115 !el if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2))
6116 !----------------------
6117 ! commom.hairpin in CSA module
6118 !----------------------
6119 ! common.mce in MCM_MD module
6120 !----------------------
6122 ! common /mdgrad/ in module.energy
6123 ! common /back_constr/ in module.energy
6124 ! common /qmeas/ in module.energy
6127 allocate(potEcomp(0:n_ene+4)) !(0:n_ene+4)
6129 allocate(d_t(3,0:nres2),d_a(3,0:nres2),d_t_old(3,0:nres2)) !(3,0:MAXRES2)
6130 allocate(d_a_work(nres6)) !(6*MAXRES)
6132 allocate(DM(nres2),DU1(nres2),DU2(nres2))
6133 allocate(DMorig(nres2),DU1orig(nres2),DU2orig(nres2))
6135 allocate(Gmat(nres2,nres2),A(nres2,nres2))
6136 if(.not.allocated(Ginv)) allocate(Ginv(nres2,nres2)) !in control: ergastulum
6137 allocate(Gsqrp(nres2,nres2),Gsqrm(nres2,nres2),Gvec(nres2,nres2)) !(maxres2,maxres2)
6139 allocate(Geigen(nres2)) !(maxres2)
6140 if(.not.allocated(vtot)) allocate(vtot(nres2)) !(maxres2)
6141 ! common /inertia/ in io_conf: parmread
6142 ! real(kind=8),dimension(:),allocatable :: ISC,msc !(ntyp+1)
6143 ! common /langevin/in io read_MDpar
6144 ! real(kind=8),dimension(:),allocatable :: gamsc !(ntyp1)
6145 ! real(kind=8),dimension(:),allocatable :: stdfsc !(ntyp)
6146 ! in io_conf: parmread
6147 ! real(kind=8),dimension(:),allocatable :: restok !(ntyp+1)
6148 ! common /mdpmpi/ in control: ergastulum
6149 if(.not.allocated(ng_start)) allocate(ng_start(0:nfgtasks-1))
6150 if(.not.allocated(ng_counts)) allocate(ng_counts(0:nfgtasks-1))
6151 if(.not.allocated(nginv_counts)) allocate(nginv_counts(0:nfgtasks-1)) !(0:MaxProcs-1)
6152 if(.not.allocated(nginv_start)) allocate(nginv_start(0:nfgtasks)) !(0:MaxProcs)
6153 !----------------------
6154 ! common.muca in read_muca
6155 ! common /double_muca/
6156 ! real(kind=8) :: elow,ehigh,factor,hbin,factor_min
6157 ! real(kind=8),dimension(:),allocatable :: emuca,nemuca,&
6158 ! nemuca2,hist !(4*maxres)
6159 ! common /integer_muca/
6160 ! integer :: nmuca,imtime,muca_smooth
6162 ! real(kind=8),dimension(:),allocatable :: elowi,ehighi !(maxprocs)
6163 !----------------------
6165 ! common /mdgrad/ in module.energy
6166 ! common /back_constr/ in module.energy
6167 ! common /qmeas/ in module.energy
6171 allocate(d_t_work(nres6),d_t_work_new(nres6),d_af_work(nres6))
6172 allocate(d_as_work(nres6),kinetic_force(nres6)) !(MAXRES6)
6173 allocate(d_t_new(3,0:nres2),d_a_old(3,0:nres2),d_a_short(3,0:nres2)) !,d_a !(3,0:MAXRES2)
6174 allocate(stdforcp(nres),stdforcsc(nres)) !(MAXRES)
6175 !----------------------
6177 allocate(D_ban(nres6)) !(MAXRES6) maxres6=6*maxres
6178 ! common /stochcalc/ stochforcvec
6179 allocate(stochforcvec(nres6)) !(MAXRES6) maxres6=6*maxres
6180 !----------------------
6182 end subroutine alloc_MD_arrays
6183 !-----------------------------------------------------------------------------
6184 !-----------------------------------------------------------------------------