2 !-----------------------------------------------------------------------------
15 !-----------------------------------------------------------------------------
17 ! common /mdgrad/ in module.energy
18 ! common /back_constr/ in module.energy
19 ! common /qmeas/ in module.energy
23 real(kind=8),dimension(:),allocatable :: d_t_work,&
24 d_t_work_new,d_af_work,d_as_work,kinetic_force !(MAXRES6)
25 real(kind=8),dimension(:,:),allocatable :: d_t_new,&
26 d_a_old,d_a_short!,d_a !(3,0:MAXRES2)
27 ! real(kind=8),dimension(:),allocatable :: d_a_work !(6*MAXRES)
28 ! real(kind=8),dimension(:,:),allocatable :: Gmat,Ginv,A,&
29 ! Gsqrp,Gsqrm,Gvec !(maxres2,maxres2)
30 ! real(kind=8),dimension(:),allocatable :: Geigen !(maxres2)
31 ! integer :: dimen,dimen1,dimen3
32 ! integer :: lang,count_reset_moment,count_reset_vel
33 ! logical :: reset_moment,reset_vel,rattle,RESPA
36 ! real(kind=8) :: rwat,etawat,stdfp,cPoise
37 ! real(kind=8),dimension(:),allocatable :: gamsc !(ntyp1)
38 ! real(kind=8),dimension(:),allocatable :: stdfsc !(ntyp)
39 real(kind=8),dimension(:),allocatable :: stdforcp,stdforcsc !(MAXRES)
40 !-----------------------------------------------------------------------------
44 ! ###################################################
45 ! ## COPYRIGHT (C) 1992 by Jay William Ponder ##
46 ! ## All Rights Reserved ##
47 ! ###################################################
49 ! #############################################################
51 ! ## sizes.i -- parameter values to set array dimensions ##
53 ! #############################################################
56 ! "sizes.i" sets values for critical array dimensions used
57 ! throughout the software; these parameters will fix the size
58 ! of the largest systems that can be handled; values too large
59 ! for the computer's memory and/or swap space to accomodate
60 ! will result in poor performance or outright failure
62 ! parameter: maximum allowed number of:
64 ! maxatm atoms in the molecular system
65 ! maxval atoms directly bonded to an atom
66 ! maxgrp ! user-defined groups of atoms
67 ! maxtyp force field atom type definitions
68 ! maxclass force field atom class definitions
69 ! maxkey lines in the keyword file
70 ! maxrot bonds for torsional rotation
71 ! maxvar optimization variables (vector storage)
72 ! maxopt optimization variables (matrix storage)
73 ! maxhess off-diagonal Hessian elements
74 ! maxlight sites for method of lights neighbors
75 ! maxvib vibrational frequencies
76 ! maxgeo distance geometry points
77 ! maxcell unit cells in replicated crystal
78 ! maxring 3-, 4-, or 5-membered rings
79 ! maxfix geometric restraints
80 ! maxbio biopolymer atom definitions
81 ! maxres residues in the macromolecule
82 ! maxamino amino acid residue types
83 ! maxnuc nucleic acid residue types
84 ! maxbnd covalent bonds in molecular system
85 ! maxang bond angles in molecular system
86 ! maxtors torsional angles in molecular system
87 ! maxpi atoms in conjugated pisystem
88 ! maxpib covalent bonds involving pisystem
89 ! maxpit torsional angles involving pisystem
92 !el integer maxatm,maxval,maxgrp
93 !el integer maxtyp,maxclass,maxkey
94 !el integer maxrot,maxopt
95 !el integer maxhess,maxlight,maxvib
96 !el integer maxgeo,maxcell,maxring
97 !el integer maxfix,maxbio
98 !el integer maxamino,maxnuc,maxbnd
99 !el integer maxang,maxtors,maxpi
100 !el integer maxpib,maxpit
101 ! integer :: maxatm !=2*nres !maxres2 maxres2=2*maxres
102 ! integer,parameter :: maxval=8
103 ! integer,parameter :: maxgrp=1000
104 ! integer,parameter :: maxtyp=3000
105 ! integer,parameter :: maxclass=500
106 ! integer,parameter :: maxkey=10000
107 ! integer,parameter :: maxrot=1000
108 ! integer,parameter :: maxopt=1000
109 ! integer,parameter :: maxhess=1000000
110 ! integer :: maxlight !=8*maxatm
111 ! integer,parameter :: maxvib=1000
112 ! integer,parameter :: maxgeo=1000
113 ! integer,parameter :: maxcell=10000
114 ! integer,parameter :: maxring=10000
115 ! integer,parameter :: maxfix=10000
116 ! integer,parameter :: maxbio=10000
117 ! integer,parameter :: maxamino=31
118 ! integer,parameter :: maxnuc=12
119 ! integer :: maxbnd !=2*maxatm
120 ! integer :: maxang !=3*maxatm
121 ! integer :: maxtors !=4*maxatm
122 ! integer,parameter :: maxpi=100
123 ! integer,parameter :: maxpib=2*maxpi
124 ! integer,parameter :: maxpit=4*maxpi
125 !-----------------------------------------------------------------------------
126 ! Maximum number of seed
127 ! integer,parameter :: max_seed=1
128 !-----------------------------------------------------------------------------
129 real(kind=8),dimension(:),allocatable :: stochforcvec !(MAXRES6) maxres6=6*maxres
130 ! common /stochcalc/ stochforcvec
131 !-----------------------------------------------------------------------------
132 ! common /przechowalnia/ subroutines: rattle1,rattle2,rattle_brown
133 real(kind=8),dimension(:,:),allocatable :: GGinv !(2*nres,2*nres) maxres2=2*maxres
134 real(kind=8),dimension(:,:,:),allocatable :: gdc !(3,2*nres,2*nres) maxres2=2*maxres
135 real(kind=8),dimension(:,:),allocatable :: Cmat !(2*nres,2*nres) maxres2=2*maxres
136 !-----------------------------------------------------------------------------
137 ! common /syfek/ subroutines: friction_force,setup_fricmat
138 !el real(kind=8),dimension(:),allocatable :: gamvec !(MAXRES6) or (MAXRES2)
139 !-----------------------------------------------------------------------------
140 ! common /przechowalnia/ subroutines: friction_force,setup_fricmat
141 real(kind=8),dimension(:,:),allocatable :: ginvfric !(2*nres,2*nres) !maxres2=2*maxres
142 !-----------------------------------------------------------------------------
143 ! common /przechowalnia/ subroutine: setup_fricmat
144 real(kind=8),dimension(:,:),allocatable :: fcopy !(2*nres,2*nres)
145 !-----------------------------------------------------------------------------
148 !-----------------------------------------------------------------------------
150 !-----------------------------------------------------------------------------
152 !-----------------------------------------------------------------------------
153 subroutine brown_step(itime)
154 !------------------------------------------------
155 ! Perform a single Euler integration step of Brownian dynamics
156 !------------------------------------------------
157 ! implicit real*8 (a-h,o-z)
159 use control, only: tcpu
162 ! use io_conf, only:cartprint
163 ! include 'DIMENSIONS'
167 ! include 'COMMON.CONTROL'
168 ! include 'COMMON.VAR'
169 ! include 'COMMON.MD'
171 ! include 'COMMON.LANGEVIN'
173 ! include 'COMMON.LANGEVIN.lang0'
175 ! include 'COMMON.CHAIN'
176 ! include 'COMMON.DERIV'
177 ! include 'COMMON.GEO'
178 ! include 'COMMON.LOCAL'
179 ! include 'COMMON.INTERACT'
180 ! include 'COMMON.IOUNITS'
181 ! include 'COMMON.NAMES'
182 ! include 'COMMON.TIME1'
183 real(kind=8),dimension(6*nres) :: zapas !(MAXRES6) maxres6=6*maxres
184 integer :: rstcount !ilen,
186 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
187 real(kind=8),dimension(6*nres,2*nres) :: Bmat,GBmat,Tmat !(MAXRES6,MAXRES2) (maxres2=2*maxres,maxres6=6*maxres)
188 real(kind=8),dimension(2*nres,2*nres) :: Cmat_,Cinv !(maxres2,maxres2) maxres2=2*maxres
189 real(kind=8),dimension(6*nres,6*nres) :: Pmat !(maxres6,maxres6) maxres6=6*maxres
190 ! real(kind=8),dimension(:,:),allocatable :: Bmat,GBmat,Tmat !(MAXRES6,MAXRES2) (maxres2=2*maxres,maxres6=6*maxres)
191 ! real(kind=8),dimension(:,:),allocatable :: Cmat_,Cinv !(maxres2,maxres2) maxres2=2*maxres
192 ! real(kind=8),dimension(:,:),allocatable :: Pmat !(maxres6,maxres6) maxres6=6*maxres
193 real(kind=8),dimension(6*nres) :: Td !(maxres6) maxres6=6*maxres
194 real(kind=8),dimension(2*nres) :: ppvec !(maxres2) maxres2=2*maxres
195 !el common /stochcalc/ stochforcvec
196 !el real(kind=8),dimension(3) :: cm !el
197 !el common /gucio/ cm
199 logical :: lprn = .false.,lprn1 = .false.
200 integer :: maxiter = 5
201 real(kind=8) :: difftol = 1.0d-5
202 real(kind=8) :: xx,diffmax,blen2,diffbond,tt0
203 integer :: i,j,nbond,k,ind,ind1,iter
204 integer :: nres2,nres6
208 ! if (.not.allocated(Bmat)) allocate(Bmat(nres6,nres2))
209 ! if (.not.allocated(GBmat)) allocate (GBmat(nres6,nres2))
210 ! if (.not.allocated(Tmat)) allocate (Tmat(nres6,nres2))
211 ! if (.not.allocated(Cmat_)) allocate(Cmat_(nres2,nres2))
212 ! if (.not.allocated(Cinv)) allocate (Cinv(nres2,nres2))
213 ! if (.not.allocated(Pmat)) allocate(Pmat(6*nres,6*nres))
215 if (.not.allocated(stochforcvec)) allocate(stochforcvec(nres6)) !(MAXRES6) maxres6=6*maxres
219 if (itype(i,1).ne.10) nbond=nbond+1
223 write (iout,*) "Generalized inverse of fricmat"
224 call matout(dimen,dimen,nres6,nres6,fricmat)
236 Bmat(ind+j,ind1)=dC_norm(j,i)
241 if (itype(i,1).ne.10) then
244 Bmat(ind+j,ind1)=dC_norm(j,i+nres)
250 write (iout,*) "Matrix Bmat"
251 call MATOUT(nbond,dimen,nres6,nres6,Bmat)
257 GBmat(i,j)=GBmat(i,j)+fricmat(i,k)*Bmat(k,j)
262 write (iout,*) "Matrix GBmat"
263 call MATOUT(nbond,dimen,nres6,nres2,Gbmat)
269 Cmat_(i,j)=Cmat_(i,j)+Bmat(k,i)*GBmat(k,j)
274 write (iout,*) "Matrix Cmat"
275 call MATOUT(nbond,nbond,nres2,nres2,Cmat_)
277 call matinvert(nbond,nres2,Cmat_,Cinv,osob)
279 write (iout,*) "Matrix Cinv"
280 call MATOUT(nbond,nbond,nres2,nres2,Cinv)
286 Tmat(i,j)=Tmat(i,j)+GBmat(i,k)*Cinv(k,j)
291 write (iout,*) "Matrix Tmat"
292 call MATOUT(nbond,dimen,nres6,nres2,Tmat)
302 Pmat(i,j)=Pmat(i,j)-Tmat(i,k)*Bmat(j,k)
307 write (iout,*) "Matrix Pmat"
308 call MATOUT(dimen,dimen,nres6,nres6,Pmat)
315 Td(i)=Td(i)+vbl*Tmat(i,ind)
318 if (itype(k,1).ne.10) then
320 Td(i)=Td(i)+vbldsc0(1,itype(k,1))*Tmat(i,ind)
325 write (iout,*) "Vector Td"
327 write (iout,'(i5,f10.5)') i,Td(i)
330 call stochastic_force(stochforcvec)
332 write (iout,*) "stochforcvec"
334 write (iout,*) i,stochforcvec(i)
338 zapas(j)=-gcart(j,0)+stochforcvec(j)
340 dC_work(j)=dC_old(j,0)
346 zapas(ind)=-gcart(j,i)+stochforcvec(ind)
347 dC_work(ind)=dC_old(j,i)
351 if (itype(i,1).ne.10) then
354 zapas(ind)=-gxcart(j,i)+stochforcvec(ind)
355 dC_work(ind)=dC_old(j,i+nres)
361 write (iout,*) "Initial d_t_work"
363 write (iout,*) i,d_t_work(i)
370 d_t_work(i)=d_t_work(i)+fricmat(i,j)*zapas(j)
377 zapas(i)=zapas(i)+Pmat(i,j)*(dC_work(j)+d_t_work(j)*d_time)
381 write (iout,*) "Final d_t_work and zapas"
383 write (iout,*) i,d_t_work(i),zapas(i)
397 dc_work(ind+j)=dc(j,i)
403 d_t(j,i+nres)=d_t_work(ind+j)
404 dc(j,i+nres)=zapas(ind+j)
405 dc_work(ind+j)=dc(j,i+nres)
411 write (iout,*) "Before correction for rotational lengthening"
412 write (iout,*) "New coordinates",&
413 " and differences between actual and standard bond lengths"
418 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
422 if (itype(i,1).ne.10) then
424 xx=vbld(i+nres)-vbldsc0(1,itype(i,1))
425 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
426 i,(dC(j,i+nres),j=1,3),xx
430 ! Second correction (rotational lengthening)
436 blen2 = scalar(dc(1,i),dc(1,i))
437 ppvec(ind)=2*vbl**2-blen2
438 diffbond=dabs(vbl-dsqrt(blen2))
439 if (diffbond.gt.diffmax) diffmax=diffbond
440 if (ppvec(ind).gt.0.0d0) then
441 ppvec(ind)=dsqrt(ppvec(ind))
446 write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind)
450 if (itype(i,1).ne.10) then
452 blen2 = scalar(dc(1,i+nres),dc(1,i+nres))
453 ppvec(ind)=2*vbldsc0(1,itype(i,1))**2-blen2
454 diffbond=dabs(vbldsc0(1,itype(i,1))-dsqrt(blen2))
455 if (diffbond.gt.diffmax) diffmax=diffbond
456 if (ppvec(ind).gt.0.0d0) then
457 ppvec(ind)=dsqrt(ppvec(ind))
462 write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind)
466 if (lprn) write (iout,*) "iter",iter," diffmax",diffmax
467 if (diffmax.lt.difftol) goto 10
471 Td(i)=Td(i)+ppvec(j)*Tmat(i,j)
477 zapas(i)=zapas(i)+Pmat(i,j)*dc_work(j)
488 dc_work(ind+j)=zapas(ind+j)
493 if (itype(i,1).ne.10) then
495 dc(j,i+nres)=zapas(ind+j)
496 dc_work(ind+j)=zapas(ind+j)
501 ! Building the chain from the newly calculated coordinates
504 if (large.and. mod(itime,ntwe).eq.0) then
505 write (iout,*) "Cartesian and internal coordinates: step 1"
508 write (iout,'(a)') "Potential forces"
510 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(-gcart(j,i),j=1,3),&
513 write (iout,'(a)') "Stochastic forces"
515 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(stochforc(j,i),j=1,3),&
516 (stochforc(j,i+nres),j=1,3)
518 write (iout,'(a)') "Velocities"
520 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
521 (d_t(j,i+nres),j=1,3)
526 write (iout,*) "After correction for rotational lengthening"
527 write (iout,*) "New coordinates",&
528 " and differences between actual and standard bond lengths"
533 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
537 if (itype(i,1).ne.10) then
539 xx=vbld(i+nres)-vbldsc0(1,itype(i,1))
540 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
541 i,(dC(j,i+nres),j=1,3),xx
546 ! write (iout,*) "Too many attempts at correcting the bonds"
554 ! Calculate energy and forces
556 call etotal(potEcomp)
557 potE=potEcomp(0)-potEcomp(20)
561 ! Calculate the kinetic and total energy and the kinetic temperature
564 t_enegrad=t_enegrad+MPI_Wtime()-tt0
566 t_enegrad=t_enegrad+tcpu()-tt0
569 kinetic_T=2.0d0/(dimen*Rb)*EK
571 end subroutine brown_step
572 !-----------------------------------------------------------------------------
574 !-----------------------------------------------------------------------------
575 subroutine gauss(RO,AP,MT,M,N,*)
577 ! CALCULATES (RO**(-1))*AP BY GAUSS ELIMINATION
578 ! RO IS A SQUARE MATRIX
579 ! THE CALCULATED PRODUCT IS STORED IN AP
580 ! ABNORMAL EXIT IF RO IS SINGULAR
582 integer :: MT, M, N, M1,I,J,IM,&
584 real(kind=8) :: RO(MT,M),AP(MT,N),X,RM,PR,Y
590 if(dabs(X).le.1.0D-13) return 1
602 if(DABS(RO(J,I)).LE.RM) goto 2
616 if(dabs(X).le.1.0E-13) return 1
625 8 AP(J,K)=AP(J,K)-Y*AP(I,K)
627 9 RO(J,K)=RO(J,K)-Y*RO(I,K)
631 if(dabs(X).le.1.0E-13) return 1
641 15 X=X-AP(K,J)*RO(MI,K)
646 !-----------------------------------------------------------------------------
648 !-----------------------------------------------------------------------------
649 subroutine kinetic(KE_total)
650 !----------------------------------------------------------------
651 ! This subroutine calculates the total kinetic energy of the chain
652 !-----------------------------------------------------------------
654 ! implicit real*8 (a-h,o-z)
655 ! include 'DIMENSIONS'
656 ! include 'COMMON.VAR'
657 ! include 'COMMON.CHAIN'
658 ! include 'COMMON.DERIV'
659 ! include 'COMMON.GEO'
660 ! include 'COMMON.LOCAL'
661 ! include 'COMMON.INTERACT'
662 ! include 'COMMON.MD'
663 ! include 'COMMON.IOUNITS'
664 real(kind=8) :: KE_total,mscab
666 integer :: i,j,k,iti,mnum
667 real(kind=8) :: KEt_p,KEt_sc,KEr_p,KEr_sc,incr(3),&
670 write (iout,*) "Velocities, kietic"
672 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
673 (d_t(j,i+nres),j=1,3)
678 ! write (iout,*) "ISC",(isc(itype(i,1)),i=1,nres)
679 ! The translational part for peptide virtual bonds
685 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
686 ! write (iout,*) "Kinetic trp:",i,(incr(j),j=1,3)
688 v(j)=incr(j)+0.5d0*d_t(j,i)
690 vtot(i)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
691 KEt_p=KEt_p+mp(mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
693 incr(j)=incr(j)+d_t(j,i)
696 ! write(iout,*) 'KEt_p', KEt_p
697 ! The translational part for the side chain virtual bond
698 ! Only now we can initialize incr with zeros. It must be equal
699 ! to the velocities of the first Calpha.
705 iti=iabs(itype(i,mnum))
711 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
712 if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
713 .or.(mnum.eq.5)) then
719 v(j)=incr(j)+d_t(j,nres+i)
722 ! write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3)
723 ! write (iout,*) "i",i," msc",msc(iti)," v",(v(j),j=1,3)
724 KEt_sc=KEt_sc+mscab*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
725 vtot(i+nres)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
727 incr(j)=incr(j)+d_t(j,i)
731 ! write(iout,*) 'KEt_sc', KEt_sc
732 ! The part due to stretching and rotation of the peptide groups
736 ! write (iout,*) "i",i
737 ! write (iout,*) "i",i," mag1",mag1," mag2",mag2
741 ! write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3)
742 KEr_p=KEr_p+Ip(mnum)*(incr(1)*incr(1)+incr(2)*incr(2) &
746 ! write(iout,*) 'KEr_p', KEr_p
747 ! The rotational part of the side chain virtual bond
751 iti=iabs(itype(i,mnum))
752 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
753 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
754 .and.(mnum.ne.5)) then
756 incr(j)=d_t(j,nres+i)
758 ! write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3)
759 KEr_sc=KEr_sc+Isc(iti,mnum)*(incr(1)*incr(1)+incr(2)*incr(2)+ &
763 ! The total kinetic energy
765 ! write(iout,*) 'KEr_sc', KEr_sc
766 KE_total=0.5d0*(KEt_p+KEt_sc+0.25d0*KEr_p+KEr_sc)
767 ! write (iout,*) "KE_total",KE_total
769 end subroutine kinetic
770 !-----------------------------------------------------------------------------
772 !-----------------------------------------------------------------------------
774 !------------------------------------------------
775 ! The driver for molecular dynamics subroutines
776 !------------------------------------------------
779 use control, only:tcpu,ovrtim
780 ! use io_comm, only:ilen
782 use compare, only:secondary2,hairpin
783 use io, only:cartout,statout
784 ! implicit real*8 (a-h,o-z)
785 ! include 'DIMENSIONS'
788 integer :: IERROR,ERRCODE
790 ! include 'COMMON.SETUP'
791 ! include 'COMMON.CONTROL'
792 ! include 'COMMON.VAR'
793 ! include 'COMMON.MD'
795 ! include 'COMMON.LANGEVIN'
797 ! include 'COMMON.LANGEVIN.lang0'
799 ! include 'COMMON.CHAIN'
800 ! include 'COMMON.DERIV'
801 ! include 'COMMON.GEO'
802 ! include 'COMMON.LOCAL'
803 ! include 'COMMON.INTERACT'
804 ! include 'COMMON.IOUNITS'
805 ! include 'COMMON.NAMES'
806 ! include 'COMMON.TIME1'
807 ! include 'COMMON.HAIRPIN'
808 real(kind=8),dimension(3) :: L,vcm
810 real(kind=8),dimension(6*nres) :: v_work,v_transf !(maxres6) maxres6=6*maxres
812 integer :: rstcount !ilen,
814 character(len=50) :: tytul
815 !el common /gucio/ cm
816 integer :: itime,i,j,nharp
817 integer,dimension(4,nres/3) :: iharp !(4,nres/3)(4,maxres/3)
819 real(kind=8) :: tt0,scalfac
825 print *,"MY tmpdir",tmpdir,ilen(tmpdir)
826 if (ilen(tmpdir).gt.0) &
827 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_" &
828 //liczba(:ilen(liczba))//'.rst')
830 if (ilen(tmpdir).gt.0) &
831 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"//'.rst')
838 write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
844 print *,"just befor setup matix",nres
845 ! Determine the inverse of the inertia matrix.
846 call setup_MD_matrices
848 print *,"AFTER SETUP MATRICES"
850 print *,"AFTER INIT MD"
853 t_MDsetup = MPI_Wtime()-tt0
855 t_MDsetup = tcpu()-tt0
858 ! Entering the MD loop
864 if (lang.eq.2 .or. lang.eq.3) then
868 call sd_verlet_p_setup
870 call sd_verlet_ciccotti_setup
874 pfric0_mat(i,j,0)=pfric_mat(i,j)
875 afric0_mat(i,j,0)=afric_mat(i,j)
876 vfric0_mat(i,j,0)=vfric_mat(i,j)
877 prand0_mat(i,j,0)=prand_mat(i,j)
878 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
879 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
884 flag_stoch(i)=.false.
888 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
890 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
894 else if (lang.eq.1 .or. lang.eq.4) then
895 print *,"before setup_fricmat"
897 print *,"after setup_fricmat"
900 t_langsetup=MPI_Wtime()-tt0
903 t_langsetup=tcpu()-tt0
906 do itime=1,n_timestep
908 if (large.and. mod(itime,ntwe).eq.0) &
909 write (iout,*) "itime",itime
911 if (lang.gt.0 .and. surfarea .and. &
912 mod(itime,reset_fricmat).eq.0) then
913 if (lang.eq.2 .or. lang.eq.3) then
917 call sd_verlet_p_setup
919 call sd_verlet_ciccotti_setup
923 pfric0_mat(i,j,0)=pfric_mat(i,j)
924 afric0_mat(i,j,0)=afric_mat(i,j)
925 vfric0_mat(i,j,0)=vfric_mat(i,j)
926 prand0_mat(i,j,0)=prand_mat(i,j)
927 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
928 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
933 flag_stoch(i)=.false.
936 else if (lang.eq.1 .or. lang.eq.4) then
939 write (iout,'(a,i10)') &
940 "Friction matrix reset based on surface area, itime",itime
942 if (reset_vel .and. tbf .and. lang.eq.0 &
943 .and. mod(itime,count_reset_vel).eq.0) then
945 write(iout,'(a,f20.2)') &
946 "Velocities reset to random values, time",totT
949 d_t_old(j,i)=d_t(j,i)
953 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
957 d_t(j,0)=d_t(j,0)-vcm(j)
960 kinetic_T=2.0d0/(dimen3*Rb)*EK
961 scalfac=dsqrt(T_bath/kinetic_T)
962 write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
965 d_t_old(j,i)=scalfac*d_t(j,i)
971 ! Time-reversible RESPA algorithm
972 ! (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
973 call RESPA_step(itime)
975 ! Variable time step algorithm.
976 call velverlet_step(itime)
980 call brown_step(itime)
982 print *,"Brown dynamics not here!"
984 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
990 if (mod(itime,ntwe).eq.0) then
1002 v_work(ind)=d_t(j,i)
1007 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1.and.mnum.ne.5) then
1010 v_work(ind)=d_t(j,i+nres)
1015 write (66,'(80f10.5)') &
1016 ((d_t(j,i),j=1,3),i=0,nres-1),((d_t(j,i+nres),j=1,3),i=1,nres)
1020 v_transf(i)=v_transf(i)+gvec(j,i)*v_work(j)
1022 v_transf(i)= v_transf(i)*dsqrt(geigen(i))
1024 write (67,'(80f10.5)') (v_transf(i),i=1,ind)
1027 if (mod(itime,ntwx).eq.0) then
1029 write (tytul,'("time",f8.2)') totT
1031 call hairpin(.true.,nharp,iharp)
1032 call secondary2(.true.)
1033 call pdbout(potE,tytul,ipdb)
1038 if (rstcount.eq.1000.or.itime.eq.n_timestep) then
1039 open(irest2,file=rest2name,status='unknown')
1040 write(irest2,*) totT,EK,potE,totE,t_bath
1042 ! AL 4/17/17: Now writing d_t(0,:) too
1044 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
1046 ! AL 4/17/17: Now writing d_c(0,:) too
1048 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
1056 t_MD=MPI_Wtime()-tt0
1060 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))') &
1062 'MD calculations setup:',t_MDsetup,&
1063 'Energy & gradient evaluation:',t_enegrad,&
1064 'Stochastic MD setup:',t_langsetup,&
1065 'Stochastic MD step setup:',t_sdsetup,&
1067 write (iout,'(/28(1h=),a25,27(1h=))') &
1068 ' End of MD calculation '
1070 write (iout,*) "time for etotal",t_etotal," elong",t_elong,&
1072 write (iout,*) "time_fric",time_fric," time_stoch",time_stoch,&
1073 " time_fricmatmult",time_fricmatmult," time_fsample ",&
1078 !-----------------------------------------------------------------------------
1079 subroutine velverlet_step(itime)
1080 !-------------------------------------------------------------------------------
1081 ! Perform a single velocity Verlet step; the time step can be rescaled if
1082 ! increments in accelerations exceed the threshold
1083 !-------------------------------------------------------------------------------
1084 ! implicit real*8 (a-h,o-z)
1085 ! include 'DIMENSIONS'
1087 use control, only:tcpu
1091 integer :: ierror,ierrcode
1092 real(kind=8) :: errcode
1094 ! include 'COMMON.SETUP'
1095 ! include 'COMMON.VAR'
1096 ! include 'COMMON.MD'
1098 ! include 'COMMON.LANGEVIN'
1100 ! include 'COMMON.LANGEVIN.lang0'
1102 ! include 'COMMON.CHAIN'
1103 ! include 'COMMON.DERIV'
1104 ! include 'COMMON.GEO'
1105 ! include 'COMMON.LOCAL'
1106 ! include 'COMMON.INTERACT'
1107 ! include 'COMMON.IOUNITS'
1108 ! include 'COMMON.NAMES'
1109 ! include 'COMMON.TIME1'
1110 ! include 'COMMON.MUCA'
1111 real(kind=8),dimension(3) :: vcm,incr
1112 real(kind=8),dimension(3) :: L
1113 integer :: count,rstcount !ilen,
1115 character(len=50) :: tytul
1116 integer :: maxcount_scale = 20
1117 !el common /gucio/ cm
1118 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1119 !el common /stochcalc/ stochforcvec
1120 integer :: itime,icount_scale,itime_scal,i,j,ifac_time
1122 real(kind=8) :: epdrift,tt0,fac_time
1124 if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres)) !(MAXRES6) maxres6=6*maxres
1130 else if (lang.eq.2 .or. lang.eq.3) then
1132 call stochastic_force(stochforcvec)
1135 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
1137 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1144 icount_scale=icount_scale+1
1145 if (icount_scale.gt.maxcount_scale) then
1147 "ERROR: too many attempts at scaling down the time step. ",&
1148 "amax=",amax,"epdrift=",epdrift,&
1149 "damax=",damax,"edriftmax=",edriftmax,&
1153 call MPI_Abort(MPI_COMM_WORLD,IERROR,IERRCODE)
1157 ! First step of the velocity Verlet algorithm
1162 else if (lang.eq.3) then
1164 call sd_verlet1_ciccotti
1166 else if (lang.eq.1) then
1171 ! Build the chain from the newly calculated coordinates
1172 call chainbuild_cart
1173 if (rattle) call rattle1
1175 if (large.and. mod(itime,ntwe).eq.0) then
1176 write (iout,*) "Cartesian and internal coordinates: step 1"
1181 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),&
1182 (dc(j,i+nres),j=1,3)
1184 write (iout,*) "Accelerations"
1186 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1187 (d_a(j,i+nres),j=1,3)
1189 write (iout,*) "Velocities, step 1"
1191 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1192 (d_t(j,i+nres),j=1,3)
1201 ! Calculate energy and forces
1203 call etotal(potEcomp)
1204 ! AL 4/17/17: Reduce the steps if NaNs occurred.
1205 if (potEcomp(0).gt.0.99e20 .or. isnan(potEcomp(0)).gt.0) then
1210 if (large.and. mod(itime,ntwe).eq.0) &
1211 call enerprint(potEcomp)
1214 t_etotal=t_etotal+MPI_Wtime()-tt0
1216 t_etotal=t_etotal+tcpu()-tt0
1219 potE=potEcomp(0)-potEcomp(20)
1221 ! Get the new accelerations
1224 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1226 t_enegrad=t_enegrad+tcpu()-tt0
1228 ! Determine maximum acceleration and scale down the timestep if needed
1230 amax=amax/(itime_scal+1)**2
1231 call predict_edrift(epdrift)
1232 if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then
1233 ! Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
1235 ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax)) &
1237 itime_scal=itime_scal+ifac_time
1238 ! fac_time=dmin1(damax/amax,0.5d0)
1239 fac_time=0.5d0**ifac_time
1240 d_time=d_time*fac_time
1241 if (lang.eq.2 .or. lang.eq.3) then
1243 ! write (iout,*) "Calling sd_verlet_setup: 1"
1244 ! Rescale the stochastic forces and recalculate or restore
1245 ! the matrices of tinker integrator
1246 if (itime_scal.gt.maxflag_stoch) then
1247 if (large) write (iout,'(a,i5,a)') &
1248 "Calculate matrices for stochastic step;",&
1249 " itime_scal ",itime_scal
1251 call sd_verlet_p_setup
1253 call sd_verlet_ciccotti_setup
1255 write (iout,'(2a,i3,a,i3,1h.)') &
1256 "Warning: cannot store matrices for stochastic",&
1257 " integration because the index",itime_scal,&
1258 " is greater than",maxflag_stoch
1259 write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",&
1260 " integration Langevin algorithm for better efficiency."
1261 else if (flag_stoch(itime_scal)) then
1262 if (large) write (iout,'(a,i5,a,l1)') &
1263 "Restore matrices for stochastic step; itime_scal ",&
1264 itime_scal," flag ",flag_stoch(itime_scal)
1267 pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
1268 afric_mat(i,j)=afric0_mat(i,j,itime_scal)
1269 vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
1270 prand_mat(i,j)=prand0_mat(i,j,itime_scal)
1271 vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
1272 vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
1276 if (large) write (iout,'(2a,i5,a,l1)') &
1277 "Calculate & store matrices for stochastic step;",&
1278 " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
1280 call sd_verlet_p_setup
1282 call sd_verlet_ciccotti_setup
1284 flag_stoch(ifac_time)=.true.
1287 pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
1288 afric0_mat(i,j,itime_scal)=afric_mat(i,j)
1289 vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
1290 prand0_mat(i,j,itime_scal)=prand_mat(i,j)
1291 vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
1292 vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
1296 fac_time=1.0d0/dsqrt(fac_time)
1298 stochforcvec(i)=fac_time*stochforcvec(i)
1301 else if (lang.eq.1) then
1302 ! Rescale the accelerations due to stochastic forces
1303 fac_time=1.0d0/dsqrt(fac_time)
1305 d_as_work(i)=d_as_work(i)*fac_time
1308 if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)') &
1309 "itime",itime," Timestep scaled down to ",&
1310 d_time," ifac_time",ifac_time," itime_scal",itime_scal
1312 ! Second step of the velocity Verlet algorithm
1317 else if (lang.eq.3) then
1319 call sd_verlet2_ciccotti
1321 else if (lang.eq.1) then
1326 if (rattle) call rattle2
1329 if (d_time.ne.d_time0) then
1332 if (lang.eq.2 .or. lang.eq.3) then
1333 if (large) write (iout,'(a)') &
1334 "Restore original matrices for stochastic step"
1335 ! write (iout,*) "Calling sd_verlet_setup: 2"
1336 ! Restore the matrices of tinker integrator if the time step has been restored
1339 pfric_mat(i,j)=pfric0_mat(i,j,0)
1340 afric_mat(i,j)=afric0_mat(i,j,0)
1341 vfric_mat(i,j)=vfric0_mat(i,j,0)
1342 prand_mat(i,j)=prand0_mat(i,j,0)
1343 vrand_mat1(i,j)=vrand0_mat1(i,j,0)
1344 vrand_mat2(i,j)=vrand0_mat2(i,j,0)
1353 ! Calculate the kinetic and the total energy and the kinetic temperature
1357 ! call kinetic1(EK1)
1358 ! write (iout,*) "step",itime," EK",EK," EK1",EK1
1360 ! Couple the system to Berendsen bath if needed
1361 if (tbf .and. lang.eq.0) then
1364 kinetic_T=2.0d0/(dimen3*Rb)*EK
1365 ! Backup the coordinates, velocities, and accelerations
1369 d_t_old(j,i)=d_t(j,i)
1370 d_a_old(j,i)=d_a(j,i)
1374 if (mod(itime,ntwe).eq.0 .and. large) then
1375 write (iout,*) "Velocities, step 2"
1377 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1378 (d_t(j,i+nres),j=1,3)
1383 end subroutine velverlet_step
1384 !-----------------------------------------------------------------------------
1385 subroutine RESPA_step(itime)
1386 !-------------------------------------------------------------------------------
1387 ! Perform a single RESPA step.
1388 !-------------------------------------------------------------------------------
1389 ! implicit real*8 (a-h,o-z)
1390 ! include 'DIMENSIONS'
1394 use control, only:tcpu
1396 ! use io_conf, only:cartprint
1399 integer :: IERROR,ERRCODE
1401 ! include 'COMMON.SETUP'
1402 ! include 'COMMON.CONTROL'
1403 ! include 'COMMON.VAR'
1404 ! include 'COMMON.MD'
1406 ! include 'COMMON.LANGEVIN'
1408 ! include 'COMMON.LANGEVIN.lang0'
1410 ! include 'COMMON.CHAIN'
1411 ! include 'COMMON.DERIV'
1412 ! include 'COMMON.GEO'
1413 ! include 'COMMON.LOCAL'
1414 ! include 'COMMON.INTERACT'
1415 ! include 'COMMON.IOUNITS'
1416 ! include 'COMMON.NAMES'
1417 ! include 'COMMON.TIME1'
1418 real(kind=8),dimension(0:n_ene) :: energia_short,energia_long
1419 real(kind=8),dimension(3) :: L,vcm,incr
1420 real(kind=8),dimension(3,0:2*nres) :: dc_old0,d_t_old0,d_a_old0 !(3,0:maxres2) maxres2=2*maxres
1421 logical :: PRINT_AMTS_MSG = .false.
1422 integer :: count,rstcount !ilen,
1424 character(len=50) :: tytul
1425 integer :: maxcount_scale = 10
1426 !el common /gucio/ cm
1427 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1428 !el common /stochcalc/ stochforcvec
1429 integer :: itime,itt,i,j,itsplit
1431 !el common /cipiszcze/ itt
1433 real(kind=8) :: epdrift,tt0,epdriftmax
1436 if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres)) !(MAXRES6) maxres6=6*maxres
1440 if (large.and. mod(itime,ntwe).eq.0) then
1441 write (iout,*) "***************** RESPA itime",itime
1442 write (iout,*) "Cartesian and internal coordinates: step 0"
1444 call pdbout(0.0d0,"cipiszcze",iout)
1446 write (iout,*) "Accelerations from long-range forces"
1448 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1449 (d_a(j,i+nres),j=1,3)
1451 write (iout,*) "Velocities, step 0"
1453 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1454 (d_t(j,i+nres),j=1,3)
1459 ! Perform the initial RESPA step (increment velocities)
1460 ! write (iout,*) "*********************** RESPA ini"
1463 if (mod(itime,ntwe).eq.0 .and. large) then
1464 write (iout,*) "Velocities, end"
1466 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1467 (d_t(j,i+nres),j=1,3)
1471 ! Compute the short-range forces
1477 ! 7/2/2009 commented out
1479 ! call etotal_short(energia_short)
1482 ! 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
1485 d_a(j,i)=d_a_short(j,i)
1489 if (large.and. mod(itime,ntwe).eq.0) then
1490 write (iout,*) "energia_short",energia_short(0)
1491 write (iout,*) "Accelerations from short-range forces"
1493 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1494 (d_a(j,i+nres),j=1,3)
1499 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1501 t_enegrad=t_enegrad+tcpu()-tt0
1506 d_t_old(j,i)=d_t(j,i)
1507 d_a_old(j,i)=d_a(j,i)
1510 ! 6/30/08 A-MTS: attempt at increasing the split number
1513 dc_old0(j,i)=dc_old(j,i)
1514 d_t_old0(j,i)=d_t_old(j,i)
1515 d_a_old0(j,i)=d_a_old(j,i)
1518 if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
1519 if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
1526 ! write (iout,*) "itime",itime," ntime_split",ntime_split
1527 ! Split the time step
1528 d_time=d_time0/ntime_split
1529 ! Perform the short-range RESPA steps (velocity Verlet increments of
1530 ! positions and velocities using short-range forces)
1531 ! write (iout,*) "*********************** RESPA split"
1532 do itsplit=1,ntime_split
1535 else if (lang.eq.2 .or. lang.eq.3) then
1537 call stochastic_force(stochforcvec)
1540 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
1542 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1547 ! First step of the velocity Verlet algorithm
1552 else if (lang.eq.3) then
1554 call sd_verlet1_ciccotti
1556 else if (lang.eq.1) then
1561 ! Build the chain from the newly calculated coordinates
1562 call chainbuild_cart
1563 if (rattle) call rattle1
1565 if (large.and. mod(itime,ntwe).eq.0) then
1566 write (iout,*) "***** ITSPLIT",itsplit
1567 write (iout,*) "Cartesian and internal coordinates: step 1"
1568 call pdbout(0.0d0,"cipiszcze",iout)
1571 write (iout,*) "Velocities, step 1"
1573 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1574 (d_t(j,i+nres),j=1,3)
1583 ! Calculate energy and forces
1585 call etotal_short(energia_short)
1586 ! AL 4/17/17: Exit itime_split loop when energy goes infinite
1587 if (energia_short(0).gt.0.99e20 .or. isnan(energia_short(0)) ) then
1588 if (PRINT_AMTS_MSG) &
1589 write (iout,*) "Infinities/NaNs in energia_short",energia_short(0),"; increasing ntime_split to",ntime_split
1590 ntime_split=ntime_split*2
1591 if (ntime_split.gt.maxtime_split) then
1594 "Cannot rescue the run; aborting job. Retry with a smaller time step"
1596 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1599 "Cannot rescue the run; terminating. Retry with a smaller time step"
1605 if (large.and. mod(itime,ntwe).eq.0) &
1606 call enerprint(energia_short)
1609 t_eshort=t_eshort+MPI_Wtime()-tt0
1611 t_eshort=t_eshort+tcpu()-tt0
1615 ! Get the new accelerations
1617 ! 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
1620 d_a_short(j,i)=d_a(j,i)
1624 if (large.and. mod(itime,ntwe).eq.0) then
1625 write (iout,*)"energia_short",energia_short(0)
1626 write (iout,*) "Accelerations from short-range forces"
1628 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1629 (d_a(j,i+nres),j=1,3)
1634 ! Determine maximum acceleration and scale down the timestep if needed
1636 amax=amax/ntime_split**2
1637 call predict_edrift(epdrift)
1638 if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0) &
1639 write (iout,*) "amax",amax," damax",damax,&
1640 " epdrift",epdrift," epdriftmax",epdriftmax
1641 ! Exit loop and try with increased split number if the change of
1642 ! acceleration is too big
1643 if (amax.gt.damax .or. epdrift.gt.edriftmax) then
1644 if (ntime_split.lt.maxtime_split) then
1646 ntime_split=ntime_split*2
1647 ! AL 4/17/17: We should exit the itime_split loop when acceleration change is too big
1651 dc_old(j,i)=dc_old0(j,i)
1652 d_t_old(j,i)=d_t_old0(j,i)
1653 d_a_old(j,i)=d_a_old0(j,i)
1656 if (PRINT_AMTS_MSG) then
1657 write (iout,*) "acceleration/energy drift too large",amax,&
1658 epdrift," split increased to ",ntime_split," itime",itime,&
1664 "Uh-hu. Bumpy landscape. Maximum splitting number",&
1666 " already reached!!! Trying to carry on!"
1670 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1672 t_enegrad=t_enegrad+tcpu()-tt0
1674 ! Second step of the velocity Verlet algorithm
1679 else if (lang.eq.3) then
1681 call sd_verlet2_ciccotti
1683 else if (lang.eq.1) then
1688 if (rattle) call rattle2
1689 ! Backup the coordinates, velocities, and accelerations
1693 d_t_old(j,i)=d_t(j,i)
1694 d_a_old(j,i)=d_a(j,i)
1701 ! Restore the time step
1703 ! Compute long-range forces
1710 call etotal_long(energia_long)
1711 if (energia_long(0).gt.0.99e20 .or. isnan(energia_long(0))) then
1714 "Infinitied/NaNs in energia_long, Aborting MPI job."
1716 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1718 write (iout,*) "Infinitied/NaNs in energia_long, terminating."
1722 if (large.and. mod(itime,ntwe).eq.0) &
1723 call enerprint(energia_long)
1726 t_elong=t_elong+MPI_Wtime()-tt0
1728 t_elong=t_elong+tcpu()-tt0
1734 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1736 t_enegrad=t_enegrad+tcpu()-tt0
1738 ! Compute accelerations from long-range forces
1740 if (large.and. mod(itime,ntwe).eq.0) then
1741 write (iout,*) "energia_long",energia_long(0)
1742 write (iout,*) "Cartesian and internal coordinates: step 2"
1744 call pdbout(0.0d0,"cipiszcze",iout)
1746 write (iout,*) "Accelerations from long-range forces"
1748 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1749 (d_a(j,i+nres),j=1,3)
1751 write (iout,*) "Velocities, step 2"
1753 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1754 (d_t(j,i+nres),j=1,3)
1758 ! Compute the final RESPA step (increment velocities)
1759 ! write (iout,*) "*********************** RESPA fin"
1761 ! Compute the complete potential energy
1763 potEcomp(i)=energia_short(i)+energia_long(i)
1765 potE=potEcomp(0)-potEcomp(20)
1766 ! potE=energia_short(0)+energia_long(0)
1769 ! Calculate the kinetic and the total energy and the kinetic temperature
1772 ! Couple the system to Berendsen bath if needed
1773 if (tbf .and. lang.eq.0) then
1776 kinetic_T=2.0d0/(dimen3*Rb)*EK
1777 ! Backup the coordinates, velocities, and accelerations
1779 if (mod(itime,ntwe).eq.0 .and. large) then
1780 write (iout,*) "Velocities, end"
1782 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1783 (d_t(j,i+nres),j=1,3)
1788 end subroutine RESPA_step
1789 !-----------------------------------------------------------------------------
1790 subroutine RESPA_vel
1791 ! First and last RESPA step (incrementing velocities using long-range
1794 ! implicit real*8 (a-h,o-z)
1795 ! include 'DIMENSIONS'
1796 ! include 'COMMON.CONTROL'
1797 ! include 'COMMON.VAR'
1798 ! include 'COMMON.MD'
1799 ! include 'COMMON.CHAIN'
1800 ! include 'COMMON.DERIV'
1801 ! include 'COMMON.GEO'
1802 ! include 'COMMON.LOCAL'
1803 ! include 'COMMON.INTERACT'
1804 ! include 'COMMON.IOUNITS'
1805 ! include 'COMMON.NAMES'
1806 integer :: i,j,inres,mnum
1809 d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
1813 d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
1818 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
1819 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
1820 .and.(mnum.ne.5)) then
1823 d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
1828 end subroutine RESPA_vel
1829 !-----------------------------------------------------------------------------
1831 ! Applying velocity Verlet algorithm - step 1 to coordinates
1833 ! implicit real*8 (a-h,o-z)
1834 ! include 'DIMENSIONS'
1835 ! include 'COMMON.CONTROL'
1836 ! include 'COMMON.VAR'
1837 ! include 'COMMON.MD'
1838 ! include 'COMMON.CHAIN'
1839 ! include 'COMMON.DERIV'
1840 ! include 'COMMON.GEO'
1841 ! include 'COMMON.LOCAL'
1842 ! include 'COMMON.INTERACT'
1843 ! include 'COMMON.IOUNITS'
1844 ! include 'COMMON.NAMES'
1845 real(kind=8) :: adt,adt2
1846 integer :: i,j,inres,mnum
1849 write (iout,*) "VELVERLET1 START: DC"
1851 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),&
1852 (dc(j,i+nres),j=1,3)
1856 adt=d_a_old(j,0)*d_time
1858 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1859 d_t_new(j,0)=d_t_old(j,0)+adt2
1860 d_t(j,0)=d_t_old(j,0)+adt
1864 adt=d_a_old(j,i)*d_time
1866 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1867 d_t_new(j,i)=d_t_old(j,i)+adt2
1868 d_t(j,i)=d_t_old(j,i)+adt
1873 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
1874 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
1875 .and.(mnum.ne.5)) then
1878 adt=d_a_old(j,inres)*d_time
1880 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1881 d_t_new(j,inres)=d_t_old(j,inres)+adt2
1882 d_t(j,inres)=d_t_old(j,inres)+adt
1887 write (iout,*) "VELVERLET1 END: DC"
1889 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),&
1890 (dc(j,i+nres),j=1,3)
1894 end subroutine verlet1
1895 !-----------------------------------------------------------------------------
1897 ! Step 2 of the velocity Verlet algorithm: update velocities
1899 ! implicit real*8 (a-h,o-z)
1900 ! include 'DIMENSIONS'
1901 ! include 'COMMON.CONTROL'
1902 ! include 'COMMON.VAR'
1903 ! include 'COMMON.MD'
1904 ! include 'COMMON.CHAIN'
1905 ! include 'COMMON.DERIV'
1906 ! include 'COMMON.GEO'
1907 ! include 'COMMON.LOCAL'
1908 ! include 'COMMON.INTERACT'
1909 ! include 'COMMON.IOUNITS'
1910 ! include 'COMMON.NAMES'
1911 integer :: i,j,inres,mnum
1914 d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
1918 d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
1923 ! iti=iabs(itype(i,mnum))
1924 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
1925 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
1926 .and.(mnum.ne.5)) then
1929 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
1934 end subroutine verlet2
1935 !-----------------------------------------------------------------------------
1936 subroutine sddir_precalc
1937 ! Applying velocity Verlet algorithm - step 1 to coordinates
1938 ! implicit real*8 (a-h,o-z)
1939 ! include 'DIMENSIONS'
1945 ! include 'COMMON.CONTROL'
1946 ! include 'COMMON.VAR'
1947 ! include 'COMMON.MD'
1949 ! include 'COMMON.LANGEVIN'
1951 ! include 'COMMON.LANGEVIN.lang0'
1953 ! include 'COMMON.CHAIN'
1954 ! include 'COMMON.DERIV'
1955 ! include 'COMMON.GEO'
1956 ! include 'COMMON.LOCAL'
1957 ! include 'COMMON.INTERACT'
1958 ! include 'COMMON.IOUNITS'
1959 ! include 'COMMON.NAMES'
1960 ! include 'COMMON.TIME1'
1961 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1962 !el common /stochcalc/ stochforcvec
1963 real(kind=8) :: time00
1965 ! Compute friction and stochastic forces
1970 time_fric=time_fric+MPI_Wtime()-time00
1972 call stochastic_force(stochforcvec)
1973 time_stoch=time_stoch+MPI_Wtime()-time00
1976 ! Compute the acceleration due to friction forces (d_af_work) and stochastic
1977 ! forces (d_as_work)
1979 call ginv_mult(fric_work, d_af_work)
1980 call ginv_mult(stochforcvec, d_as_work)
1982 end subroutine sddir_precalc
1983 !-----------------------------------------------------------------------------
1984 subroutine sddir_verlet1
1985 ! Applying velocity Verlet algorithm - step 1 to velocities
1988 ! implicit real*8 (a-h,o-z)
1989 ! include 'DIMENSIONS'
1990 ! include 'COMMON.CONTROL'
1991 ! include 'COMMON.VAR'
1992 ! include 'COMMON.MD'
1994 ! include 'COMMON.LANGEVIN'
1996 ! include 'COMMON.LANGEVIN.lang0'
1998 ! include 'COMMON.CHAIN'
1999 ! include 'COMMON.DERIV'
2000 ! include 'COMMON.GEO'
2001 ! include 'COMMON.LOCAL'
2002 ! include 'COMMON.INTERACT'
2003 ! include 'COMMON.IOUNITS'
2004 ! include 'COMMON.NAMES'
2005 ! Revised 3/31/05 AL: correlation between random contributions to
2006 ! position and velocity increments included.
2007 real(kind=8) :: sqrt13 = 0.57735026918962576451d0 ! 1/sqrt(3)
2008 real(kind=8) :: adt,adt2
2009 integer :: i,j,ind,inres,mnum
2011 ! Add the contribution from BOTH friction and stochastic force to the
2012 ! coordinates, but ONLY the contribution from the friction forces to velocities
2015 adt=(d_a_old(j,0)+d_af_work(j))*d_time
2016 adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
2017 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
2018 d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
2019 d_t(j,0)=d_t_old(j,0)+adt
2024 adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
2025 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
2026 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
2027 d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
2028 d_t(j,i)=d_t_old(j,i)+adt
2034 ! iti=iabs(itype(i,mnum))
2035 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2036 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2037 .and.(mnum.ne.5)) then
2040 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
2041 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
2042 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
2043 d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
2044 d_t(j,inres)=d_t_old(j,inres)+adt
2050 end subroutine sddir_verlet1
2051 !-----------------------------------------------------------------------------
2052 subroutine sddir_verlet2
2053 ! Calculating the adjusted velocities for accelerations
2056 ! implicit real*8 (a-h,o-z)
2057 ! include 'DIMENSIONS'
2058 ! include 'COMMON.CONTROL'
2059 ! include 'COMMON.VAR'
2060 ! include 'COMMON.MD'
2062 ! include 'COMMON.LANGEVIN'
2064 ! include 'COMMON.LANGEVIN.lang0'
2066 ! include 'COMMON.CHAIN'
2067 ! include 'COMMON.DERIV'
2068 ! include 'COMMON.GEO'
2069 ! include 'COMMON.LOCAL'
2070 ! include 'COMMON.INTERACT'
2071 ! include 'COMMON.IOUNITS'
2072 ! include 'COMMON.NAMES'
2073 real(kind=8),dimension(6*nres) :: stochforcvec,d_as_work1 !(MAXRES6) maxres6=6*maxres
2074 real(kind=8) :: cos60 = 0.5d0, sin60 = 0.86602540378443864676d0
2075 integer :: i,j,ind,inres,mnum
2076 ! Revised 3/31/05 AL: correlation between random contributions to
2077 ! position and velocity increments included.
2078 ! The correlation coefficients are calculated at low-friction limit.
2079 ! Also, friction forces are now not calculated with new velocities.
2081 ! call friction_force
2082 call stochastic_force(stochforcvec)
2084 ! Compute the acceleration due to friction forces (d_af_work) and stochastic
2085 ! forces (d_as_work)
2087 call ginv_mult(stochforcvec, d_as_work1)
2093 d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j)) &
2094 +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
2099 d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j)) &
2100 +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
2106 ! iti=iabs(itype(i,mnum))
2107 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2108 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2109 .and.(mnum.ne.5)) then
2112 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres) &
2113 +d_af_work(ind+j))+sin60*d_as_work(ind+j) &
2114 +cos60*d_as_work1(ind+j))*d_time
2120 end subroutine sddir_verlet2
2121 !-----------------------------------------------------------------------------
2122 subroutine max_accel
2124 ! Find the maximum difference in the accelerations of the the sites
2125 ! at the beginning and the end of the time step.
2128 ! implicit real*8 (a-h,o-z)
2129 ! include 'DIMENSIONS'
2130 ! include 'COMMON.CONTROL'
2131 ! include 'COMMON.VAR'
2132 ! include 'COMMON.MD'
2133 ! include 'COMMON.CHAIN'
2134 ! include 'COMMON.DERIV'
2135 ! include 'COMMON.GEO'
2136 ! include 'COMMON.LOCAL'
2137 ! include 'COMMON.INTERACT'
2138 ! include 'COMMON.IOUNITS'
2139 real(kind=8),dimension(3) :: aux,accel,accel_old
2140 real(kind=8) :: dacc
2144 ! aux(j)=d_a(j,0)-d_a_old(j,0)
2145 accel_old(j)=d_a_old(j,0)
2152 ! 7/3/08 changed to asymmetric difference
2154 ! accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
2155 accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
2156 accel(j)=accel(j)+0.5d0*d_a(j,i)
2157 ! if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
2158 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
2159 dacc=dabs(accel(j)-accel_old(j))
2160 ! write (iout,*) i,dacc
2161 if (dacc.gt.amax) amax=dacc
2169 accel_old(j)=d_a_old(j,0)
2174 accel_old(j)=accel_old(j)+d_a_old(j,1)
2175 accel(j)=accel(j)+d_a(j,1)
2180 ! iti=iabs(itype(i,mnum))
2181 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2182 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2183 .and.(mnum.ne.5)) then
2185 ! accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
2186 accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
2187 accel(j)=accel(j)+d_a(j,i+nres)
2191 ! if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
2192 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
2193 dacc=dabs(accel(j)-accel_old(j))
2194 ! write (iout,*) "side-chain",i,dacc
2195 if (dacc.gt.amax) amax=dacc
2199 accel_old(j)=accel_old(j)+d_a_old(j,i)
2200 accel(j)=accel(j)+d_a(j,i)
2201 ! aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
2205 end subroutine max_accel
2206 !-----------------------------------------------------------------------------
2207 subroutine predict_edrift(epdrift)
2209 ! Predict the drift of the potential energy
2212 use control_data, only: lmuca
2213 ! implicit real*8 (a-h,o-z)
2214 ! include 'DIMENSIONS'
2215 ! include 'COMMON.CONTROL'
2216 ! include 'COMMON.VAR'
2217 ! include 'COMMON.MD'
2218 ! include 'COMMON.CHAIN'
2219 ! include 'COMMON.DERIV'
2220 ! include 'COMMON.GEO'
2221 ! include 'COMMON.LOCAL'
2222 ! include 'COMMON.INTERACT'
2223 ! include 'COMMON.IOUNITS'
2224 ! include 'COMMON.MUCA'
2225 real(kind=8) :: epdrift,epdriftij
2227 ! Drift of the potential energy
2233 epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
2234 if (lmuca) epdriftij=epdriftij*factor
2235 ! write (iout,*) "back",i,j,epdriftij
2236 if (epdriftij.gt.epdrift) epdrift=epdriftij
2240 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1.and.&
2241 molnum(i).ne.5) then
2244 dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
2245 if (lmuca) epdriftij=epdriftij*factor
2246 ! write (iout,*) "side",i,j,epdriftij
2247 if (epdriftij.gt.epdrift) epdrift=epdriftij
2251 epdrift=0.5d0*epdrift*d_time*d_time
2252 ! write (iout,*) "epdrift",epdrift
2254 end subroutine predict_edrift
2255 !-----------------------------------------------------------------------------
2256 subroutine verlet_bath
2258 ! Coupling to the thermostat by using the Berendsen algorithm
2261 ! implicit real*8 (a-h,o-z)
2262 ! include 'DIMENSIONS'
2263 ! include 'COMMON.CONTROL'
2264 ! include 'COMMON.VAR'
2265 ! include 'COMMON.MD'
2266 ! include 'COMMON.CHAIN'
2267 ! include 'COMMON.DERIV'
2268 ! include 'COMMON.GEO'
2269 ! include 'COMMON.LOCAL'
2270 ! include 'COMMON.INTERACT'
2271 ! include 'COMMON.IOUNITS'
2272 ! include 'COMMON.NAMES'
2273 real(kind=8) :: T_half,fact
2274 integer :: i,j,inres,mnum
2276 T_half=2.0d0/(dimen3*Rb)*EK
2277 fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
2278 ! write(iout,*) "T_half", T_half
2279 ! write(iout,*) "EK", EK
2280 ! write(iout,*) "fact", fact
2282 d_t(j,0)=fact*d_t(j,0)
2286 d_t(j,i)=fact*d_t(j,i)
2291 ! iti=iabs(itype(i,mnum))
2292 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2293 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2294 .and.(mnum.ne.5)) then
2297 d_t(j,inres)=fact*d_t(j,inres)
2302 end subroutine verlet_bath
2303 !-----------------------------------------------------------------------------
2305 ! Set up the initial conditions of a MD simulation
2308 use control, only:tcpu
2309 !el use io_basic, only:ilen
2312 use minimm, only:minim_dc,minimize,sc_move
2313 use io_config, only:readrst
2314 use io, only:statout
2315 ! implicit real*8 (a-h,o-z)
2316 ! include 'DIMENSIONS'
2319 character(len=16) :: form
2320 integer :: IERROR,ERRCODE
2322 ! include 'COMMON.SETUP'
2323 ! include 'COMMON.CONTROL'
2324 ! include 'COMMON.VAR'
2325 ! include 'COMMON.MD'
2327 ! include 'COMMON.LANGEVIN'
2329 ! include 'COMMON.LANGEVIN.lang0'
2331 ! include 'COMMON.CHAIN'
2332 ! include 'COMMON.DERIV'
2333 ! include 'COMMON.GEO'
2334 ! include 'COMMON.LOCAL'
2335 ! include 'COMMON.INTERACT'
2336 ! include 'COMMON.IOUNITS'
2337 ! include 'COMMON.NAMES'
2338 ! include 'COMMON.REMD'
2339 real(kind=8),dimension(0:n_ene) :: energia_long,energia_short
2340 real(kind=8),dimension(3) :: vcm,incr,L
2341 real(kind=8) :: xv,sigv,lowb,highb
2342 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
2343 character(len=256) :: qstr
2346 character(len=50) :: tytul
2347 logical :: file_exist
2348 !el common /gucio/ cm
2349 integer :: i,j,ipos,iq,iw,nft_sc,iretcode,nfun,itime,ierr,mnum
2350 real(kind=8) :: etot,tt0
2354 ! write(iout,*) "d_time", d_time
2355 ! Compute the standard deviations of stochastic forces for Langevin dynamics
2356 ! if the friction coefficients do not depend on surface area
2357 if (lang.gt.0 .and. .not.surfarea) then
2360 stdforcp(i)=stdfp(mnum)*dsqrt(gamp(mnum))
2364 stdforcsc(i)=stdfsc(iabs(itype(i,mnum)),mnum) &
2365 *dsqrt(gamsc(iabs(itype(i,mnum)),mnum))
2368 ! Open the pdb file for snapshotshots
2371 if (ilen(tmpdir).gt.0) &
2372 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2373 liczba(:ilen(liczba))//".pdb")
2375 file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2379 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) &
2380 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2381 liczba(:ilen(liczba))//".x")
2382 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2385 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) &
2386 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2387 liczba(:ilen(liczba))//".cx")
2388 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2394 if (ilen(tmpdir).gt.0) &
2395 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
2396 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
2398 if (ilen(tmpdir).gt.0) &
2399 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
2400 cartname=prefix(:ilen(prefix))//"_MD.cx"
2404 write (qstr,'(256(1h ))')
2407 iq = qinfrag(i,iset)*10
2408 iw = wfrag(i,iset)/100
2410 if(me.eq.king.or..not.out1file) &
2411 write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
2412 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
2417 iq = qinpair(i,iset)*10
2418 iw = wpair(i,iset)/100
2420 if(me.eq.king.or..not.out1file) &
2421 write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
2422 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
2426 ! pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
2428 ! cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
2430 ! cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
2432 ! statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
2436 if (restart1file) then
2438 inquire(file=mremd_rst_name,exist=file_exist)
2439 write (*,*) me," Before broadcast: file_exist",file_exist
2441 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,&
2444 write (*,*) me," After broadcast: file_exist",file_exist
2445 ! inquire(file=mremd_rst_name,exist=file_exist)
2446 if(me.eq.king.or..not.out1file) &
2447 write(iout,*) "Initial state read by master and distributed"
2449 if (ilen(tmpdir).gt.0) &
2450 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_' &
2451 //liczba(:ilen(liczba))//'.rst')
2452 inquire(file=rest2name,exist=file_exist)
2455 if(.not.restart1file) then
2456 if(me.eq.king.or..not.out1file) &
2457 write(iout,*) "Initial state will be read from file ",&
2458 rest2name(:ilen(rest2name))
2461 call rescale_weights(t_bath)
2463 if(me.eq.king.or..not.out1file)then
2464 if (restart1file) then
2465 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),&
2468 write(iout,*) "File ",rest2name(:ilen(rest2name)),&
2471 write(iout,*) "Initial velocities randomly generated"
2478 ! Generate initial velocities
2479 if(me.eq.king.or..not.out1file) &
2480 write(iout,*) "Initial velocities randomly generated"
2485 ! rest2name = prefix(:ilen(prefix))//'.rst'
2486 if(me.eq.king.or..not.out1file)then
2487 write (iout,*) "Initial velocities"
2489 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2490 (d_t(j,i+nres),j=1,3)
2492 ! Zeroing the total angular momentum of the system
2493 write(iout,*) "Calling the zero-angular momentum subroutine"
2496 ! Getting the potential energy and forces and velocities and accelerations
2498 ! write (iout,*) "velocity of the center of the mass:"
2499 ! write (iout,*) (vcm(j),j=1,3)
2501 d_t(j,0)=d_t(j,0)-vcm(j)
2503 ! Removing the velocity of the center of mass
2505 if(me.eq.king.or..not.out1file)then
2506 write (iout,*) "vcm right after adjustment:"
2507 write (iout,*) (vcm(j),j=1,3)
2509 if ((.not.rest).and.(indpdb.eq.0)) then
2511 if(iranconf.ne.0) then
2513 print *, 'Calling OVERLAP_SC'
2514 call overlap_sc(fail)
2517 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
2518 print *,'SC_move',nft_sc,etot
2519 if(me.eq.king.or..not.out1file) &
2520 write(iout,*) 'SC_move',nft_sc,etot
2524 print *, 'Calling MINIM_DC'
2525 call minim_dc(etot,iretcode,nfun)
2527 call geom_to_var(nvar,varia)
2528 print *,'Calling MINIMIZE.'
2529 call minimize(etot,varia,iretcode,nfun)
2530 call var_to_geom(nvar,varia)
2532 if(me.eq.king.or..not.out1file) &
2533 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2536 call chainbuild_cart
2541 kinetic_T=2.0d0/(dimen3*Rb)*EK
2542 if(me.eq.king.or..not.out1file)then
2552 call etotal(potEcomp)
2553 if (large) call enerprint(potEcomp)
2556 t_etotal=t_etotal+MPI_Wtime()-tt0
2558 t_etotal=t_etotal+tcpu()-tt0
2565 if (amax*d_time .gt. dvmax) then
2566 d_time=d_time*dvmax/amax
2567 if(me.eq.king.or..not.out1file) write (iout,*) &
2568 "Time step reduced to",d_time,&
2569 " because of too large initial acceleration."
2571 if(me.eq.king.or..not.out1file)then
2572 write(iout,*) "Potential energy and its components"
2573 call enerprint(potEcomp)
2574 ! write(iout,*) (potEcomp(i),i=0,n_ene)
2576 potE=potEcomp(0)-potEcomp(20)
2579 if (ntwe.ne.0) call statout(itime)
2580 if(me.eq.king.or..not.out1file) &
2581 write (iout,'(/a/3(a25,1pe14.5/))') "Initial:", &
2582 " Kinetic energy",EK," Potential energy",potE, &
2583 " Total energy",totE," Maximum acceleration ", &
2586 write (iout,*) "Initial coordinates"
2588 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),&
2591 write (iout,*) "Initial dC"
2593 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),&
2594 (dc(j,i+nres),j=1,3)
2596 write (iout,*) "Initial velocities"
2597 write (iout,"(13x,' backbone ',23x,' side chain')")
2599 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2600 (d_t(j,i+nres),j=1,3)
2602 write (iout,*) "Initial accelerations"
2604 ! write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2605 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),&
2606 (d_a(j,i+nres),j=1,3)
2612 d_t_old(j,i)=d_t(j,i)
2613 d_a_old(j,i)=d_a(j,i)
2615 ! write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
2624 call etotal_short(energia_short)
2625 if (large) call enerprint(potEcomp)
2628 t_eshort=t_eshort+MPI_Wtime()-tt0
2630 t_eshort=t_eshort+tcpu()-tt0
2635 if(.not.out1file .and. large) then
2636 write (iout,*) "energia_long",energia_long(0),&
2637 " energia_short",energia_short(0),&
2638 " total",energia_long(0)+energia_short(0)
2639 write (iout,*) "Initial fast-force accelerations"
2641 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
2642 (d_a(j,i+nres),j=1,3)
2645 ! 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
2648 d_a_short(j,i)=d_a(j,i)
2657 call etotal_long(energia_long)
2658 if (large) call enerprint(potEcomp)
2661 t_elong=t_elong+MPI_Wtime()-tt0
2663 t_elong=t_elong+tcpu()-tt0
2668 if(.not.out1file .and. large) then
2669 write (iout,*) "energia_long",energia_long(0)
2670 write (iout,*) "Initial slow-force accelerations"
2672 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
2673 (d_a(j,i+nres),j=1,3)
2677 t_enegrad=t_enegrad+MPI_Wtime()-tt0
2679 t_enegrad=t_enegrad+tcpu()-tt0
2683 end subroutine init_MD
2684 !-----------------------------------------------------------------------------
2685 subroutine random_vel
2687 ! implicit real*8 (a-h,o-z)
2689 use random, only:anorm_distr
2691 ! include 'DIMENSIONS'
2692 ! include 'COMMON.CONTROL'
2693 ! include 'COMMON.VAR'
2694 ! include 'COMMON.MD'
2696 ! include 'COMMON.LANGEVIN'
2698 ! include 'COMMON.LANGEVIN.lang0'
2700 ! include 'COMMON.CHAIN'
2701 ! include 'COMMON.DERIV'
2702 ! include 'COMMON.GEO'
2703 ! include 'COMMON.LOCAL'
2704 ! include 'COMMON.INTERACT'
2705 ! include 'COMMON.IOUNITS'
2706 ! include 'COMMON.NAMES'
2707 ! include 'COMMON.TIME1'
2708 real(kind=8) :: xv,sigv,lowb,highb ,Ek1
2710 real(kind=8) ,allocatable, dimension(:) :: DDU1,DDU2,DL2,DL1,xsolv,DML,rs
2711 real(kind=8) :: sumx
2713 real(kind=8) ,allocatable, dimension(:) :: rsold
2714 real (kind=8),allocatable,dimension(:,:) :: matold
2717 integer :: i,j,ii,k,ind,mark,imark,mnum
2718 ! Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
2719 ! First generate velocities in the eigenspace of the G matrix
2720 ! write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
2727 sigv=dsqrt((Rb*t_bath)/geigen(i))
2730 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
2732 write (iout,*) "i",i," ii",ii," geigen",geigen(i),&
2733 " d_t_work_new",d_t_work_new(ii)
2744 Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
2747 write (iout,*) "Ek from eigenvectors",Ek1
2748 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
2752 ! Transform velocities to UNRES coordinate space
2753 allocate (DL1(2*nres))
2754 allocate (DDU1(2*nres))
2755 allocate (DL2(2*nres))
2756 allocate (DDU2(2*nres))
2757 allocate (xsolv(2*nres))
2758 allocate (DML(2*nres))
2759 allocate (rs(2*nres))
2761 allocate (rsold(2*nres))
2762 allocate (matold(2*nres,2*nres))
2764 matold(1,1)=DMorig(1)
2765 matold(1,2)=DU1orig(1)
2766 matold(1,3)=DU2orig(1)
2767 write (*,*) DMorig(1),DU1orig(1),DU2orig(1)
2772 matold(i,j)=DMorig(i)
2773 matold(i,j-1)=DU1orig(i-1)
2775 matold(i,j-2)=DU2orig(i-2)
2783 matold(i,j+1)=DU1orig(i)
2789 matold(i,j+2)=DU2orig(i)
2793 matold(dimen,dimen)=DMorig(dimen)
2794 matold(dimen,dimen-1)=DU1orig(dimen-1)
2795 matold(dimen,dimen-2)=DU2orig(dimen-2)
2796 write(iout,*) "old gmatrix"
2797 call matout(dimen,dimen,2*nres,2*nres,matold)
2801 ! Find the ith eigenvector of the pentadiagonal inertiq matrix
2805 DML(j)=DMorig(j)-geigen(i)
2808 DML(j-1)=DMorig(j)-geigen(i)
2813 DDU1(imark-1)=DU2orig(imark-1)
2814 do j=imark+1,dimen-1
2815 DDU1(j-1)=DU1orig(j)
2823 DDU2(j)=DU2orig(j+1)
2832 write (iout,*) "DL2,DL1,DML,DDU1,DDU2"
2833 write(iout,'(10f10.5)') (DL2(k),k=3,dimen-1)
2834 write(iout,'(10f10.5)') (DL1(k),k=2,dimen-1)
2835 write(iout,'(10f10.5)') (DML(k),k=1,dimen-1)
2836 write(iout,'(10f10.5)') (DDU1(k),k=1,dimen-2)
2837 write(iout,'(10f10.5)') (DDU2(k),k=1,dimen-3)
2840 if (imark.gt.2) rs(imark-2)=-DU2orig(imark-2)
2841 if (imark.gt.1) rs(imark-1)=-DU1orig(imark-1)
2842 if (imark.lt.dimen) rs(imark)=-DU1orig(imark)
2843 if (imark.lt.dimen-1) rs(imark+1)=-DU2orig(imark)
2847 ! write (iout,*) "Vector rs"
2849 ! write (iout,*) j,rs(j)
2852 call FDIAG(dimen-1,DL2,DL1,DML,DDU1,DDU2,rs,xsolv,mark)
2859 sumx=-geigen(i)*xsolv(j)
2861 sumx=sumx+matold(j,k)*xsolv(k)
2864 sumx=sumx+matold(j,k)*xsolv(k-1)
2866 write(iout,'(i5,3f10.5)') j,sumx,rsold(j),sumx-rsold(j)
2869 sumx=-geigen(i)*xsolv(j-1)
2871 sumx=sumx+matold(j,k)*xsolv(k)
2874 sumx=sumx+matold(j,k)*xsolv(k-1)
2876 write(iout,'(i5,3f10.5)') j-1,sumx,rsold(j-1),sumx-rsold(j-1)
2880 "Solution of equations system",i," for eigenvalue",geigen(i)
2882 write(iout,'(i5,f10.5)') j,xsolv(j)
2885 do j=dimen-1,imark,-1
2890 write (iout,*) "Un-normalized eigenvector",i," for eigenvalue",geigen(i)
2892 write(iout,'(i5,f10.5)') j,xsolv(j)
2895 ! Normalize ith eigenvector
2898 sumx=sumx+xsolv(j)**2
2902 xsolv(j)=xsolv(j)/sumx
2905 write (iout,*) "Eigenvector",i," for eigenvalue",geigen(i)
2907 write(iout,'(i5,3f10.5)') j,xsolv(j)
2910 ! All done at this point for eigenvector i, exit loop
2918 write (iout,*) "Unable to find eigenvector",i
2921 ! write (iout,*) "i=",i
2923 ! write (iout,*) "k=",k
2926 ! write(iout,*) "ind",ind," ind1",3*(i-1)+k
2927 d_t_work(ind)=d_t_work(ind) &
2928 +xsolv(j)*d_t_work_new(3*(i-1)+k)
2931 enddo ! i (loop over the eigenvectors)
2934 write (iout,*) "d_t_work"
2936 write (iout,"(i5,f10.5)") i,d_t_work(i)
2941 ! if (itype(i,1).eq.10) then
2943 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
2944 iti=iabs(itype(i,mnum))
2945 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2946 if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
2947 .or.(mnum.eq.5)) then
2954 ! write (iout,*) "k",k," ii+k",ii+k," ii+j+k",ii+j+k,"EK1",Ek1
2955 Ek1=Ek1+0.5d0*mp(mnum)*((d_t_work(ii+j+k)+d_t_work_new(ii+k))/2)**2+&
2956 0.5d0*0.25d0*IP(mnum)*(d_t_work(ii+j+k)-d_t_work_new(ii+k))**2
2959 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2960 .and.(mnum.ne.5)) ii=ii+3
2961 write (iout,*) "i",i," itype",itype(i,mnum)," mass",msc(itype(i,mnum),mnum)
2962 write (iout,*) "ii",ii
2965 write (iout,*) "k",k," ii",ii,"EK1",EK1
2966 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2968 Ek1=Ek1+0.5d0*Isc(iabs(itype(i,mnum),mnum))*(d_t_work(ii)-d_t_work(ii-3))**2
2969 Ek1=Ek1+0.5d0*msc(iabs(itype(i,mnum)),mnum)*d_t_work(ii)**2
2971 write (iout,*) "i",i," ii",ii
2973 write (iout,*) "Ek from d_t_work",Ek1
2974 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
2990 d_t(k,j)=d_t_work(ind)
2994 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2995 .and.(mnum.ne.5)) then
2997 d_t(k,j+nres)=d_t_work(ind)
3003 write (iout,*) "Random velocities in the Calpha,SC space"
3005 write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3)
3008 write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3)
3015 ! if (itype(i,1).eq.10) then
3017 if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
3018 .or.(mnum.eq.5)) then
3020 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
3024 d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
3025 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
3030 write (iout,*)"Random velocities in the virtual-bond-vector space"
3032 write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3)
3035 write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3)
3038 write (iout,*) "Ek from d_t_work",Ek1
3039 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
3047 d_t_work(ind)=d_t_work(ind) &
3048 +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
3050 ! write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
3054 ! Transfer to the d_t vector
3056 d_t(j,0)=d_t_work(j)
3062 d_t(j,i)=d_t_work(ind)
3067 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3068 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3069 .and.(mnum.ne.5)) then
3072 d_t(j,i+nres)=d_t_work(ind)
3078 ! write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",&
3079 ! 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
3083 end subroutine random_vel
3084 !-----------------------------------------------------------------------------
3086 subroutine sd_verlet_p_setup
3087 ! Sets up the parameters of stochastic Verlet algorithm
3088 ! implicit real*8 (a-h,o-z)
3089 ! include 'DIMENSIONS'
3090 use control, only: tcpu
3095 ! include 'COMMON.CONTROL'
3096 ! include 'COMMON.VAR'
3097 ! include 'COMMON.MD'
3099 ! include 'COMMON.LANGEVIN'
3101 ! include 'COMMON.LANGEVIN.lang0'
3103 ! include 'COMMON.CHAIN'
3104 ! include 'COMMON.DERIV'
3105 ! include 'COMMON.GEO'
3106 ! include 'COMMON.LOCAL'
3107 ! include 'COMMON.INTERACT'
3108 ! include 'COMMON.IOUNITS'
3109 ! include 'COMMON.NAMES'
3110 ! include 'COMMON.TIME1'
3111 real(kind=8),dimension(6*nres) :: emgdt !(MAXRES6) maxres6=6*maxres
3112 real(kind=8) :: pterm,vterm,rho,rhoc,vsig
3113 real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
3114 prand_vec,vrand_vec1,vrand_vec2 !(MAXRES6) maxres6=6*maxres
3115 logical :: lprn = .false.
3116 real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
3117 real(kind=8) :: ktm,gdt,egdt,gdt2,gdt3,gdt4,gdt5,gdt6,gdt7,gdt8,&
3119 integer :: i,maxres2
3126 ! AL 8/17/04 Code adapted from tinker
3128 ! Get the frictional and random terms for stochastic dynamics in the
3129 ! eigenspace of mass-scaled UNRES friction matrix
3133 gdt = fricgam(i) * d_time
3135 ! Stochastic dynamics reduces to simple MD for zero friction
3137 if (gdt .le. zero) then
3138 pfric_vec(i) = 1.0d0
3139 vfric_vec(i) = d_time
3140 afric_vec(i) = 0.5d0 * d_time * d_time
3141 prand_vec(i) = 0.0d0
3142 vrand_vec1(i) = 0.0d0
3143 vrand_vec2(i) = 0.0d0
3145 ! Analytical expressions when friction coefficient is large
3148 if (gdt .ge. gdt_radius) then
3151 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
3152 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
3153 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
3154 vterm = 1.0d0 - egdt**2
3155 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
3157 ! Use series expansions when friction coefficient is small
3168 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0 &
3169 - gdt5/120.0d0 + gdt6/720.0d0 &
3170 - gdt7/5040.0d0 + gdt8/40320.0d0 &
3171 - gdt9/362880.0d0) / fricgam(i)**2
3172 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
3173 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
3174 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0 &
3175 + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0 &
3176 + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0 &
3177 + 127.0d0*gdt9/90720.0d0
3178 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0 &
3179 - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0 &
3180 - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0 &
3181 - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
3182 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0 &
3183 - 17.0d0*gdt2/1280.0d0 &
3184 + 17.0d0*gdt3/6144.0d0 &
3185 + 40967.0d0*gdt4/34406400.0d0 &
3186 - 57203.0d0*gdt5/275251200.0d0 &
3187 - 1429487.0d0*gdt6/13212057600.0d0)
3190 ! Compute the scaling factors of random terms for the nonzero friction case
3192 ktm = 0.5d0*d_time/fricgam(i)
3193 psig = dsqrt(ktm*pterm) / fricgam(i)
3194 vsig = dsqrt(ktm*vterm)
3195 rhoc = dsqrt(1.0d0 - rho*rho)
3197 vrand_vec1(i) = vsig * rho
3198 vrand_vec2(i) = vsig * rhoc
3203 "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",&
3206 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),&
3207 afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
3211 ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
3214 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
3215 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
3216 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
3217 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
3218 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
3219 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
3222 t_sdsetup=t_sdsetup+MPI_Wtime()
3224 t_sdsetup=t_sdsetup+tcpu()-tt0
3227 end subroutine sd_verlet_p_setup
3228 !-----------------------------------------------------------------------------
3229 subroutine eigtransf1(n,ndim,ab,d,c)
3233 real(kind=8) :: ab(ndim,ndim,n),c(ndim,n),d(ndim)
3239 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
3244 end subroutine eigtransf1
3245 !-----------------------------------------------------------------------------
3246 subroutine eigtransf(n,ndim,a,b,d,c)
3250 real(kind=8) :: a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
3256 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
3261 end subroutine eigtransf
3262 !-----------------------------------------------------------------------------
3263 subroutine sd_verlet1
3265 ! Applying stochastic velocity Verlet algorithm - step 1 to velocities
3267 ! implicit real*8 (a-h,o-z)
3268 ! include 'DIMENSIONS'
3269 ! include 'COMMON.CONTROL'
3270 ! include 'COMMON.VAR'
3271 ! include 'COMMON.MD'
3273 ! include 'COMMON.LANGEVIN'
3275 ! include 'COMMON.LANGEVIN.lang0'
3277 ! include 'COMMON.CHAIN'
3278 ! include 'COMMON.DERIV'
3279 ! include 'COMMON.GEO'
3280 ! include 'COMMON.LOCAL'
3281 ! include 'COMMON.INTERACT'
3282 ! include 'COMMON.IOUNITS'
3283 ! include 'COMMON.NAMES'
3284 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
3285 !el common /stochcalc/ stochforcvec
3286 logical :: lprn = .false.
3287 real(kind=8) :: ddt1,ddt2
3288 integer :: i,j,ind,inres
3290 ! write (iout,*) "dc_old"
3292 ! write (iout,'(i5,3f10.5,5x,3f10.5)')
3293 ! & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3296 dc_work(j)=dc_old(j,0)
3297 d_t_work(j)=d_t_old(j,0)
3298 d_a_work(j)=d_a_old(j,0)
3303 dc_work(ind+j)=dc_old(j,i)
3304 d_t_work(ind+j)=d_t_old(j,i)
3305 d_a_work(ind+j)=d_a_old(j,i)
3311 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3312 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3313 .and.(mnum.ne.5)) then
3315 dc_work(ind+j)=dc_old(j,i+nres)
3316 d_t_work(ind+j)=d_t_old(j,i+nres)
3317 d_a_work(ind+j)=d_a_old(j,i+nres)
3325 "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",&
3329 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
3330 vfric_mat(i,j),afric_mat(i,j),&
3331 prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3339 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) &
3340 +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3341 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3342 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3344 d_t_work_new(i)=ddt1+0.5d0*ddt2
3345 d_t_work(i)=ddt1+ddt2
3350 d_t(j,0)=d_t_work(j)
3355 dc(j,i)=dc_work(ind+j)
3356 d_t(j,i)=d_t_work(ind+j)
3362 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3363 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3364 .and.(mnum.ne.5)) then
3367 dc(j,inres)=dc_work(ind+j)
3368 d_t(j,inres)=d_t_work(ind+j)
3374 end subroutine sd_verlet1
3375 !-----------------------------------------------------------------------------
3376 subroutine sd_verlet2
3378 ! Calculating the adjusted velocities for accelerations
3380 ! implicit real*8 (a-h,o-z)
3381 ! include 'DIMENSIONS'
3382 ! include 'COMMON.CONTROL'
3383 ! include 'COMMON.VAR'
3384 ! include 'COMMON.MD'
3386 ! include 'COMMON.LANGEVIN'
3388 ! include 'COMMON.LANGEVIN.lang0'
3390 ! include 'COMMON.CHAIN'
3391 ! include 'COMMON.DERIV'
3392 ! include 'COMMON.GEO'
3393 ! include 'COMMON.LOCAL'
3394 ! include 'COMMON.INTERACT'
3395 ! include 'COMMON.IOUNITS'
3396 ! include 'COMMON.NAMES'
3397 !el real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV !(MAXRES6) maxres6=6*maxres
3398 real(kind=8),dimension(6*nres) :: stochforcvecV !(MAXRES6) maxres6=6*maxres
3399 !el common /stochcalc/ stochforcvec
3401 real(kind=8) :: ddt1,ddt2
3402 integer :: i,j,ind,inres
3403 ! Compute the stochastic forces which contribute to velocity change
3405 call stochastic_force(stochforcvecV)
3412 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3413 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+ &
3414 vrand_mat2(i,j)*stochforcvecV(j)
3416 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
3420 d_t(j,0)=d_t_work(j)
3425 d_t(j,i)=d_t_work(ind+j)
3431 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3432 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3433 .and.(mnum.ne.5)) then
3436 d_t(j,inres)=d_t_work(ind+j)
3442 end subroutine sd_verlet2
3443 !-----------------------------------------------------------------------------
3444 subroutine sd_verlet_ciccotti_setup
3446 ! Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
3448 ! implicit real*8 (a-h,o-z)
3449 ! include 'DIMENSIONS'
3450 use control, only: tcpu
3455 ! include 'COMMON.CONTROL'
3456 ! include 'COMMON.VAR'
3457 ! include 'COMMON.MD'
3459 ! include 'COMMON.LANGEVIN'
3461 ! include 'COMMON.LANGEVIN.lang0'
3463 ! include 'COMMON.CHAIN'
3464 ! include 'COMMON.DERIV'
3465 ! include 'COMMON.GEO'
3466 ! include 'COMMON.LOCAL'
3467 ! include 'COMMON.INTERACT'
3468 ! include 'COMMON.IOUNITS'
3469 ! include 'COMMON.NAMES'
3470 ! include 'COMMON.TIME1'
3471 real(kind=8),dimension(6*nres) :: emgdt !(MAXRES6) maxres6=6*maxres
3472 real(kind=8) :: pterm,vterm,rho,rhoc,vsig
3473 real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
3474 prand_vec,vrand_vec1,vrand_vec2 !(MAXRES6) maxres6=6*maxres
3475 logical :: lprn = .false.
3476 real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
3477 real(kind=8) :: ktm,gdt,egdt,tt0
3478 integer :: i,maxres2
3485 ! AL 8/17/04 Code adapted from tinker
3487 ! Get the frictional and random terms for stochastic dynamics in the
3488 ! eigenspace of mass-scaled UNRES friction matrix
3492 write (iout,*) "i",i," fricgam",fricgam(i)
3493 gdt = fricgam(i) * d_time
3495 ! Stochastic dynamics reduces to simple MD for zero friction
3497 if (gdt .le. zero) then
3498 pfric_vec(i) = 1.0d0
3499 vfric_vec(i) = d_time
3500 afric_vec(i) = 0.5d0*d_time*d_time
3501 prand_vec(i) = afric_vec(i)
3502 vrand_vec2(i) = vfric_vec(i)
3504 ! Analytical expressions when friction coefficient is large
3509 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
3510 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
3511 prand_vec(i) = afric_vec(i)
3512 vrand_vec2(i) = vfric_vec(i)
3514 ! Compute the scaling factors of random terms for the nonzero friction case
3516 ! ktm = 0.5d0*d_time/fricgam(i)
3517 ! psig = dsqrt(ktm*pterm) / fricgam(i)
3518 ! vsig = dsqrt(ktm*vterm)
3519 ! prand_vec(i) = psig*afric_vec(i)
3520 ! vrand_vec2(i) = vsig*vfric_vec(i)
3525 "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",&
3528 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),&
3529 afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
3533 ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
3535 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
3536 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
3537 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
3538 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
3539 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
3541 t_sdsetup=t_sdsetup+MPI_Wtime()
3543 t_sdsetup=t_sdsetup+tcpu()-tt0
3546 end subroutine sd_verlet_ciccotti_setup
3547 !-----------------------------------------------------------------------------
3548 subroutine sd_verlet1_ciccotti
3550 ! Applying stochastic velocity Verlet algorithm - step 1 to velocities
3551 ! implicit real*8 (a-h,o-z)
3553 ! include 'DIMENSIONS'
3557 ! include 'COMMON.CONTROL'
3558 ! include 'COMMON.VAR'
3559 ! include 'COMMON.MD'
3561 ! include 'COMMON.LANGEVIN'
3563 ! include 'COMMON.LANGEVIN.lang0'
3565 ! include 'COMMON.CHAIN'
3566 ! include 'COMMON.DERIV'
3567 ! include 'COMMON.GEO'
3568 ! include 'COMMON.LOCAL'
3569 ! include 'COMMON.INTERACT'
3570 ! include 'COMMON.IOUNITS'
3571 ! include 'COMMON.NAMES'
3572 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
3573 !el common /stochcalc/ stochforcvec
3574 logical :: lprn = .false.
3575 real(kind=8) :: ddt1,ddt2
3576 integer :: i,j,ind,inres
3577 ! write (iout,*) "dc_old"
3579 ! write (iout,'(i5,3f10.5,5x,3f10.5)')
3580 ! & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3583 dc_work(j)=dc_old(j,0)
3584 d_t_work(j)=d_t_old(j,0)
3585 d_a_work(j)=d_a_old(j,0)
3590 dc_work(ind+j)=dc_old(j,i)
3591 d_t_work(ind+j)=d_t_old(j,i)
3592 d_a_work(ind+j)=d_a_old(j,i)
3597 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3599 dc_work(ind+j)=dc_old(j,i+nres)
3600 d_t_work(ind+j)=d_t_old(j,i+nres)
3601 d_a_work(ind+j)=d_a_old(j,i+nres)
3610 "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",&
3614 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
3615 vfric_mat(i,j),afric_mat(i,j),&
3616 prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3624 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) &
3625 +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3626 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3627 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3629 d_t_work_new(i)=ddt1+0.5d0*ddt2
3630 d_t_work(i)=ddt1+ddt2
3635 d_t(j,0)=d_t_work(j)
3640 dc(j,i)=dc_work(ind+j)
3641 d_t(j,i)=d_t_work(ind+j)
3647 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3648 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3649 .and.(mnum.ne.5)) then
3652 dc(j,inres)=dc_work(ind+j)
3653 d_t(j,inres)=d_t_work(ind+j)
3659 end subroutine sd_verlet1_ciccotti
3660 !-----------------------------------------------------------------------------
3661 subroutine sd_verlet2_ciccotti
3663 ! Calculating the adjusted velocities for accelerations
3665 ! implicit real*8 (a-h,o-z)
3666 ! include 'DIMENSIONS'
3667 ! include 'COMMON.CONTROL'
3668 ! include 'COMMON.VAR'
3669 ! include 'COMMON.MD'
3671 ! include 'COMMON.LANGEVIN'
3673 ! include 'COMMON.LANGEVIN.lang0'
3675 ! include 'COMMON.CHAIN'
3676 ! include 'COMMON.DERIV'
3677 ! include 'COMMON.GEO'
3678 ! include 'COMMON.LOCAL'
3679 ! include 'COMMON.INTERACT'
3680 ! include 'COMMON.IOUNITS'
3681 ! include 'COMMON.NAMES'
3682 !el real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV !(MAXRES6) maxres6=6*maxres
3683 real(kind=8),dimension(6*nres) :: stochforcvecV !(MAXRES6) maxres6=6*maxres
3684 !el common /stochcalc/ stochforcvec
3685 real(kind=8) :: ddt1,ddt2
3686 integer :: i,j,ind,inres
3688 ! Compute the stochastic forces which contribute to velocity change
3690 call stochastic_force(stochforcvecV)
3697 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3698 ! ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
3699 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
3701 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
3705 d_t(j,0)=d_t_work(j)
3710 d_t(j,i)=d_t_work(ind+j)
3716 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3718 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3721 d_t(j,inres)=d_t_work(ind+j)
3727 end subroutine sd_verlet2_ciccotti
3729 !-----------------------------------------------------------------------------
3731 !-----------------------------------------------------------------------------
3732 subroutine inertia_tensor
3734 ! Calculating the intertia tensor for the entire protein in order to
3735 ! remove the perpendicular components of velocity matrix which cause
3736 ! the molecule to rotate.
3739 ! implicit real*8 (a-h,o-z)
3740 ! include 'DIMENSIONS'
3741 ! include 'COMMON.CONTROL'
3742 ! include 'COMMON.VAR'
3743 ! include 'COMMON.MD'
3744 ! include 'COMMON.CHAIN'
3745 ! include 'COMMON.DERIV'
3746 ! include 'COMMON.GEO'
3747 ! include 'COMMON.LOCAL'
3748 ! include 'COMMON.INTERACT'
3749 ! include 'COMMON.IOUNITS'
3750 ! include 'COMMON.NAMES'
3752 real(kind=8),dimension(3,3) :: Im,Imcp,eigvec,Id
3753 real(kind=8),dimension(3) :: pr,eigval,L,vp,vrot
3754 real(kind=8) :: M_SC,mag,mag2,M_PEP
3755 real(kind=8),dimension(3,0:nres) :: vpp !(3,0:MAXRES)
3756 real(kind=8),dimension(3) :: vs_p,pp,incr,v
3757 real(kind=8),dimension(3,3) :: pr1,pr2
3759 !el common /gucio/ cm
3760 integer :: iti,inres,i,j,k,mnum
3771 ! calculating the center of the mass of the protein
3775 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
3776 if (itype(i,mnum).eq.ntyp1_molec(mnum)) cycle
3777 M_PEP=M_PEP+mp(mnum)
3779 cm(j)=cm(j)+(c(j,i)+0.5d0*dc(j,i))*mp(mnum)
3788 if (mnum.eq.5) cycle
3789 iti=iabs(itype(i,mnum))
3790 M_SC=M_SC+msc(iabs(iti),mnum)
3793 cm(j)=cm(j)+msc(iabs(iti),mnum)*c(j,inres)
3797 cm(j)=cm(j)/(M_SC+M_PEP)
3802 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
3804 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
3806 Im(1,1)=Im(1,1)+mp(mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
3807 Im(1,2)=Im(1,2)-mp(mnum)*pr(1)*pr(2)
3808 Im(1,3)=Im(1,3)-mp(mnum)*pr(1)*pr(3)
3809 Im(2,3)=Im(2,3)-mp(mnum)*pr(2)*pr(3)
3810 Im(2,2)=Im(2,2)+mp(mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
3811 Im(3,3)=Im(3,3)+mp(mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
3816 iti=iabs(itype(i,mnum))
3817 if (mnum.eq.5) cycle
3820 pr(j)=c(j,inres)-cm(j)
3822 Im(1,1)=Im(1,1)+msc(iabs(iti),mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
3823 Im(1,2)=Im(1,2)-msc(iabs(iti),mnum)*pr(1)*pr(2)
3824 Im(1,3)=Im(1,3)-msc(iabs(iti),mnum)*pr(1)*pr(3)
3825 Im(2,3)=Im(2,3)-msc(iabs(iti),mnum)*pr(2)*pr(3)
3826 Im(2,2)=Im(2,2)+msc(iabs(iti),mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
3827 Im(3,3)=Im(3,3)+msc(iabs(iti),mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
3832 Im(1,1)=Im(1,1)+Ip(mnum)*(1-dc_norm(1,i)*dc_norm(1,i))* &
3833 vbld(i+1)*vbld(i+1)*0.25d0
3834 Im(1,2)=Im(1,2)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(2,i))* &
3835 vbld(i+1)*vbld(i+1)*0.25d0
3836 Im(1,3)=Im(1,3)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(3,i))* &
3837 vbld(i+1)*vbld(i+1)*0.25d0
3838 Im(2,3)=Im(2,3)+Ip(mnum)*(-dc_norm(2,i)*dc_norm(3,i))* &
3839 vbld(i+1)*vbld(i+1)*0.25d0
3840 Im(2,2)=Im(2,2)+Ip(mnum)*(1-dc_norm(2,i)*dc_norm(2,i))* &
3841 vbld(i+1)*vbld(i+1)*0.25d0
3842 Im(3,3)=Im(3,3)+Ip(mnum)*(1-dc_norm(3,i)*dc_norm(3,i))* &
3843 vbld(i+1)*vbld(i+1)*0.25d0
3849 ! if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
3850 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3851 .and.(mnum.ne.5)) then
3852 iti=iabs(itype(i,mnum))
3854 Im(1,1)=Im(1,1)+Isc(iti,mnum)*(1-dc_norm(1,inres)* &
3855 dc_norm(1,inres))*vbld(inres)*vbld(inres)
3856 Im(1,2)=Im(1,2)-Isc(iti,mnum)*(dc_norm(1,inres)* &
3857 dc_norm(2,inres))*vbld(inres)*vbld(inres)
3858 Im(1,3)=Im(1,3)-Isc(iti,mnum)*(dc_norm(1,inres)* &
3859 dc_norm(3,inres))*vbld(inres)*vbld(inres)
3860 Im(2,3)=Im(2,3)-Isc(iti,mnum)*(dc_norm(2,inres)* &
3861 dc_norm(3,inres))*vbld(inres)*vbld(inres)
3862 Im(2,2)=Im(2,2)+Isc(iti,mnum)*(1-dc_norm(2,inres)* &
3863 dc_norm(2,inres))*vbld(inres)*vbld(inres)
3864 Im(3,3)=Im(3,3)+Isc(iti,mnum)*(1-dc_norm(3,inres)* &
3865 dc_norm(3,inres))*vbld(inres)*vbld(inres)
3870 ! write(iout,*) "The angular momentum before adjustment:"
3871 ! write(iout,*) (L(j),j=1,3)
3877 ! Copying the Im matrix for the djacob subroutine
3885 ! Finding the eigenvectors and eignvalues of the inertia tensor
3886 call djacob(3,3,10000,1.0d-10,Imcp,eigvec,eigval)
3887 ! write (iout,*) "Eigenvalues & Eigenvectors"
3888 ! write (iout,'(5x,3f10.5)') (eigval(i),i=1,3)
3891 ! write (iout,'(i5,3f10.5)') i,(eigvec(i,j),j=1,3)
3893 ! Constructing the diagonalized matrix
3895 if (dabs(eigval(i)).gt.1.0d-15) then
3896 Id(i,i)=1.0d0/eigval(i)
3903 Imcp(i,j)=eigvec(j,i)
3909 pr1(i,j)=pr1(i,j)+Id(i,k)*Imcp(k,j)
3916 pr2(i,j)=pr2(i,j)+eigvec(i,k)*pr1(k,j)
3920 ! Calculating the total rotational velocity of the molecule
3923 vrot(i)=vrot(i)+pr2(i,j)*L(j)
3926 ! Resetting the velocities
3928 call vecpr(vrot(1),dc(1,i),vp)
3930 d_t(j,i)=d_t(j,i)-vp(j)
3935 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3936 .and.(mnum.ne.5)) then
3937 ! if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
3938 ! if(itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3940 call vecpr(vrot(1),dc(1,inres),vp)
3942 d_t(j,inres)=d_t(j,inres)-vp(j)
3947 ! write(iout,*) "The angular momentum after adjustment:"
3948 ! write(iout,*) (L(j),j=1,3)
3951 end subroutine inertia_tensor
3952 !-----------------------------------------------------------------------------
3953 subroutine angmom(cm,L)
3956 ! implicit real*8 (a-h,o-z)
3957 ! include 'DIMENSIONS'
3958 ! include 'COMMON.CONTROL'
3959 ! include 'COMMON.VAR'
3960 ! include 'COMMON.MD'
3961 ! include 'COMMON.CHAIN'
3962 ! include 'COMMON.DERIV'
3963 ! include 'COMMON.GEO'
3964 ! include 'COMMON.LOCAL'
3965 ! include 'COMMON.INTERACT'
3966 ! include 'COMMON.IOUNITS'
3967 ! include 'COMMON.NAMES'
3968 real(kind=8) :: mscab
3969 real(kind=8),dimension(3) :: L,cm,pr,vp,vrot,incr,v,pp
3970 integer :: iti,inres,i,j,mnum
3971 ! Calculate the angular momentum
3980 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
3982 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
3985 v(j)=incr(j)+0.5d0*d_t(j,i)
3988 incr(j)=incr(j)+d_t(j,i)
3990 call vecpr(pr(1),v(1),vp)
3992 L(j)=L(j)+mp(mnum)*vp(j)
3996 pp(j)=0.5d0*d_t(j,i)
3998 call vecpr(pr(1),pp(1),vp)
4000 L(j)=L(j)+Ip(mnum)*vp(j)
4008 iti=iabs(itype(i,mnum))
4016 pr(j)=c(j,inres)-cm(j)
4018 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4019 .and.(mnum.ne.5)) then
4021 v(j)=incr(j)+d_t(j,inres)
4028 call vecpr(pr(1),v(1),vp)
4029 ! write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),
4030 ! & " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
4032 L(j)=L(j)+mscab*vp(j)
4034 ! write (iout,*) "L",(l(j),j=1,3)
4035 if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4036 .and.(mnum.ne.5)) then
4038 v(j)=incr(j)+d_t(j,inres)
4040 call vecpr(dc(1,inres),d_t(1,inres),vp)
4042 L(j)=L(j)+Isc(iti,mnum)*vp(j)
4046 incr(j)=incr(j)+d_t(j,i)
4050 end subroutine angmom
4051 !-----------------------------------------------------------------------------
4052 subroutine vcm_vel(vcm)
4055 ! implicit real*8 (a-h,o-z)
4056 ! include 'DIMENSIONS'
4057 ! include 'COMMON.VAR'
4058 ! include 'COMMON.MD'
4059 ! include 'COMMON.CHAIN'
4060 ! include 'COMMON.DERIV'
4061 ! include 'COMMON.GEO'
4062 ! include 'COMMON.LOCAL'
4063 ! include 'COMMON.INTERACT'
4064 ! include 'COMMON.IOUNITS'
4065 real(kind=8),dimension(3) :: vcm,vv
4066 real(kind=8) :: summas,amas
4076 if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
4078 summas=summas+mp(mnum)
4080 vcm(j)=vcm(j)+mp(mnum)*(vv(j)+0.5d0*d_t(j,i))
4084 amas=msc(iabs(itype(i,mnum)),mnum)
4089 if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4090 .and.(mnum.ne.5)) then
4092 vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres))
4096 vcm(j)=vcm(j)+amas*vv(j)
4100 vv(j)=vv(j)+d_t(j,i)
4103 ! write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas
4105 vcm(j)=vcm(j)/summas
4108 end subroutine vcm_vel
4109 !-----------------------------------------------------------------------------
4111 !-----------------------------------------------------------------------------
4113 ! RATTLE algorithm for velocity Verlet - step 1, UNRES
4117 ! implicit real*8 (a-h,o-z)
4118 ! include 'DIMENSIONS'
4120 ! include 'COMMON.CONTROL'
4121 ! include 'COMMON.VAR'
4122 ! include 'COMMON.MD'
4124 ! include 'COMMON.LANGEVIN'
4126 ! include 'COMMON.LANGEVIN.lang0'
4128 ! include 'COMMON.CHAIN'
4129 ! include 'COMMON.DERIV'
4130 ! include 'COMMON.GEO'
4131 ! include 'COMMON.LOCAL'
4132 ! include 'COMMON.INTERACT'
4133 ! include 'COMMON.IOUNITS'
4134 ! include 'COMMON.NAMES'
4135 ! include 'COMMON.TIME1'
4136 !el real(kind=8) :: gginv(2*nres,2*nres),&
4137 !el gdc(3,2*nres,2*nres)
4138 real(kind=8) :: dC_uncor(3,2*nres) !,&
4139 !el real(kind=8) :: Cmat(2*nres,2*nres)
4140 real(kind=8) :: x(2*nres),xcorr(3,2*nres) !maxres2=2*maxres
4141 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
4142 !el common /przechowalnia/ nbond
4143 integer :: max_rattle = 5
4144 logical :: lprn = .false., lprn1 = .false., not_done
4145 real(kind=8) :: tol_rattle = 1.0d-5
4147 integer :: ii,i,j,jj,l,ind,ind1,nres2
4150 !el /common/ przechowalnia
4152 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4153 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4154 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4156 if (lprn) write (iout,*) "RATTLE1"
4160 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4161 .and.(mnum.ne.5)) nbond=nbond+1
4163 ! Make a folded form of the Ginv-matrix
4176 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4181 if (itype(k,1).ne.10 .and. itype(k,mnum).ne.ntyp1_molec(mnum)&
4182 .and.(mnum.ne.5)) then
4186 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4194 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4205 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4209 if (itype(k,1).ne.10) then
4213 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4221 write (iout,*) "Matrix GGinv"
4222 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
4228 if (iter.gt.max_rattle) then
4229 write (iout,*) "Error - too many iterations in RATTLE."
4232 ! Calculate the matrix C = GG**(-1) dC_old o dC
4237 dC_uncor(j,ind1)=dC(j,i)
4241 if (itype(i,1).ne.10) then
4244 dC_uncor(j,ind1)=dC(j,i+nres)
4253 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
4257 if (itype(k,1).ne.10) then
4260 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
4265 ! Calculate deviations from standard virtual-bond lengths
4269 x(ind)=vbld(i+1)**2-vbl**2
4272 if (itype(i,1).ne.10) then
4274 x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4278 write (iout,*) "Coordinates and violations"
4280 write(iout,'(i5,3f10.5,5x,e15.5)') &
4281 i,(dC_uncor(j,i),j=1,3),x(i)
4283 write (iout,*) "Velocities and violations"
4287 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4288 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4292 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4293 .and.(mnum.ne.5)) then
4296 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4297 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4298 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4301 ! write (iout,*) "gdc"
4303 ! write (iout,*) "i",i
4305 ! write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4311 if (dabs(x(i)).gt.xmax) then
4315 if (xmax.lt.tol_rattle) then
4319 ! Calculate the matrix of the system of equations
4324 Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
4329 write (iout,*) "Matrix Cmat"
4330 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4332 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4333 ! Add constraint term to positions
4340 xx = xx+x(ii)*gdc(j,ind,ii)
4344 d_t_new(j,i)=d_t_new(j,i)-xx/d_time
4348 if (itype(i,1).ne.10) then
4353 xx = xx+x(ii)*gdc(j,ind,ii)
4356 dC(j,i+nres)=dC(j,i+nres)-xx
4357 d_t_new(j,i+nres)=d_t_new(j,i+nres)-xx/d_time
4361 ! Rebuild the chain using the new coordinates
4362 call chainbuild_cart
4364 write (iout,*) "New coordinates, Lagrange multipliers,",&
4365 " and differences between actual and standard bond lengths"
4369 xx=vbld(i+1)**2-vbl**2
4370 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4371 i,(dC(j,i),j=1,3),x(ind),xx
4375 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4378 xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4379 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4380 i,(dC(j,i+nres),j=1,3),x(ind),xx
4383 write (iout,*) "Velocities and violations"
4387 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4388 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4391 if (itype(i,1).ne.10) then
4393 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4394 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4395 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4402 10 write (iout,*) "Error - singularity in solving the system",&
4403 " of equations for Lagrange multipliers."
4407 "RATTLE inactive; use -DRATTLE switch at compile time."
4410 end subroutine rattle1
4411 !-----------------------------------------------------------------------------
4413 ! RATTLE algorithm for velocity Verlet - step 2, UNRES
4417 ! implicit real*8 (a-h,o-z)
4418 ! include 'DIMENSIONS'
4420 ! include 'COMMON.CONTROL'
4421 ! include 'COMMON.VAR'
4422 ! include 'COMMON.MD'
4424 ! include 'COMMON.LANGEVIN'
4426 ! include 'COMMON.LANGEVIN.lang0'
4428 ! include 'COMMON.CHAIN'
4429 ! include 'COMMON.DERIV'
4430 ! include 'COMMON.GEO'
4431 ! include 'COMMON.LOCAL'
4432 ! include 'COMMON.INTERACT'
4433 ! include 'COMMON.IOUNITS'
4434 ! include 'COMMON.NAMES'
4435 ! include 'COMMON.TIME1'
4436 !el real(kind=8) :: gginv(2*nres,2*nres),&
4437 !el gdc(3,2*nres,2*nres)
4438 real(kind=8) :: dC_uncor(3,2*nres) !,&
4439 !el Cmat(2*nres,2*nres)
4440 real(kind=8) :: x(2*nres) !maxres2=2*maxres
4441 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
4442 !el common /przechowalnia/ nbond
4443 integer :: max_rattle = 5
4444 logical :: lprn = .false., lprn1 = .false., not_done
4445 real(kind=8) :: tol_rattle = 1.0d-5
4449 !el /common/ przechowalnia
4450 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4451 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4452 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4454 if (lprn) write (iout,*) "RATTLE2"
4455 if (lprn) write (iout,*) "Velocity correction"
4456 ! Calculate the matrix G dC
4462 gdc(j,i,ind)=GGinv(i,ind)*dC(j,k)
4467 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4468 .and.(mnum.ne.5)) then
4471 gdc(j,i,ind)=GGinv(i,ind)*dC(j,k+nres)
4477 ! write (iout,*) "gdc"
4479 ! write (iout,*) "i",i
4481 ! write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4485 ! Calculate the matrix of the system of equations
4492 Cmat(ind,j)=Cmat(ind,j)+dC(k,i)*gdc(k,ind,j)
4498 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4499 .and.(mnum.ne.5)) then
4504 Cmat(ind,j)=Cmat(ind,j)+dC(k,i+nres)*gdc(k,ind,j)
4509 ! Calculate the scalar product dC o d_t_new
4513 x(ind)=scalar(d_t(1,i),dC(1,i))
4517 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4518 .and.(mnum.ne.5)) then
4520 x(ind)=scalar(d_t(1,i+nres),dC(1,i+nres))
4524 write (iout,*) "Velocities and violations"
4528 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4529 i,ind,(d_t(j,i),j=1,3),x(ind)
4533 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4534 .and.(mnum.ne.5)) then
4536 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4537 i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind)
4543 if (dabs(x(i)).gt.xmax) then
4547 if (xmax.lt.tol_rattle) then
4552 write (iout,*) "Matrix Cmat"
4553 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4555 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4556 ! Add constraint term to velocities
4563 xx = xx+x(ii)*gdc(j,ind,ii)
4565 d_t(j,i)=d_t(j,i)-xx
4570 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4571 .and.(mnum.ne.5)) then
4576 xx = xx+x(ii)*gdc(j,ind,ii)
4578 d_t(j,i+nres)=d_t(j,i+nres)-xx
4584 "New velocities, Lagrange multipliers violations"
4588 if (lprn) write (iout,'(2i5,3f10.5,5x,2e15.5)') &
4589 i,ind,(d_t(j,i),j=1,3),x(ind),scalar(d_t(1,i),dC(1,i))
4593 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4596 write (iout,'(2i5,3f10.5,5x,2e15.5)') &
4597 i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind),&
4598 scalar(d_t(1,i+nres),dC(1,i+nres))
4604 10 write (iout,*) "Error - singularity in solving the system",&
4605 " of equations for Lagrange multipliers."
4609 "RATTLE inactive; use -DRATTLE option at compile time."
4612 end subroutine rattle2
4613 !-----------------------------------------------------------------------------
4614 subroutine rattle_brown
4615 ! RATTLE/LINCS algorithm for Brownian dynamics, UNRES
4619 ! implicit real*8 (a-h,o-z)
4620 ! include 'DIMENSIONS'
4622 ! include 'COMMON.CONTROL'
4623 ! include 'COMMON.VAR'
4624 ! include 'COMMON.MD'
4626 ! include 'COMMON.LANGEVIN'
4628 ! include 'COMMON.LANGEVIN.lang0'
4630 ! include 'COMMON.CHAIN'
4631 ! include 'COMMON.DERIV'
4632 ! include 'COMMON.GEO'
4633 ! include 'COMMON.LOCAL'
4634 ! include 'COMMON.INTERACT'
4635 ! include 'COMMON.IOUNITS'
4636 ! include 'COMMON.NAMES'
4637 ! include 'COMMON.TIME1'
4638 !el real(kind=8) :: gginv(2*nres,2*nres),&
4639 !el gdc(3,2*nres,2*nres)
4640 real(kind=8) :: dC_uncor(3,2*nres) !,&
4641 !el real(kind=8) :: Cmat(2*nres,2*nres)
4642 real(kind=8) :: x(2*nres) !maxres2=2*maxres
4643 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
4644 !el common /przechowalnia/ nbond
4645 integer :: max_rattle = 5
4646 logical :: lprn = .false., lprn1 = .false., not_done
4647 real(kind=8) :: tol_rattle = 1.0d-5
4651 !el /common/ przechowalnia
4652 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4653 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4654 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4657 if (lprn) write (iout,*) "RATTLE_BROWN"
4660 if (itype(i,1).ne.10) nbond=nbond+1
4662 ! Make a folded form of the Ginv-matrix
4675 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4679 if (itype(k,1).ne.10) then
4683 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4690 if (itype(i,1).ne.10) then
4700 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4704 if (itype(k,1).ne.10) then
4708 if (j.eq.1 .and. l.eq.1)GGinv(ii,jj)=fricmat(ind,ind1)
4716 write (iout,*) "Matrix GGinv"
4717 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
4723 if (iter.gt.max_rattle) then
4724 write (iout,*) "Error - too many iterations in RATTLE."
4727 ! Calculate the matrix C = GG**(-1) dC_old o dC
4732 dC_uncor(j,ind1)=dC(j,i)
4736 if (itype(i,1).ne.10) then
4739 dC_uncor(j,ind1)=dC(j,i+nres)
4748 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
4752 if (itype(k,1).ne.10) then
4755 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
4760 ! Calculate deviations from standard virtual-bond lengths
4764 x(ind)=vbld(i+1)**2-vbl**2
4767 if (itype(i,1).ne.10) then
4769 x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4773 write (iout,*) "Coordinates and violations"
4775 write(iout,'(i5,3f10.5,5x,e15.5)') &
4776 i,(dC_uncor(j,i),j=1,3),x(i)
4778 write (iout,*) "Velocities and violations"
4782 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4783 i,ind,(d_t(j,i),j=1,3),scalar(d_t(1,i),dC_old(1,i))
4786 if (itype(i,1).ne.10) then
4788 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4789 i+nres,ind,(d_t(j,i+nres),j=1,3),&
4790 scalar(d_t(1,i+nres),dC_old(1,i+nres))
4793 write (iout,*) "gdc"
4795 write (iout,*) "i",i
4797 write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4803 if (dabs(x(i)).gt.xmax) then
4807 if (xmax.lt.tol_rattle) then
4811 ! Calculate the matrix of the system of equations
4816 Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
4821 write (iout,*) "Matrix Cmat"
4822 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4824 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4825 ! Add constraint term to positions
4832 xx = xx+x(ii)*gdc(j,ind,ii)
4835 d_t(j,i)=d_t(j,i)+xx/d_time
4840 if (itype(i,1).ne.10) then
4845 xx = xx+x(ii)*gdc(j,ind,ii)
4848 d_t(j,i+nres)=d_t(j,i+nres)+xx/d_time
4849 dC(j,i+nres)=dC(j,i+nres)+xx
4853 ! Rebuild the chain using the new coordinates
4854 call chainbuild_cart
4856 write (iout,*) "New coordinates, Lagrange multipliers,",&
4857 " and differences between actual and standard bond lengths"
4861 xx=vbld(i+1)**2-vbl**2
4862 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4863 i,(dC(j,i),j=1,3),x(ind),xx
4866 if (itype(i,1).ne.10) then
4868 xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4869 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4870 i,(dC(j,i+nres),j=1,3),x(ind),xx
4873 write (iout,*) "Velocities and violations"
4877 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4878 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4881 if (itype(i,1).ne.10) then
4883 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4884 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4885 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4892 10 write (iout,*) "Error - singularity in solving the system",&
4893 " of equations for Lagrange multipliers."
4897 "RATTLE inactive; use -DRATTLE option at compile time"
4900 end subroutine rattle_brown
4901 !-----------------------------------------------------------------------------
4903 !-----------------------------------------------------------------------------
4904 subroutine friction_force
4909 ! implicit real*8 (a-h,o-z)
4910 ! include 'DIMENSIONS'
4911 ! include 'COMMON.VAR'
4912 ! include 'COMMON.CHAIN'
4913 ! include 'COMMON.DERIV'
4914 ! include 'COMMON.GEO'
4915 ! include 'COMMON.LOCAL'
4916 ! include 'COMMON.INTERACT'
4917 ! include 'COMMON.MD'
4919 ! include 'COMMON.LANGEVIN'
4921 ! include 'COMMON.LANGEVIN.lang0'
4923 ! include 'COMMON.IOUNITS'
4924 !el real(kind=8),dimension(6*nres) :: gamvec !(MAXRES6) maxres6=6*maxres
4925 !el common /syfek/ gamvec
4926 real(kind=8) :: vv(3),vvtot(3,nres),v_work(6*nres) !,&
4927 !el ginvfric(2*nres,2*nres) !maxres2=2*maxres
4928 !el common /przechowalnia/ ginvfric
4930 logical :: lprn = .false., checkmode = .false.
4931 integer :: i,j,ind,k,nres2,nres6,mnum
4935 if(.not.allocated(gamvec)) allocate(gamvec(nres6)) !(MAXRES6)
4936 if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
4944 d_t_work(j)=d_t(j,0)
4949 d_t_work(ind+j)=d_t(j,i)
4955 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
4956 .and.(mnum.ne.5)) then
4958 d_t_work(ind+j)=d_t(j,i+nres)
4964 call fricmat_mult(d_t_work,fric_work)
4966 if (.not.checkmode) return
4969 write (iout,*) "d_t_work and fric_work"
4971 write (iout,'(i3,2e15.5)') i,d_t_work(i),fric_work(i)
4975 friction(j,0)=fric_work(j)
4980 friction(j,i)=fric_work(ind+j)
4986 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
4987 .and.(mnum.ne.5)) then
4988 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
4990 friction(j,i+nres)=fric_work(ind+j)
4996 write(iout,*) "Friction backbone"
4998 write(iout,'(i5,3e15.5,5x,3e15.5)') &
4999 i,(friction(j,i),j=1,3),(d_t(j,i),j=1,3)
5001 write(iout,*) "Friction side chain"
5003 write(iout,'(i5,3e15.5,5x,3e15.5)') &
5004 i,(friction(j,i+nres),j=1,3),(d_t(j,i+nres),j=1,3)
5013 vvtot(j,i)=vv(j)+0.5d0*d_t(j,i)
5014 vvtot(j,i+nres)=vv(j)+d_t(j,i+nres)
5015 vv(j)=vv(j)+d_t(j,i)
5018 write (iout,*) "vvtot backbone and sidechain"
5020 write (iout,'(i5,3e15.5,5x,3e15.5)') i,(vvtot(j,i),j=1,3),&
5021 (vvtot(j,i+nres),j=1,3)
5026 v_work(ind+j)=vvtot(j,i)
5032 v_work(ind+j)=vvtot(j,i+nres)
5036 write (iout,*) "v_work gamvec and site-based friction forces"
5038 write (iout,'(i5,3e15.5)') i,v_work(i),gamvec(i),&
5042 ! fric_work1(i)=0.0d0
5044 ! fric_work1(i)=fric_work1(i)-A(j,i)*gamvec(j)*v_work(j)
5047 ! write (iout,*) "fric_work and fric_work1"
5049 ! write (iout,'(i5,2e15.5)') i,fric_work(i),fric_work1(i)
5055 ginvfric(i,j)=ginvfric(i,j)+ginv(i,k)*fricmat(k,j)
5059 write (iout,*) "ginvfric"
5061 write (iout,'(i5,100f8.3)') i,(ginvfric(i,j),j=1,dimen)
5063 write (iout,*) "symmetry check"
5066 write (iout,*) i,j,ginvfric(i,j)-ginvfric(j,i)
5071 end subroutine friction_force
5072 !-----------------------------------------------------------------------------
5073 subroutine setup_fricmat
5077 use control_data, only:time_Bcast
5078 use control, only:tcpu
5080 ! implicit real*8 (a-h,o-z)
5084 real(kind=8) :: time00
5086 ! include 'DIMENSIONS'
5087 ! include 'COMMON.VAR'
5088 ! include 'COMMON.CHAIN'
5089 ! include 'COMMON.DERIV'
5090 ! include 'COMMON.GEO'
5091 ! include 'COMMON.LOCAL'
5092 ! include 'COMMON.INTERACT'
5093 ! include 'COMMON.MD'
5094 ! include 'COMMON.SETUP'
5095 ! include 'COMMON.TIME1'
5096 ! integer licznik /0/
5099 ! include 'COMMON.LANGEVIN'
5101 ! include 'COMMON.LANGEVIN.lang0'
5103 ! include 'COMMON.IOUNITS'
5105 integer :: i,j,ind,ind1,m
5106 logical :: lprn = .false.
5107 real(kind=8) :: dtdi !el ,gamvec(2*nres)
5108 !el real(kind=8),dimension(2*nres,2*nres) :: ginvfric,fcopy
5109 ! real(kind=8),allocatable,dimension(:,:) :: fcopy
5110 !el real(kind=8),dimension(2*nres*(2*nres+1)/2) :: Ghalf !(mmaxres2) (mmaxres2=(maxres2*(maxres2+1)/2))
5111 !el common /syfek/ gamvec
5112 real(kind=8) :: work(8*2*nres)
5113 integer :: iwork(2*nres)
5114 !el common /przechowalnia/ ginvfric,Ghalf,fcopy
5115 integer :: ii,iti,k,l,nzero,nres2,nres6,ierr,mnum
5119 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
5120 if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
5121 if (fg_rank.ne.king) goto 10
5126 if(.not.allocated(gamvec)) allocate(gamvec(nres2)) !(MAXRES2)
5127 if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
5128 if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
5129 !el allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
5130 if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !maxres2=2*maxres
5132 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
5133 ! Zeroing out fricmat
5139 ! Load the friction coefficients corresponding to peptide groups
5144 gamvec(ind1)=gamp(mnum)
5146 ! Load the friction coefficients corresponding to side chains
5150 gamsc(ntyp1_molec(j),j)=1.0d0
5157 gamvec(ii)=gamsc(iabs(iti),mnum)
5159 if (surfarea) call sdarea(gamvec)
5161 ! write (iout,*) "Matrix A and vector gamma"
5163 ! write (iout,'(i2,$)') i
5165 ! write (iout,'(f4.1,$)') A(i,j)
5167 ! write (iout,'(f8.3)') gamvec(i)
5171 write (iout,*) "Vector gamvec"
5173 write (iout,'(i5,f10.5)') i, gamvec(i)
5177 ! The friction matrix
5182 dtdi=dtdi+A(j,k)*A(j,i)*gamvec(j)
5189 write (iout,'(//a)') "Matrix fricmat"
5190 call matout2(dimen,dimen,nres2,nres2,fricmat)
5192 if (lang.eq.2 .or. lang.eq.3) then
5193 ! Mass-scale the friction matrix if non-direct integration will be performed
5199 Ginvfric(i,j)=Ginvfric(i,j)+ &
5200 Gsqrm(i,k)*Gsqrm(l,j)*fricmat(k,l)
5205 ! Diagonalize the friction matrix
5210 Ghalf(ind)=Ginvfric(i,j)
5213 call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
5216 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
5217 " mass-scaled friction matrix"
5218 call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
5220 ! Precompute matrices for tinker stochastic integrator
5227 mt1(i,j)=mt1(i,j)+fricvec(k,i)*gsqrm(k,j)
5228 mt2(i,j)=mt2(i,j)+fricvec(k,i)*gsqrp(k,j)
5234 else if (lang.eq.4) then
5235 ! Diagonalize the friction matrix
5240 Ghalf(ind)=fricmat(i,j)
5243 call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
5246 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
5248 call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
5250 ! Determine the number of zero eigenvalues of the friction matrix
5251 nzero=max0(dimen-dimen1,0)
5252 ! do while (fricgam(nzero+1).le.1.0d-5 .and. nzero.lt.dimen)
5255 write (iout,*) "Number of zero eigenvalues:",nzero
5260 fricmat(i,j)=fricmat(i,j) &
5261 +fricvec(i,k)*fricvec(j,k)/fricgam(k)
5266 write (iout,'(//a)') "Generalized inverse of fricmat"
5267 call matout(dimen,dimen,nres6,nres6,fricmat)
5272 if (nfgtasks.gt.1) then
5273 if (fg_rank.eq.0) then
5274 ! The matching BROADCAST for fg processors is called in ERGASTULUM
5280 call MPI_Bcast(10,1,MPI_INTEGER,king,FG_COMM,IERROR)
5282 time_Bcast=time_Bcast+MPI_Wtime()-time00
5284 time_Bcast=time_Bcast+tcpu()-time00
5286 ! print *,"Processor",myrank,
5287 ! & " BROADCAST iorder in SETUP_FRICMAT"
5290 write (iout,*) "setup_fricmat licznik"!,licznik !sp
5296 ! Scatter the friction matrix
5297 call MPI_Scatterv(fricmat(1,1),nginv_counts(0),&
5298 nginv_start(0),MPI_DOUBLE_PRECISION,fcopy(1,1),&
5299 myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
5302 time_scatter=time_scatter+MPI_Wtime()-time00
5303 time_scatter_fmat=time_scatter_fmat+MPI_Wtime()-time00
5305 time_scatter=time_scatter+tcpu()-time00
5306 time_scatter_fmat=time_scatter_fmat+tcpu()-time00
5310 do j=1,2*my_ng_count
5311 fricmat(j,i)=fcopy(i,j)
5314 ! write (iout,*) "My chunk of fricmat"
5315 ! call MATOUT2(my_ng_count,dimen,maxres2,maxres2,fcopy)
5319 end subroutine setup_fricmat
5320 !-----------------------------------------------------------------------------
5321 subroutine stochastic_force(stochforcvec)
5324 use random, only:anorm_distr
5325 ! implicit real*8 (a-h,o-z)
5326 ! include 'DIMENSIONS'
5327 use control, only: tcpu
5332 ! include 'COMMON.VAR'
5333 ! include 'COMMON.CHAIN'
5334 ! include 'COMMON.DERIV'
5335 ! include 'COMMON.GEO'
5336 ! include 'COMMON.LOCAL'
5337 ! include 'COMMON.INTERACT'
5338 ! include 'COMMON.MD'
5339 ! include 'COMMON.TIME1'
5341 ! include 'COMMON.LANGEVIN'
5343 ! include 'COMMON.LANGEVIN.lang0'
5345 ! include 'COMMON.IOUNITS'
5347 real(kind=8) :: x,sig,lowb,highb
5348 real(kind=8) :: ff(3),force(3,0:2*nres),zeta2,lowb2
5349 real(kind=8) :: highb2,sig2,forcvec(6*nres),stochforcvec(6*nres)
5350 real(kind=8) :: time00
5351 logical :: lprn = .false.
5352 integer :: i,j,ind,mnum
5356 stochforc(j,i)=0.0d0
5366 ! Compute the stochastic forces acting on bodies. Store in force.
5372 force(j,i)=anorm_distr(x,sig,lowb,highb)
5380 force(j,i+nres)=anorm_distr(x,sig2,lowb2,highb2)
5384 time_fsample=time_fsample+MPI_Wtime()-time00
5386 time_fsample=time_fsample+tcpu()-time00
5388 ! Compute the stochastic forces acting on virtual-bond vectors.
5394 stochforc(j,i)=ff(j)+0.5d0*force(j,i)
5397 ff(j)=ff(j)+force(j,i)
5399 ! if (itype(i+1,1).ne.ntyp1) then
5401 if (itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
5403 stochforc(j,i)=stochforc(j,i)+force(j,i+nres+1)
5404 ff(j)=ff(j)+force(j,i+nres+1)
5409 stochforc(j,0)=ff(j)+force(j,nnt+nres)
5413 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5414 .and.(mnum.ne.5)) then
5415 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
5417 stochforc(j,i+nres)=force(j,i+nres)
5423 stochforcvec(j)=stochforc(j,0)
5428 stochforcvec(ind+j)=stochforc(j,i)
5434 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5435 .and.(mnum.ne.5)) then
5436 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
5438 stochforcvec(ind+j)=stochforc(j,i+nres)
5444 write (iout,*) "stochforcvec"
5446 write(iout,'(i5,e15.5)') i,stochforcvec(i)
5448 write(iout,*) "Stochastic forces backbone"
5450 write(iout,'(i5,3e15.5)') i,(stochforc(j,i),j=1,3)
5452 write(iout,*) "Stochastic forces side chain"
5454 write(iout,'(i5,3e15.5)') &
5455 i,(stochforc(j,i+nres),j=1,3)
5463 write (iout,*) i,ind
5465 forcvec(ind+j)=force(j,i)
5470 write (iout,*) i,ind
5472 forcvec(j+ind)=force(j,i+nres)
5477 write (iout,*) "forcvec"
5481 write (iout,'(2i3,2f10.5)') i,j,force(j,i),&
5488 write (iout,'(2i3,2f10.5)') i,j,force(j,i+nres),&
5497 end subroutine stochastic_force
5498 !-----------------------------------------------------------------------------
5499 subroutine sdarea(gamvec)
5501 ! Scale the friction coefficients according to solvent accessible surface areas
5502 ! Code adapted from TINKER
5506 ! implicit real*8 (a-h,o-z)
5507 ! include 'DIMENSIONS'
5508 ! include 'COMMON.CONTROL'
5509 ! include 'COMMON.VAR'
5510 ! include 'COMMON.MD'
5512 ! include 'COMMON.LANGEVIN'
5514 ! include 'COMMON.LANGEVIN.lang0'
5516 ! include 'COMMON.CHAIN'
5517 ! include 'COMMON.DERIV'
5518 ! include 'COMMON.GEO'
5519 ! include 'COMMON.LOCAL'
5520 ! include 'COMMON.INTERACT'
5521 ! include 'COMMON.IOUNITS'
5522 ! include 'COMMON.NAMES'
5523 real(kind=8),dimension(2*nres) :: radius,gamvec !(maxres2)
5524 real(kind=8),parameter :: twosix = 1.122462048309372981d0
5525 logical :: lprn = .false.
5526 real(kind=8) :: probe,area,ratio
5527 integer :: i,j,ind,iti,mnum
5529 ! determine new friction coefficients every few SD steps
5531 ! set the atomic radii to estimates of sigma values
5533 ! print *,"Entered sdarea"
5539 ! Load peptide group radii
5542 radius(i)=pstok(mnum)
5544 ! Load side chain radii
5548 radius(i+nres)=restok(iti,mnum)
5551 ! write (iout,*) "i",i," radius",radius(i)
5554 radius(i) = radius(i) / twosix
5555 if (radius(i) .ne. 0.0d0) radius(i) = radius(i) + probe
5558 ! scale atomic friction coefficients by accessible area
5560 if (lprn) write (iout,*) &
5561 "Original gammas, surface areas, scaling factors, new gammas, ",&
5562 "std's of stochastic forces"
5565 if (radius(i).gt.0.0d0) then
5566 call surfatom (i,area,radius)
5567 ratio = dmax1(area/(4.0d0*pi*radius(i)**2),1.0d-1)
5568 if (lprn) write (iout,'(i5,3f10.5,$)') &
5569 i,gamvec(ind+1),area,ratio
5572 gamvec(ind) = ratio * gamvec(ind)
5575 stdforcp(i)=stdfp(mnum)*dsqrt(gamvec(ind))
5576 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcp(i)
5580 if (radius(i+nres).gt.0.0d0) then
5581 call surfatom (i+nres,area,radius)
5582 ratio = dmax1(area/(4.0d0*pi*radius(i+nres)**2),1.0d-1)
5583 if (lprn) write (iout,'(i5,3f10.5,$)') &
5584 i,gamvec(ind+1),area,ratio
5587 gamvec(ind) = ratio * gamvec(ind)
5590 stdforcsc(i)=stdfsc(itype(i,mnum),mnum)*dsqrt(gamvec(ind))
5591 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcsc(i)
5596 end subroutine sdarea
5597 !-----------------------------------------------------------------------------
5599 !-----------------------------------------------------------------------------
5602 ! ###################################################
5603 ! ## COPYRIGHT (C) 1996 by Jay William Ponder ##
5604 ! ## All Rights Reserved ##
5605 ! ###################################################
5607 ! ################################################################
5609 ! ## subroutine surfatom -- exposed surface area of an atom ##
5611 ! ################################################################
5614 ! "surfatom" performs an analytical computation of the surface
5615 ! area of a specified atom; a simplified version of "surface"
5617 ! literature references:
5619 ! T. J. Richmond, "Solvent Accessible Surface Area and
5620 ! Excluded Volume in Proteins", Journal of Molecular Biology,
5623 ! L. Wesson and D. Eisenberg, "Atomic Solvation Parameters
5624 ! Applied to Molecular Dynamics of Proteins in Solution",
5625 ! Protein Science, 1, 227-235 (1992)
5627 ! variables and parameters:
5629 ! ir number of atom for which area is desired
5630 ! area accessible surface area of the atom
5631 ! radius radii of each of the individual atoms
5634 subroutine surfatom(ir,area,radius)
5636 ! implicit real*8 (a-h,o-z)
5637 ! include 'DIMENSIONS'
5639 ! include 'COMMON.GEO'
5640 ! include 'COMMON.IOUNITS'
5642 integer :: nsup,nstart_sup
5643 ! double precision c,dc,dc_old,d_c_work,xloc,xrot,dc_norm
5644 ! common /chain/ c(3,maxres2+2),dc(3,0:maxres2),dc_old(3,0:maxres2),
5645 ! & xloc(3,maxres),xrot(3,maxres),dc_norm(3,0:maxres2),
5646 ! & dc_work(MAXRES6),nres,nres0
5647 integer,parameter :: maxarc=300
5651 integer :: mi,ni,narc
5652 integer :: key(maxarc)
5653 integer :: intag(maxarc)
5654 integer :: intag1(maxarc)
5655 real(kind=8) :: area,arcsum
5656 real(kind=8) :: arclen,exang
5657 real(kind=8) :: delta,delta2
5658 real(kind=8) :: eps,rmove
5659 real(kind=8) :: xr,yr,zr
5660 real(kind=8) :: rr,rrsq
5661 real(kind=8) :: rplus,rminus
5662 real(kind=8) :: axx,axy,axz
5663 real(kind=8) :: ayx,ayy
5664 real(kind=8) :: azx,azy,azz
5665 real(kind=8) :: uxj,uyj,uzj
5666 real(kind=8) :: tx,ty,tz
5667 real(kind=8) :: txb,tyb,td
5668 real(kind=8) :: tr2,tr,txr,tyr
5669 real(kind=8) :: tk1,tk2
5670 real(kind=8) :: thec,the,t,tb
5671 real(kind=8) :: txk,tyk,tzk
5672 real(kind=8) :: t1,ti,tf,tt
5673 real(kind=8) :: txj,tyj,tzj
5674 real(kind=8) :: ccsq,cc,xysq
5675 real(kind=8) :: bsqk,bk,cosine
5676 real(kind=8) :: dsqj,gi,pix2
5677 real(kind=8) :: therk,dk,gk
5678 real(kind=8) :: risqk,rik
5679 real(kind=8) :: radius(2*nres) !(maxatm) (maxatm=maxres2)
5680 real(kind=8) :: ri(maxarc),risq(maxarc)
5681 real(kind=8) :: ux(maxarc),uy(maxarc),uz(maxarc)
5682 real(kind=8) :: xc(maxarc),yc(maxarc),zc(maxarc)
5683 real(kind=8) :: xc1(maxarc),yc1(maxarc),zc1(maxarc)
5684 real(kind=8) :: dsq(maxarc),bsq(maxarc)
5685 real(kind=8) :: dsq1(maxarc),bsq1(maxarc)
5686 real(kind=8) :: arci(maxarc),arcf(maxarc)
5687 real(kind=8) :: ex(maxarc),lt(maxarc),gr(maxarc)
5688 real(kind=8) :: b(maxarc),b1(maxarc),bg(maxarc)
5689 real(kind=8) :: kent(maxarc),kout(maxarc)
5690 real(kind=8) :: ther(maxarc)
5691 logical :: moved,top
5692 logical :: omit(maxarc)
5695 ! maxatm = 2*nres !maxres2 maxres2=2*maxres
5696 ! maxlight = 8*maxatm
5699 ! maxtors = 4*maxatm
5701 ! zero out the surface area for the sphere of interest
5704 ! write (2,*) "ir",ir," radius",radius(ir)
5705 if (radius(ir) .eq. 0.0d0) return
5707 ! set the overlap significance and connectivity shift
5711 delta2 = delta * delta
5716 ! store coordinates and radius of the sphere of interest
5724 ! initialize values of some counters and summations
5733 ! test each sphere to see if it overlaps the sphere of interest
5736 if (i.eq.ir .or. radius(i).eq.0.0d0) goto 30
5737 rplus = rr + radius(i)
5739 if (abs(tx) .ge. rplus) goto 30
5741 if (abs(ty) .ge. rplus) goto 30
5743 if (abs(tz) .ge. rplus) goto 30
5745 ! check for sphere overlap by testing distance against radii
5747 xysq = tx*tx + ty*ty
5748 if (xysq .lt. delta2) then
5755 if (rplus-cc .le. delta) goto 30
5756 rminus = rr - radius(i)
5758 ! check to see if sphere of interest is completely buried
5760 if (cc-abs(rminus) .le. delta) then
5761 if (rminus .le. 0.0d0) goto 170
5765 ! check for too many overlaps with sphere of interest
5767 if (io .ge. maxarc) then
5769 20 format (/,' SURFATOM -- Increase the Value of MAXARC')
5773 ! get overlap between current sphere and sphere of interest
5782 gr(io) = (ccsq+rplus*rminus) / (2.0d0*rr*b1(io))
5788 ! case where no other spheres overlap the sphere of interest
5791 area = 4.0d0 * pi * rrsq
5795 ! case where only one sphere overlaps the sphere of interest
5798 area = pix2 * (1.0d0 + gr(1))
5799 area = mod(area,4.0d0*pi) * rrsq
5803 ! case where many spheres intersect the sphere of interest;
5804 ! sort the intersecting spheres by their degree of overlap
5806 call sort2 (io,gr,key)
5809 intag(i) = intag1(k)
5818 ! get radius of each overlap circle on surface of the sphere
5823 risq(i) = rrsq - gi*gi
5824 ri(i) = sqrt(risq(i))
5825 ther(i) = 0.5d0*pi - asin(min(1.0d0,max(-1.0d0,gr(i))))
5828 ! find boundary of inaccessible area on sphere of interest
5831 if (.not. omit(k)) then
5838 ! check to see if J circle is intersecting K circle;
5839 ! get distance between circle centers and sum of radii
5842 if (omit(j)) goto 60
5843 cc = (txk*xc(j)+tyk*yc(j)+tzk*zc(j))/(bk*b(j))
5844 cc = acos(min(1.0d0,max(-1.0d0,cc)))
5845 td = therk + ther(j)
5847 ! check to see if circles enclose separate regions
5849 if (cc .ge. td) goto 60
5851 ! check for circle J completely inside circle K
5853 if (cc+ther(j) .lt. therk) goto 40
5855 ! check for circles that are essentially parallel
5857 if (cc .gt. delta) goto 50
5862 ! check to see if sphere of interest is completely buried
5865 if (pix2-cc .le. td) goto 170
5871 ! find T value of circle intersections
5874 if (omit(k)) goto 110
5889 ! rotation matrix elements
5901 if (.not. omit(j)) then
5906 ! rotate spheres so K vector colinear with z-axis
5908 uxj = txj*axx + tyj*axy - tzj*axz
5909 uyj = tyj*ayy - txj*ayx
5910 uzj = txj*azx + tyj*azy + tzj*azz
5911 cosine = min(1.0d0,max(-1.0d0,uzj/b(j)))
5912 if (acos(cosine) .lt. therk+ther(j)) then
5913 dsqj = uxj*uxj + uyj*uyj
5918 tr2 = risqk*dsqj - tb*tb
5924 ! get T values of intersection for K circle
5927 tb = min(1.0d0,max(-1.0d0,tb))
5929 if (tyb-txr .lt. 0.0d0) tk1 = pix2 - tk1
5931 tb = min(1.0d0,max(-1.0d0,tb))
5933 if (tyb+txr .lt. 0.0d0) tk2 = pix2 - tk2
5934 thec = (rrsq*uzj-gk*bg(j)) / (rik*ri(j)*b(j))
5935 if (abs(thec) .lt. 1.0d0) then
5937 else if (thec .ge. 1.0d0) then
5939 else if (thec .le. -1.0d0) then
5943 ! see if "tk1" is entry or exit point; check t=0 point;
5944 ! "ti" is exit point, "tf" is entry point
5946 cosine = min(1.0d0,max(-1.0d0, &
5947 (uzj*gk-uxj*rik)/(b(j)*rr)))
5948 if ((acos(cosine)-ther(j))*(tk2-tk1) .le. 0.0d0) then
5956 if (narc .ge. maxarc) then
5958 70 format (/,' SURFATOM -- Increase the Value',&
5962 if (tf .le. ti) then
5983 ! special case; K circle without intersections
5985 if (narc .le. 0) goto 90
5987 ! general case; sum up arclength and set connectivity code
5989 call sort2 (narc,arci,key)
5994 if (narc .gt. 1) then
5997 if (t .lt. arci(j)) then
5998 arcsum = arcsum + arci(j) - t
5999 exang = exang + ex(ni)
6001 if (jb .ge. maxarc) then
6003 80 format (/,' SURFATOM -- Increase the Value',&
6008 kent(jb) = maxarc*i + k
6010 kout(jb) = maxarc*k + i
6019 arcsum = arcsum + pix2 - t
6021 exang = exang + ex(ni)
6024 kent(jb) = maxarc*i + k
6026 kout(jb) = maxarc*k + i
6033 arclen = arclen + gr(k)*arcsum
6036 if (arclen .eq. 0.0d0) goto 170
6037 if (jb .eq. 0) goto 150
6039 ! find number of independent boundaries and check connectivity
6043 if (kout(k) .ne. 0) then
6050 if (m .eq. kent(ii)) then
6053 if (j .eq. jb) goto 150
6065 ! attempt to fix connectivity error by moving atom slightly
6069 140 format (/,' SURFATOM -- Connectivity Error at Atom',i6)
6078 ! compute the exposed surface area for the sphere of interest
6081 area = ib*pix2 + exang + arclen
6082 area = mod(area,4.0d0*pi) * rrsq
6084 ! attempt to fix negative area by moving atom slightly
6086 if (area .lt. 0.0d0) then
6089 160 format (/,' SURFATOM -- Negative Area at Atom',i6)
6100 end subroutine surfatom
6101 !----------------------------------------------------------------
6102 !----------------------------------------------------------------
6103 subroutine alloc_MD_arrays
6104 !EL Allocation of arrays used by MD module
6106 integer :: nres2,nres6
6109 !----------------------
6113 allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
6114 allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
6115 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
6116 allocate(fricvec(nres2,nres2))
6117 allocate(pfric_mat(nres2,nres2),vfric_mat(nres2,nres2))
6118 allocate(afric_mat(nres2,nres2),prand_mat(nres2,nres2))
6119 allocate(vrand_mat1(nres2,nres2),vrand_mat2(nres2,nres2)) !(MAXRES2,MAXRES2)
6120 allocate(pfric0_mat(nres2,nres2,0:maxflag_stoch))
6121 allocate(afric0_mat(nres2,nres2,0:maxflag_stoch))
6122 allocate(vfric0_mat(nres2,nres2,0:maxflag_stoch))
6123 allocate(prand0_mat(nres2,nres2,0:maxflag_stoch))
6124 allocate(vrand0_mat1(nres2,nres2,0:maxflag_stoch))
6125 allocate(vrand0_mat2(nres2,nres2,0:maxflag_stoch)) !(MAXRES2,MAXRES2,0:maxflag_stoch)
6126 allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
6128 allocate(mt1(nres2,nres2),mt2(nres2,nres2),mt3(nres2,nres2)) !(maxres2,maxres2)
6129 !----------------------
6131 ! commom.langevin.lang0
6133 allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
6134 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
6135 allocate(fricvec(nres2,nres2)) !(MAXRES2,MAXRES2)
6136 allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
6137 allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
6140 if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2))
6141 !----------------------
6142 ! commom.hairpin in CSA module
6143 !----------------------
6144 ! common.mce in MCM_MD module
6145 !----------------------
6147 ! common /mdgrad/ in module.energy
6148 ! common /back_constr/ in module.energy
6149 ! common /qmeas/ in module.energy
6152 allocate(potEcomp(0:n_ene+4)) !(0:n_ene+4)
6154 allocate(d_t(3,0:nres2),d_a(3,0:nres2),d_t_old(3,0:nres2)) !(3,0:MAXRES2)
6155 allocate(d_a_work(nres6)) !(6*MAXRES)
6157 allocate(DM(nres2),DU1(nres2),DU2(nres2))
6158 allocate(DMorig(nres2),DU1orig(nres2),DU2orig(nres2))
6160 allocate(Gmat(nres2,nres2),A(nres2,nres2))
6161 if(.not.allocated(Ginv)) allocate(Ginv(nres2,nres2)) !in control: ergastulum
6162 allocate(Gsqrp(nres2,nres2),Gsqrm(nres2,nres2),Gvec(nres2,nres2)) !(maxres2,maxres2)
6164 allocate(Geigen(nres2)) !(maxres2)
6165 if(.not.allocated(vtot)) allocate(vtot(nres2)) !(maxres2)
6166 ! common /inertia/ in io_conf: parmread
6167 ! real(kind=8),dimension(:),allocatable :: ISC,msc !(ntyp+1)
6168 ! common /langevin/in io read_MDpar
6169 ! real(kind=8),dimension(:),allocatable :: gamsc !(ntyp1)
6170 ! real(kind=8),dimension(:),allocatable :: stdfsc !(ntyp)
6171 ! in io_conf: parmread
6172 ! real(kind=8),dimension(:),allocatable :: restok !(ntyp+1)
6173 ! common /mdpmpi/ in control: ergastulum
6174 if(.not.allocated(ng_start)) allocate(ng_start(0:nfgtasks-1))
6175 if(.not.allocated(ng_counts)) allocate(ng_counts(0:nfgtasks-1))
6176 if(.not.allocated(nginv_counts)) allocate(nginv_counts(0:nfgtasks-1)) !(0:MaxProcs-1)
6177 if(.not.allocated(nginv_start)) allocate(nginv_start(0:nfgtasks)) !(0:MaxProcs)
6178 !----------------------
6179 ! common.muca in read_muca
6180 ! common /double_muca/
6181 ! real(kind=8) :: elow,ehigh,factor,hbin,factor_min
6182 ! real(kind=8),dimension(:),allocatable :: emuca,nemuca,&
6183 ! nemuca2,hist !(4*maxres)
6184 ! common /integer_muca/
6185 ! integer :: nmuca,imtime,muca_smooth
6187 ! real(kind=8),dimension(:),allocatable :: elowi,ehighi !(maxprocs)
6188 !----------------------
6190 ! common /mdgrad/ in module.energy
6191 ! common /back_constr/ in module.energy
6192 ! common /qmeas/ in module.energy
6196 allocate(d_t_work(nres6),d_t_work_new(nres6),d_af_work(nres6))
6197 allocate(d_as_work(nres6),kinetic_force(nres6)) !(MAXRES6)
6198 allocate(d_t_new(3,0:nres2),d_a_old(3,0:nres2),d_a_short(3,0:nres2)) !,d_a !(3,0:MAXRES2)
6199 allocate(stdforcp(nres),stdforcsc(nres)) !(MAXRES)
6200 !----------------------
6202 allocate(D_ban(nres6)) !(MAXRES6) maxres6=6*maxres
6203 ! common /stochcalc/ stochforcvec
6204 allocate(stochforcvec(nres6)) !(MAXRES6) maxres6=6*maxres
6205 !----------------------
6207 end subroutine alloc_MD_arrays
6208 !-----------------------------------------------------------------------------
6209 !-----------------------------------------------------------------------------