2 !-----------------------------------------------------------------------------
15 !-----------------------------------------------------------------------------
17 ! common /mdgrad/ in module.energy
18 ! common /back_constr/ in module.energy
19 ! common /qmeas/ in module.energy
23 real(kind=8),dimension(:),allocatable :: d_t_work,&
24 d_t_work_new,d_af_work,d_as_work,kinetic_force !(MAXRES6)
25 real(kind=8),dimension(:,:),allocatable :: d_t_new,&
26 d_a_old,d_a_short!,d_a !(3,0:MAXRES2)
27 ! real(kind=8),dimension(:),allocatable :: d_a_work !(6*MAXRES)
28 ! real(kind=8),dimension(:,:),allocatable :: Gmat,Ginv,A,&
29 ! Gsqrp,Gsqrm,Gvec !(maxres2,maxres2)
30 ! real(kind=8),dimension(:),allocatable :: Geigen !(maxres2)
31 ! integer :: dimen,dimen1,dimen3
32 ! integer :: lang,count_reset_moment,count_reset_vel
33 ! logical :: reset_moment,reset_vel,rattle,RESPA
36 ! real(kind=8) :: rwat,etawat,stdfp,cPoise
37 ! real(kind=8),dimension(:),allocatable :: gamsc !(ntyp1)
38 ! real(kind=8),dimension(:),allocatable :: stdfsc !(ntyp)
39 real(kind=8),dimension(:),allocatable :: stdforcp,stdforcsc !(MAXRES)
40 !-----------------------------------------------------------------------------
44 ! ###################################################
45 ! ## COPYRIGHT (C) 1992 by Jay William Ponder ##
46 ! ## All Rights Reserved ##
47 ! ###################################################
49 ! #############################################################
51 ! ## sizes.i -- parameter values to set array dimensions ##
53 ! #############################################################
56 ! "sizes.i" sets values for critical array dimensions used
57 ! throughout the software; these parameters will fix the size
58 ! of the largest systems that can be handled; values too large
59 ! for the computer's memory and/or swap space to accomodate
60 ! will result in poor performance or outright failure
62 ! parameter: maximum allowed number of:
64 ! maxatm atoms in the molecular system
65 ! maxval atoms directly bonded to an atom
66 ! maxgrp ! user-defined groups of atoms
67 ! maxtyp force field atom type definitions
68 ! maxclass force field atom class definitions
69 ! maxkey lines in the keyword file
70 ! maxrot bonds for torsional rotation
71 ! maxvar optimization variables (vector storage)
72 ! maxopt optimization variables (matrix storage)
73 ! maxhess off-diagonal Hessian elements
74 ! maxlight sites for method of lights neighbors
75 ! maxvib vibrational frequencies
76 ! maxgeo distance geometry points
77 ! maxcell unit cells in replicated crystal
78 ! maxring 3-, 4-, or 5-membered rings
79 ! maxfix geometric restraints
80 ! maxbio biopolymer atom definitions
81 ! maxres residues in the macromolecule
82 ! maxamino amino acid residue types
83 ! maxnuc nucleic acid residue types
84 ! maxbnd covalent bonds in molecular system
85 ! maxang bond angles in molecular system
86 ! maxtors torsional angles in molecular system
87 ! maxpi atoms in conjugated pisystem
88 ! maxpib covalent bonds involving pisystem
89 ! maxpit torsional angles involving pisystem
92 !el integer maxatm,maxval,maxgrp
93 !el integer maxtyp,maxclass,maxkey
94 !el integer maxrot,maxopt
95 !el integer maxhess,maxlight,maxvib
96 !el integer maxgeo,maxcell,maxring
97 !el integer maxfix,maxbio
98 !el integer maxamino,maxnuc,maxbnd
99 !el integer maxang,maxtors,maxpi
100 !el integer maxpib,maxpit
101 ! integer :: maxatm !=2*nres !maxres2 maxres2=2*maxres
102 ! integer,parameter :: maxval=8
103 ! integer,parameter :: maxgrp=1000
104 ! integer,parameter :: maxtyp=3000
105 ! integer,parameter :: maxclass=500
106 ! integer,parameter :: maxkey=10000
107 ! integer,parameter :: maxrot=1000
108 ! integer,parameter :: maxopt=1000
109 ! integer,parameter :: maxhess=1000000
110 ! integer :: maxlight !=8*maxatm
111 ! integer,parameter :: maxvib=1000
112 ! integer,parameter :: maxgeo=1000
113 ! integer,parameter :: maxcell=10000
114 ! integer,parameter :: maxring=10000
115 ! integer,parameter :: maxfix=10000
116 ! integer,parameter :: maxbio=10000
117 ! integer,parameter :: maxamino=31
118 ! integer,parameter :: maxnuc=12
119 ! integer :: maxbnd !=2*maxatm
120 ! integer :: maxang !=3*maxatm
121 ! integer :: maxtors !=4*maxatm
122 ! integer,parameter :: maxpi=100
123 ! integer,parameter :: maxpib=2*maxpi
124 ! integer,parameter :: maxpit=4*maxpi
125 !-----------------------------------------------------------------------------
126 ! Maximum number of seed
127 ! integer,parameter :: max_seed=1
128 !-----------------------------------------------------------------------------
129 real(kind=8),dimension(:),allocatable :: stochforcvec !(MAXRES6) maxres6=6*maxres
130 ! common /stochcalc/ stochforcvec
131 !-----------------------------------------------------------------------------
132 ! common /przechowalnia/ subroutines: rattle1,rattle2,rattle_brown
133 real(kind=8),dimension(:,:),allocatable :: GGinv !(2*nres,2*nres) maxres2=2*maxres
134 real(kind=8),dimension(:,:,:),allocatable :: gdc !(3,2*nres,2*nres) maxres2=2*maxres
135 real(kind=8),dimension(:,:),allocatable :: Cmat !(2*nres,2*nres) maxres2=2*maxres
136 !-----------------------------------------------------------------------------
137 ! common /syfek/ subroutines: friction_force,setup_fricmat
138 !el real(kind=8),dimension(:),allocatable :: gamvec !(MAXRES6) or (MAXRES2)
139 !-----------------------------------------------------------------------------
140 ! common /przechowalnia/ subroutines: friction_force,setup_fricmat
141 real(kind=8),dimension(:,:),allocatable :: ginvfric !(2*nres,2*nres) !maxres2=2*maxres
142 !-----------------------------------------------------------------------------
143 ! common /przechowalnia/ subroutine: setup_fricmat
144 !el real(kind=8),dimension(:,:),allocatable :: fcopy !(2*nres,2*nres)
145 !-----------------------------------------------------------------------------
148 !-----------------------------------------------------------------------------
150 !-----------------------------------------------------------------------------
152 !-----------------------------------------------------------------------------
153 subroutine brown_step(itime)
154 !------------------------------------------------
155 ! Perform a single Euler integration step of Brownian dynamics
156 !------------------------------------------------
157 ! implicit real*8 (a-h,o-z)
159 use control, only: tcpu
162 ! use io_conf, only:cartprint
163 ! include 'DIMENSIONS'
167 ! include 'COMMON.CONTROL'
168 ! include 'COMMON.VAR'
169 ! include 'COMMON.MD'
171 ! include 'COMMON.LANGEVIN'
173 ! include 'COMMON.LANGEVIN.lang0'
175 ! include 'COMMON.CHAIN'
176 ! include 'COMMON.DERIV'
177 ! include 'COMMON.GEO'
178 ! include 'COMMON.LOCAL'
179 ! include 'COMMON.INTERACT'
180 ! include 'COMMON.IOUNITS'
181 ! include 'COMMON.NAMES'
182 ! include 'COMMON.TIME1'
183 real(kind=8),dimension(6*nres) :: zapas !(MAXRES6) maxres6=6*maxres
184 integer :: rstcount !ilen,
186 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
187 real(kind=8),dimension(6*nres,2*nres) :: Bmat,GBmat,Tmat !(MAXRES6,MAXRES2) (maxres2=2*maxres,maxres6=6*maxres)
188 real(kind=8),dimension(2*nres,2*nres) :: Cmat_,Cinv !(maxres2,maxres2) maxres2=2*maxres
189 real(kind=8),dimension(6*nres,6*nres) :: Pmat !(maxres6,maxres6) maxres6=6*maxres
190 real(kind=8),dimension(6*nres) :: Td !(maxres6) maxres6=6*maxres
191 real(kind=8),dimension(2*nres) :: ppvec !(maxres2) maxres2=2*maxres
192 !el common /stochcalc/ stochforcvec
193 !el real(kind=8),dimension(3) :: cm !el
194 !el common /gucio/ cm
196 logical :: lprn = .false.,lprn1 = .false.
197 integer :: maxiter = 5
198 real(kind=8) :: difftol = 1.0d-5
199 real(kind=8) :: xx,diffmax,blen2,diffbond,tt0
200 integer :: i,j,nbond,k,ind,ind1,iter
201 integer :: nres2,nres6
206 if (.not.allocated(stochforcvec)) allocate(stochforcvec(nres6)) !(MAXRES6) maxres6=6*maxres
210 if (itype(i,1).ne.10) nbond=nbond+1
214 write (iout,*) "Generalized inverse of fricmat"
215 call matout(dimen,dimen,nres6,nres6,fricmat)
227 Bmat(ind+j,ind1)=dC_norm(j,i)
232 if (itype(i,1).ne.10) then
235 Bmat(ind+j,ind1)=dC_norm(j,i+nres)
241 write (iout,*) "Matrix Bmat"
242 call MATOUT(nbond,dimen,nres6,nres6,Bmat)
248 GBmat(i,j)=GBmat(i,j)+fricmat(i,k)*Bmat(k,j)
253 write (iout,*) "Matrix GBmat"
254 call MATOUT(nbond,dimen,nres6,nres2,Gbmat)
260 Cmat_(i,j)=Cmat_(i,j)+Bmat(k,i)*GBmat(k,j)
265 write (iout,*) "Matrix Cmat"
266 call MATOUT(nbond,nbond,nres2,nres2,Cmat_)
268 call matinvert(nbond,nres2,Cmat_,Cinv,osob)
270 write (iout,*) "Matrix Cinv"
271 call MATOUT(nbond,nbond,nres2,nres2,Cinv)
277 Tmat(i,j)=Tmat(i,j)+GBmat(i,k)*Cinv(k,j)
282 write (iout,*) "Matrix Tmat"
283 call MATOUT(nbond,dimen,nres6,nres2,Tmat)
293 Pmat(i,j)=Pmat(i,j)-Tmat(i,k)*Bmat(j,k)
298 write (iout,*) "Matrix Pmat"
299 call MATOUT(dimen,dimen,nres6,nres6,Pmat)
306 Td(i)=Td(i)+vbl*Tmat(i,ind)
309 if (itype(k,1).ne.10) then
311 Td(i)=Td(i)+vbldsc0(1,itype(k,1))*Tmat(i,ind)
316 write (iout,*) "Vector Td"
318 write (iout,'(i5,f10.5)') i,Td(i)
321 call stochastic_force(stochforcvec)
323 write (iout,*) "stochforcvec"
325 write (iout,*) i,stochforcvec(i)
329 zapas(j)=-gcart(j,0)+stochforcvec(j)
331 dC_work(j)=dC_old(j,0)
337 zapas(ind)=-gcart(j,i)+stochforcvec(ind)
338 dC_work(ind)=dC_old(j,i)
342 if (itype(i,1).ne.10) then
345 zapas(ind)=-gxcart(j,i)+stochforcvec(ind)
346 dC_work(ind)=dC_old(j,i+nres)
352 write (iout,*) "Initial d_t_work"
354 write (iout,*) i,d_t_work(i)
361 d_t_work(i)=d_t_work(i)+fricmat(i,j)*zapas(j)
368 zapas(i)=zapas(i)+Pmat(i,j)*(dC_work(j)+d_t_work(j)*d_time)
372 write (iout,*) "Final d_t_work and zapas"
374 write (iout,*) i,d_t_work(i),zapas(i)
388 dc_work(ind+j)=dc(j,i)
394 d_t(j,i+nres)=d_t_work(ind+j)
395 dc(j,i+nres)=zapas(ind+j)
396 dc_work(ind+j)=dc(j,i+nres)
402 write (iout,*) "Before correction for rotational lengthening"
403 write (iout,*) "New coordinates",&
404 " and differences between actual and standard bond lengths"
409 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
413 if (itype(i,1).ne.10) then
415 xx=vbld(i+nres)-vbldsc0(1,itype(i,1))
416 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
417 i,(dC(j,i+nres),j=1,3),xx
421 ! Second correction (rotational lengthening)
427 blen2 = scalar(dc(1,i),dc(1,i))
428 ppvec(ind)=2*vbl**2-blen2
429 diffbond=dabs(vbl-dsqrt(blen2))
430 if (diffbond.gt.diffmax) diffmax=diffbond
431 if (ppvec(ind).gt.0.0d0) then
432 ppvec(ind)=dsqrt(ppvec(ind))
437 write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind)
441 if (itype(i,1).ne.10) then
443 blen2 = scalar(dc(1,i+nres),dc(1,i+nres))
444 ppvec(ind)=2*vbldsc0(1,itype(i,1))**2-blen2
445 diffbond=dabs(vbldsc0(1,itype(i,1))-dsqrt(blen2))
446 if (diffbond.gt.diffmax) diffmax=diffbond
447 if (ppvec(ind).gt.0.0d0) then
448 ppvec(ind)=dsqrt(ppvec(ind))
453 write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind)
457 if (lprn) write (iout,*) "iter",iter," diffmax",diffmax
458 if (diffmax.lt.difftol) goto 10
462 Td(i)=Td(i)+ppvec(j)*Tmat(i,j)
468 zapas(i)=zapas(i)+Pmat(i,j)*dc_work(j)
479 dc_work(ind+j)=zapas(ind+j)
484 if (itype(i,1).ne.10) then
486 dc(j,i+nres)=zapas(ind+j)
487 dc_work(ind+j)=zapas(ind+j)
492 ! Building the chain from the newly calculated coordinates
495 if (large.and. mod(itime,ntwe).eq.0) then
496 write (iout,*) "Cartesian and internal coordinates: step 1"
499 write (iout,'(a)') "Potential forces"
501 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(-gcart(j,i),j=1,3),&
504 write (iout,'(a)') "Stochastic forces"
506 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(stochforc(j,i),j=1,3),&
507 (stochforc(j,i+nres),j=1,3)
509 write (iout,'(a)') "Velocities"
511 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
512 (d_t(j,i+nres),j=1,3)
517 write (iout,*) "After correction for rotational lengthening"
518 write (iout,*) "New coordinates",&
519 " and differences between actual and standard bond lengths"
524 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
528 if (itype(i,1).ne.10) then
530 xx=vbld(i+nres)-vbldsc0(1,itype(i,1))
531 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
532 i,(dC(j,i+nres),j=1,3),xx
537 ! write (iout,*) "Too many attempts at correcting the bonds"
545 ! Calculate energy and forces
547 call etotal(potEcomp)
548 potE=potEcomp(0)-potEcomp(20)
552 ! Calculate the kinetic and total energy and the kinetic temperature
555 t_enegrad=t_enegrad+MPI_Wtime()-tt0
557 t_enegrad=t_enegrad+tcpu()-tt0
560 kinetic_T=2.0d0/(dimen*Rb)*EK
562 end subroutine brown_step
563 !-----------------------------------------------------------------------------
565 !-----------------------------------------------------------------------------
566 subroutine gauss(RO,AP,MT,M,N,*)
568 ! CALCULATES (RO**(-1))*AP BY GAUSS ELIMINATION
569 ! RO IS A SQUARE MATRIX
570 ! THE CALCULATED PRODUCT IS STORED IN AP
571 ! ABNORMAL EXIT IF RO IS SINGULAR
573 integer :: MT, M, N, M1,I,J,IM,&
575 real(kind=8) :: RO(MT,M),AP(MT,N),X,RM,PR,Y
581 if(dabs(X).le.1.0D-13) return 1
593 if(DABS(RO(J,I)).LE.RM) goto 2
607 if(dabs(X).le.1.0E-13) return 1
616 8 AP(J,K)=AP(J,K)-Y*AP(I,K)
618 9 RO(J,K)=RO(J,K)-Y*RO(I,K)
622 if(dabs(X).le.1.0E-13) return 1
632 15 X=X-AP(K,J)*RO(MI,K)
637 !-----------------------------------------------------------------------------
639 !-----------------------------------------------------------------------------
640 subroutine kinetic(KE_total)
641 !----------------------------------------------------------------
642 ! This subroutine calculates the total kinetic energy of the chain
643 !-----------------------------------------------------------------
645 ! implicit real*8 (a-h,o-z)
646 ! include 'DIMENSIONS'
647 ! include 'COMMON.VAR'
648 ! include 'COMMON.CHAIN'
649 ! include 'COMMON.DERIV'
650 ! include 'COMMON.GEO'
651 ! include 'COMMON.LOCAL'
652 ! include 'COMMON.INTERACT'
653 ! include 'COMMON.MD'
654 ! include 'COMMON.IOUNITS'
655 real(kind=8) :: KE_total
657 integer :: i,j,k,iti,mnum
658 real(kind=8) :: KEt_p,KEt_sc,KEr_p,KEr_sc,incr(3),&
661 write (iout,*) "Velocities, kietic"
663 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
664 (d_t(j,i+nres),j=1,3)
669 ! write (iout,*) "ISC",(isc(itype(i,1)),i=1,nres)
670 ! The translational part for peptide virtual bonds
676 ! write (iout,*) "Kinetic trp:",i,(incr(j),j=1,3)
678 v(j)=incr(j)+0.5d0*d_t(j,i)
680 vtot(i)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
681 KEt_p=KEt_p+mp(mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
683 incr(j)=incr(j)+d_t(j,i)
686 ! write(iout,*) 'KEt_p', KEt_p
687 ! The translational part for the side chain virtual bond
688 ! Only now we can initialize incr with zeros. It must be equal
689 ! to the velocities of the first Calpha.
695 iti=iabs(itype(i,mnum))
696 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
697 if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
698 .or.(mnum.eq.5)) then
704 v(j)=incr(j)+d_t(j,nres+i)
707 ! write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3)
708 ! write (iout,*) "i",i," msc",msc(iti)," v",(v(j),j=1,3)
709 KEt_sc=KEt_sc+msc(iti,mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
710 vtot(i+nres)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
712 incr(j)=incr(j)+d_t(j,i)
716 ! write(iout,*) 'KEt_sc', KEt_sc
717 ! The part due to stretching and rotation of the peptide groups
721 ! write (iout,*) "i",i
722 ! write (iout,*) "i",i," mag1",mag1," mag2",mag2
726 ! write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3)
727 KEr_p=KEr_p+Ip(mnum)*(incr(1)*incr(1)+incr(2)*incr(2) &
731 ! write(iout,*) 'KEr_p', KEr_p
732 ! The rotational part of the side chain virtual bond
736 iti=iabs(itype(i,mnum))
737 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
738 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
739 .and.(mnum.ne.5)) then
741 incr(j)=d_t(j,nres+i)
743 ! write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3)
744 KEr_sc=KEr_sc+Isc(iti,mnum)*(incr(1)*incr(1)+incr(2)*incr(2)+ &
748 ! The total kinetic energy
750 ! write(iout,*) 'KEr_sc', KEr_sc
751 KE_total=0.5d0*(KEt_p+KEt_sc+0.25d0*KEr_p+KEr_sc)
752 ! write (iout,*) "KE_total",KE_total
754 end subroutine kinetic
755 !-----------------------------------------------------------------------------
757 !-----------------------------------------------------------------------------
759 !------------------------------------------------
760 ! The driver for molecular dynamics subroutines
761 !------------------------------------------------
764 use control, only:tcpu,ovrtim
765 ! use io_comm, only:ilen
767 use compare, only:secondary2,hairpin
768 use io, only:cartout,statout
769 ! implicit real*8 (a-h,o-z)
770 ! include 'DIMENSIONS'
773 integer :: IERROR,ERRCODE
775 ! include 'COMMON.SETUP'
776 ! include 'COMMON.CONTROL'
777 ! include 'COMMON.VAR'
778 ! include 'COMMON.MD'
780 ! include 'COMMON.LANGEVIN'
782 ! include 'COMMON.LANGEVIN.lang0'
784 ! include 'COMMON.CHAIN'
785 ! include 'COMMON.DERIV'
786 ! include 'COMMON.GEO'
787 ! include 'COMMON.LOCAL'
788 ! include 'COMMON.INTERACT'
789 ! include 'COMMON.IOUNITS'
790 ! include 'COMMON.NAMES'
791 ! include 'COMMON.TIME1'
792 ! include 'COMMON.HAIRPIN'
793 real(kind=8),dimension(3) :: L,vcm
795 real(kind=8),dimension(6*nres) :: v_work,v_transf !(maxres6) maxres6=6*maxres
797 integer :: rstcount !ilen,
799 character(len=50) :: tytul
800 !el common /gucio/ cm
801 integer :: itime,i,j,nharp
802 integer,dimension(4,nres/3) :: iharp !(4,nres/3)(4,maxres/3)
804 real(kind=8) :: tt0,scalfac
809 if (ilen(tmpdir).gt.0) &
810 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_" &
811 //liczba(:ilen(liczba))//'.rst')
813 if (ilen(tmpdir).gt.0) &
814 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"//'.rst')
821 write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
827 ! Determine the inverse of the inertia matrix.
828 call setup_MD_matrices
832 t_MDsetup = MPI_Wtime()-tt0
834 t_MDsetup = tcpu()-tt0
837 ! Entering the MD loop
843 if (lang.eq.2 .or. lang.eq.3) then
847 call sd_verlet_p_setup
849 call sd_verlet_ciccotti_setup
853 pfric0_mat(i,j,0)=pfric_mat(i,j)
854 afric0_mat(i,j,0)=afric_mat(i,j)
855 vfric0_mat(i,j,0)=vfric_mat(i,j)
856 prand0_mat(i,j,0)=prand_mat(i,j)
857 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
858 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
863 flag_stoch(i)=.false.
867 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
869 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
873 else if (lang.eq.1 .or. lang.eq.4) then
877 t_langsetup=MPI_Wtime()-tt0
880 t_langsetup=tcpu()-tt0
883 do itime=1,n_timestep
885 if (large.and. mod(itime,ntwe).eq.0) &
886 write (iout,*) "itime",itime
888 if (lang.gt.0 .and. surfarea .and. &
889 mod(itime,reset_fricmat).eq.0) then
890 if (lang.eq.2 .or. lang.eq.3) then
894 call sd_verlet_p_setup
896 call sd_verlet_ciccotti_setup
900 pfric0_mat(i,j,0)=pfric_mat(i,j)
901 afric0_mat(i,j,0)=afric_mat(i,j)
902 vfric0_mat(i,j,0)=vfric_mat(i,j)
903 prand0_mat(i,j,0)=prand_mat(i,j)
904 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
905 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
910 flag_stoch(i)=.false.
913 else if (lang.eq.1 .or. lang.eq.4) then
916 write (iout,'(a,i10)') &
917 "Friction matrix reset based on surface area, itime",itime
919 if (reset_vel .and. tbf .and. lang.eq.0 &
920 .and. mod(itime,count_reset_vel).eq.0) then
922 write(iout,'(a,f20.2)') &
923 "Velocities reset to random values, time",totT
926 d_t_old(j,i)=d_t(j,i)
930 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
934 d_t(j,0)=d_t(j,0)-vcm(j)
937 kinetic_T=2.0d0/(dimen3*Rb)*EK
938 scalfac=dsqrt(T_bath/kinetic_T)
939 write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
942 d_t_old(j,i)=scalfac*d_t(j,i)
948 ! Time-reversible RESPA algorithm
949 ! (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
950 call RESPA_step(itime)
952 ! Variable time step algorithm.
953 call velverlet_step(itime)
957 call brown_step(itime)
959 print *,"Brown dynamics not here!"
961 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
967 if (mod(itime,ntwe).eq.0) call statout(itime)
980 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
983 v_work(ind)=d_t(j,i+nres)
988 write (66,'(80f10.5)') &
989 ((d_t(j,i),j=1,3),i=0,nres-1),((d_t(j,i+nres),j=1,3),i=1,nres)
993 v_transf(i)=v_transf(i)+gvec(j,i)*v_work(j)
995 v_transf(i)= v_transf(i)*dsqrt(geigen(i))
997 write (67,'(80f10.5)') (v_transf(i),i=1,ind)
1000 if (mod(itime,ntwx).eq.0) then
1001 write (tytul,'("time",f8.2)') totT
1003 call hairpin(.true.,nharp,iharp)
1004 call secondary2(.true.)
1005 call pdbout(potE,tytul,ipdb)
1010 if (rstcount.eq.1000.or.itime.eq.n_timestep) then
1011 open(irest2,file=rest2name,status='unknown')
1012 write(irest2,*) totT,EK,potE,totE,t_bath
1014 ! AL 4/17/17: Now writing d_t(0,:) too
1016 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
1018 ! AL 4/17/17: Now writing d_c(0,:) too
1020 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
1028 t_MD=MPI_Wtime()-tt0
1032 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))') &
1034 'MD calculations setup:',t_MDsetup,&
1035 'Energy & gradient evaluation:',t_enegrad,&
1036 'Stochastic MD setup:',t_langsetup,&
1037 'Stochastic MD step setup:',t_sdsetup,&
1039 write (iout,'(/28(1h=),a25,27(1h=))') &
1040 ' End of MD calculation '
1042 write (iout,*) "time for etotal",t_etotal," elong",t_elong,&
1044 write (iout,*) "time_fric",time_fric," time_stoch",time_stoch,&
1045 " time_fricmatmult",time_fricmatmult," time_fsample ",&
1050 !-----------------------------------------------------------------------------
1051 subroutine velverlet_step(itime)
1052 !-------------------------------------------------------------------------------
1053 ! Perform a single velocity Verlet step; the time step can be rescaled if
1054 ! increments in accelerations exceed the threshold
1055 !-------------------------------------------------------------------------------
1056 ! implicit real*8 (a-h,o-z)
1057 ! include 'DIMENSIONS'
1059 use control, only:tcpu
1063 integer :: ierror,ierrcode
1064 real(kind=8) :: errcode
1066 ! include 'COMMON.SETUP'
1067 ! include 'COMMON.VAR'
1068 ! include 'COMMON.MD'
1070 ! include 'COMMON.LANGEVIN'
1072 ! include 'COMMON.LANGEVIN.lang0'
1074 ! include 'COMMON.CHAIN'
1075 ! include 'COMMON.DERIV'
1076 ! include 'COMMON.GEO'
1077 ! include 'COMMON.LOCAL'
1078 ! include 'COMMON.INTERACT'
1079 ! include 'COMMON.IOUNITS'
1080 ! include 'COMMON.NAMES'
1081 ! include 'COMMON.TIME1'
1082 ! include 'COMMON.MUCA'
1083 real(kind=8),dimension(3) :: vcm,incr
1084 real(kind=8),dimension(3) :: L
1085 integer :: count,rstcount !ilen,
1087 character(len=50) :: tytul
1088 integer :: maxcount_scale = 20
1089 !el common /gucio/ cm
1090 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1091 !el common /stochcalc/ stochforcvec
1092 integer :: itime,icount_scale,itime_scal,i,j,ifac_time
1094 real(kind=8) :: epdrift,tt0,fac_time
1096 if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres)) !(MAXRES6) maxres6=6*maxres
1102 else if (lang.eq.2 .or. lang.eq.3) then
1104 call stochastic_force(stochforcvec)
1107 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
1109 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1116 icount_scale=icount_scale+1
1117 if (icount_scale.gt.maxcount_scale) then
1119 "ERROR: too many attempts at scaling down the time step. ",&
1120 "amax=",amax,"epdrift=",epdrift,&
1121 "damax=",damax,"edriftmax=",edriftmax,&
1125 call MPI_Abort(MPI_COMM_WORLD,IERROR,IERRCODE)
1129 ! First step of the velocity Verlet algorithm
1134 else if (lang.eq.3) then
1136 call sd_verlet1_ciccotti
1138 else if (lang.eq.1) then
1143 ! Build the chain from the newly calculated coordinates
1144 call chainbuild_cart
1145 if (rattle) call rattle1
1147 if (large.and. mod(itime,ntwe).eq.0) then
1148 write (iout,*) "Cartesian and internal coordinates: step 1"
1153 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),&
1154 (dc(j,i+nres),j=1,3)
1156 write (iout,*) "Accelerations"
1158 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1159 (d_a(j,i+nres),j=1,3)
1161 write (iout,*) "Velocities, step 1"
1163 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1164 (d_t(j,i+nres),j=1,3)
1173 ! Calculate energy and forces
1175 call etotal(potEcomp)
1176 ! AL 4/17/17: Reduce the steps if NaNs occurred.
1177 if (potEcomp(0).gt.0.99e20 .or. isnan(potEcomp(0)).gt.0) then
1182 if (large.and. mod(itime,ntwe).eq.0) &
1183 call enerprint(potEcomp)
1186 t_etotal=t_etotal+MPI_Wtime()-tt0
1188 t_etotal=t_etotal+tcpu()-tt0
1191 potE=potEcomp(0)-potEcomp(20)
1193 ! Get the new accelerations
1196 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1198 t_enegrad=t_enegrad+tcpu()-tt0
1200 ! Determine maximum acceleration and scale down the timestep if needed
1202 amax=amax/(itime_scal+1)**2
1203 call predict_edrift(epdrift)
1204 if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then
1205 ! Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
1207 ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax)) &
1209 itime_scal=itime_scal+ifac_time
1210 ! fac_time=dmin1(damax/amax,0.5d0)
1211 fac_time=0.5d0**ifac_time
1212 d_time=d_time*fac_time
1213 if (lang.eq.2 .or. lang.eq.3) then
1215 ! write (iout,*) "Calling sd_verlet_setup: 1"
1216 ! Rescale the stochastic forces and recalculate or restore
1217 ! the matrices of tinker integrator
1218 if (itime_scal.gt.maxflag_stoch) then
1219 if (large) write (iout,'(a,i5,a)') &
1220 "Calculate matrices for stochastic step;",&
1221 " itime_scal ",itime_scal
1223 call sd_verlet_p_setup
1225 call sd_verlet_ciccotti_setup
1227 write (iout,'(2a,i3,a,i3,1h.)') &
1228 "Warning: cannot store matrices for stochastic",&
1229 " integration because the index",itime_scal,&
1230 " is greater than",maxflag_stoch
1231 write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",&
1232 " integration Langevin algorithm for better efficiency."
1233 else if (flag_stoch(itime_scal)) then
1234 if (large) write (iout,'(a,i5,a,l1)') &
1235 "Restore matrices for stochastic step; itime_scal ",&
1236 itime_scal," flag ",flag_stoch(itime_scal)
1239 pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
1240 afric_mat(i,j)=afric0_mat(i,j,itime_scal)
1241 vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
1242 prand_mat(i,j)=prand0_mat(i,j,itime_scal)
1243 vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
1244 vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
1248 if (large) write (iout,'(2a,i5,a,l1)') &
1249 "Calculate & store matrices for stochastic step;",&
1250 " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
1252 call sd_verlet_p_setup
1254 call sd_verlet_ciccotti_setup
1256 flag_stoch(ifac_time)=.true.
1259 pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
1260 afric0_mat(i,j,itime_scal)=afric_mat(i,j)
1261 vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
1262 prand0_mat(i,j,itime_scal)=prand_mat(i,j)
1263 vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
1264 vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
1268 fac_time=1.0d0/dsqrt(fac_time)
1270 stochforcvec(i)=fac_time*stochforcvec(i)
1273 else if (lang.eq.1) then
1274 ! Rescale the accelerations due to stochastic forces
1275 fac_time=1.0d0/dsqrt(fac_time)
1277 d_as_work(i)=d_as_work(i)*fac_time
1280 if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)') &
1281 "itime",itime," Timestep scaled down to ",&
1282 d_time," ifac_time",ifac_time," itime_scal",itime_scal
1284 ! Second step of the velocity Verlet algorithm
1289 else if (lang.eq.3) then
1291 call sd_verlet2_ciccotti
1293 else if (lang.eq.1) then
1298 if (rattle) call rattle2
1301 if (d_time.ne.d_time0) then
1304 if (lang.eq.2 .or. lang.eq.3) then
1305 if (large) write (iout,'(a)') &
1306 "Restore original matrices for stochastic step"
1307 ! write (iout,*) "Calling sd_verlet_setup: 2"
1308 ! Restore the matrices of tinker integrator if the time step has been restored
1311 pfric_mat(i,j)=pfric0_mat(i,j,0)
1312 afric_mat(i,j)=afric0_mat(i,j,0)
1313 vfric_mat(i,j)=vfric0_mat(i,j,0)
1314 prand_mat(i,j)=prand0_mat(i,j,0)
1315 vrand_mat1(i,j)=vrand0_mat1(i,j,0)
1316 vrand_mat2(i,j)=vrand0_mat2(i,j,0)
1325 ! Calculate the kinetic and the total energy and the kinetic temperature
1329 ! call kinetic1(EK1)
1330 ! write (iout,*) "step",itime," EK",EK," EK1",EK1
1332 ! Couple the system to Berendsen bath if needed
1333 if (tbf .and. lang.eq.0) then
1336 kinetic_T=2.0d0/(dimen3*Rb)*EK
1337 ! Backup the coordinates, velocities, and accelerations
1341 d_t_old(j,i)=d_t(j,i)
1342 d_a_old(j,i)=d_a(j,i)
1346 if (mod(itime,ntwe).eq.0 .and. large) then
1347 write (iout,*) "Velocities, step 2"
1349 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1350 (d_t(j,i+nres),j=1,3)
1355 end subroutine velverlet_step
1356 !-----------------------------------------------------------------------------
1357 subroutine RESPA_step(itime)
1358 !-------------------------------------------------------------------------------
1359 ! Perform a single RESPA step.
1360 !-------------------------------------------------------------------------------
1361 ! implicit real*8 (a-h,o-z)
1362 ! include 'DIMENSIONS'
1366 use control, only:tcpu
1368 ! use io_conf, only:cartprint
1371 integer :: IERROR,ERRCODE
1373 ! include 'COMMON.SETUP'
1374 ! include 'COMMON.CONTROL'
1375 ! include 'COMMON.VAR'
1376 ! include 'COMMON.MD'
1378 ! include 'COMMON.LANGEVIN'
1380 ! include 'COMMON.LANGEVIN.lang0'
1382 ! include 'COMMON.CHAIN'
1383 ! include 'COMMON.DERIV'
1384 ! include 'COMMON.GEO'
1385 ! include 'COMMON.LOCAL'
1386 ! include 'COMMON.INTERACT'
1387 ! include 'COMMON.IOUNITS'
1388 ! include 'COMMON.NAMES'
1389 ! include 'COMMON.TIME1'
1390 real(kind=8),dimension(0:n_ene) :: energia_short,energia_long
1391 real(kind=8),dimension(3) :: L,vcm,incr
1392 real(kind=8),dimension(3,0:2*nres) :: dc_old0,d_t_old0,d_a_old0 !(3,0:maxres2) maxres2=2*maxres
1393 logical :: PRINT_AMTS_MSG = .false.
1394 integer :: count,rstcount !ilen,
1396 character(len=50) :: tytul
1397 integer :: maxcount_scale = 10
1398 !el common /gucio/ cm
1399 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1400 !el common /stochcalc/ stochforcvec
1401 integer :: itime,itt,i,j,itsplit
1403 !el common /cipiszcze/ itt
1405 real(kind=8) :: epdrift,tt0,epdriftmax
1408 if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres)) !(MAXRES6) maxres6=6*maxres
1412 if (large.and. mod(itime,ntwe).eq.0) then
1413 write (iout,*) "***************** RESPA itime",itime
1414 write (iout,*) "Cartesian and internal coordinates: step 0"
1416 call pdbout(0.0d0,"cipiszcze",iout)
1418 write (iout,*) "Accelerations from long-range forces"
1420 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1421 (d_a(j,i+nres),j=1,3)
1423 write (iout,*) "Velocities, step 0"
1425 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1426 (d_t(j,i+nres),j=1,3)
1431 ! Perform the initial RESPA step (increment velocities)
1432 ! write (iout,*) "*********************** RESPA ini"
1435 if (mod(itime,ntwe).eq.0 .and. large) then
1436 write (iout,*) "Velocities, end"
1438 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1439 (d_t(j,i+nres),j=1,3)
1443 ! Compute the short-range forces
1449 ! 7/2/2009 commented out
1451 ! call etotal_short(energia_short)
1454 ! 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
1457 d_a(j,i)=d_a_short(j,i)
1461 if (large.and. mod(itime,ntwe).eq.0) then
1462 write (iout,*) "energia_short",energia_short(0)
1463 write (iout,*) "Accelerations from short-range forces"
1465 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1466 (d_a(j,i+nres),j=1,3)
1471 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1473 t_enegrad=t_enegrad+tcpu()-tt0
1478 d_t_old(j,i)=d_t(j,i)
1479 d_a_old(j,i)=d_a(j,i)
1482 ! 6/30/08 A-MTS: attempt at increasing the split number
1485 dc_old0(j,i)=dc_old(j,i)
1486 d_t_old0(j,i)=d_t_old(j,i)
1487 d_a_old0(j,i)=d_a_old(j,i)
1490 if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
1491 if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
1498 ! write (iout,*) "itime",itime," ntime_split",ntime_split
1499 ! Split the time step
1500 d_time=d_time0/ntime_split
1501 ! Perform the short-range RESPA steps (velocity Verlet increments of
1502 ! positions and velocities using short-range forces)
1503 ! write (iout,*) "*********************** RESPA split"
1504 do itsplit=1,ntime_split
1507 else if (lang.eq.2 .or. lang.eq.3) then
1509 call stochastic_force(stochforcvec)
1512 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
1514 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1519 ! First step of the velocity Verlet algorithm
1524 else if (lang.eq.3) then
1526 call sd_verlet1_ciccotti
1528 else if (lang.eq.1) then
1533 ! Build the chain from the newly calculated coordinates
1534 call chainbuild_cart
1535 if (rattle) call rattle1
1537 if (large.and. mod(itime,ntwe).eq.0) then
1538 write (iout,*) "***** ITSPLIT",itsplit
1539 write (iout,*) "Cartesian and internal coordinates: step 1"
1540 call pdbout(0.0d0,"cipiszcze",iout)
1543 write (iout,*) "Velocities, step 1"
1545 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1546 (d_t(j,i+nres),j=1,3)
1555 ! Calculate energy and forces
1557 call etotal_short(energia_short)
1558 ! AL 4/17/17: Exit itime_split loop when energy goes infinite
1559 if (energia_short(0).gt.0.99e20 .or. isnan(energia_short(0)) ) then
1560 if (PRINT_AMTS_MSG) &
1561 write (iout,*) "Infinities/NaNs in energia_short",energia_short(0),"; increasing ntime_split to",ntime_split
1562 ntime_split=ntime_split*2
1563 if (ntime_split.gt.maxtime_split) then
1566 "Cannot rescue the run; aborting job. Retry with a smaller time step"
1568 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1571 "Cannot rescue the run; terminating. Retry with a smaller time step"
1577 if (large.and. mod(itime,ntwe).eq.0) &
1578 call enerprint(energia_short)
1581 t_eshort=t_eshort+MPI_Wtime()-tt0
1583 t_eshort=t_eshort+tcpu()-tt0
1587 ! Get the new accelerations
1589 ! 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
1592 d_a_short(j,i)=d_a(j,i)
1596 if (large.and. mod(itime,ntwe).eq.0) then
1597 write (iout,*)"energia_short",energia_short(0)
1598 write (iout,*) "Accelerations from short-range forces"
1600 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1601 (d_a(j,i+nres),j=1,3)
1606 ! Determine maximum acceleration and scale down the timestep if needed
1608 amax=amax/ntime_split**2
1609 call predict_edrift(epdrift)
1610 if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0) &
1611 write (iout,*) "amax",amax," damax",damax,&
1612 " epdrift",epdrift," epdriftmax",epdriftmax
1613 ! Exit loop and try with increased split number if the change of
1614 ! acceleration is too big
1615 if (amax.gt.damax .or. epdrift.gt.edriftmax) then
1616 if (ntime_split.lt.maxtime_split) then
1618 ntime_split=ntime_split*2
1619 ! AL 4/17/17: We should exit the itime_split loop when acceleration change is too big
1623 dc_old(j,i)=dc_old0(j,i)
1624 d_t_old(j,i)=d_t_old0(j,i)
1625 d_a_old(j,i)=d_a_old0(j,i)
1628 if (PRINT_AMTS_MSG) then
1629 write (iout,*) "acceleration/energy drift too large",amax,&
1630 epdrift," split increased to ",ntime_split," itime",itime,&
1636 "Uh-hu. Bumpy landscape. Maximum splitting number",&
1638 " already reached!!! Trying to carry on!"
1642 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1644 t_enegrad=t_enegrad+tcpu()-tt0
1646 ! Second step of the velocity Verlet algorithm
1651 else if (lang.eq.3) then
1653 call sd_verlet2_ciccotti
1655 else if (lang.eq.1) then
1660 if (rattle) call rattle2
1661 ! Backup the coordinates, velocities, and accelerations
1665 d_t_old(j,i)=d_t(j,i)
1666 d_a_old(j,i)=d_a(j,i)
1673 ! Restore the time step
1675 ! Compute long-range forces
1682 call etotal_long(energia_long)
1683 if (energia_long(0).gt.0.99e20 .or. isnan(energia_long(0))) then
1686 "Infinitied/NaNs in energia_long, Aborting MPI job."
1688 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
1690 write (iout,*) "Infinitied/NaNs in energia_long, terminating."
1694 if (large.and. mod(itime,ntwe).eq.0) &
1695 call enerprint(energia_long)
1698 t_elong=t_elong+MPI_Wtime()-tt0
1700 t_elong=t_elong+tcpu()-tt0
1706 t_enegrad=t_enegrad+MPI_Wtime()-tt0
1708 t_enegrad=t_enegrad+tcpu()-tt0
1710 ! Compute accelerations from long-range forces
1712 if (large.and. mod(itime,ntwe).eq.0) then
1713 write (iout,*) "energia_long",energia_long(0)
1714 write (iout,*) "Cartesian and internal coordinates: step 2"
1716 call pdbout(0.0d0,"cipiszcze",iout)
1718 write (iout,*) "Accelerations from long-range forces"
1720 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
1721 (d_a(j,i+nres),j=1,3)
1723 write (iout,*) "Velocities, step 2"
1725 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1726 (d_t(j,i+nres),j=1,3)
1730 ! Compute the final RESPA step (increment velocities)
1731 ! write (iout,*) "*********************** RESPA fin"
1733 ! Compute the complete potential energy
1735 potEcomp(i)=energia_short(i)+energia_long(i)
1737 potE=potEcomp(0)-potEcomp(20)
1738 ! potE=energia_short(0)+energia_long(0)
1741 ! Calculate the kinetic and the total energy and the kinetic temperature
1744 ! Couple the system to Berendsen bath if needed
1745 if (tbf .and. lang.eq.0) then
1748 kinetic_T=2.0d0/(dimen3*Rb)*EK
1749 ! Backup the coordinates, velocities, and accelerations
1751 if (mod(itime,ntwe).eq.0 .and. large) then
1752 write (iout,*) "Velocities, end"
1754 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
1755 (d_t(j,i+nres),j=1,3)
1760 end subroutine RESPA_step
1761 !-----------------------------------------------------------------------------
1762 subroutine RESPA_vel
1763 ! First and last RESPA step (incrementing velocities using long-range
1766 ! implicit real*8 (a-h,o-z)
1767 ! include 'DIMENSIONS'
1768 ! include 'COMMON.CONTROL'
1769 ! include 'COMMON.VAR'
1770 ! include 'COMMON.MD'
1771 ! include 'COMMON.CHAIN'
1772 ! include 'COMMON.DERIV'
1773 ! include 'COMMON.GEO'
1774 ! include 'COMMON.LOCAL'
1775 ! include 'COMMON.INTERACT'
1776 ! include 'COMMON.IOUNITS'
1777 ! include 'COMMON.NAMES'
1778 integer :: i,j,inres,mnum
1781 d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
1785 d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
1790 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
1791 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
1792 .and.(mnum.ne.5)) then
1795 d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
1800 end subroutine RESPA_vel
1801 !-----------------------------------------------------------------------------
1803 ! Applying velocity Verlet algorithm - step 1 to coordinates
1805 ! implicit real*8 (a-h,o-z)
1806 ! include 'DIMENSIONS'
1807 ! include 'COMMON.CONTROL'
1808 ! include 'COMMON.VAR'
1809 ! include 'COMMON.MD'
1810 ! include 'COMMON.CHAIN'
1811 ! include 'COMMON.DERIV'
1812 ! include 'COMMON.GEO'
1813 ! include 'COMMON.LOCAL'
1814 ! include 'COMMON.INTERACT'
1815 ! include 'COMMON.IOUNITS'
1816 ! include 'COMMON.NAMES'
1817 real(kind=8) :: adt,adt2
1818 integer :: i,j,inres,mnum
1821 write (iout,*) "VELVERLET1 START: DC"
1823 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),&
1824 (dc(j,i+nres),j=1,3)
1828 adt=d_a_old(j,0)*d_time
1830 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1831 d_t_new(j,0)=d_t_old(j,0)+adt2
1832 d_t(j,0)=d_t_old(j,0)+adt
1836 adt=d_a_old(j,i)*d_time
1838 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1839 d_t_new(j,i)=d_t_old(j,i)+adt2
1840 d_t(j,i)=d_t_old(j,i)+adt
1845 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
1846 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
1847 .and.(mnum.ne.5)) then
1850 adt=d_a_old(j,inres)*d_time
1852 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
1853 d_t_new(j,inres)=d_t_old(j,inres)+adt2
1854 d_t(j,inres)=d_t_old(j,inres)+adt
1859 write (iout,*) "VELVERLET1 END: DC"
1861 write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),&
1862 (dc(j,i+nres),j=1,3)
1866 end subroutine verlet1
1867 !-----------------------------------------------------------------------------
1869 ! Step 2 of the velocity Verlet algorithm: update velocities
1871 ! implicit real*8 (a-h,o-z)
1872 ! include 'DIMENSIONS'
1873 ! include 'COMMON.CONTROL'
1874 ! include 'COMMON.VAR'
1875 ! include 'COMMON.MD'
1876 ! include 'COMMON.CHAIN'
1877 ! include 'COMMON.DERIV'
1878 ! include 'COMMON.GEO'
1879 ! include 'COMMON.LOCAL'
1880 ! include 'COMMON.INTERACT'
1881 ! include 'COMMON.IOUNITS'
1882 ! include 'COMMON.NAMES'
1883 integer :: i,j,inres,mnum
1886 d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
1890 d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
1895 ! iti=iabs(itype(i,mnum))
1896 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
1897 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
1898 .and.(mnum.ne.5)) then
1901 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
1906 end subroutine verlet2
1907 !-----------------------------------------------------------------------------
1908 subroutine sddir_precalc
1909 ! Applying velocity Verlet algorithm - step 1 to coordinates
1910 ! implicit real*8 (a-h,o-z)
1911 ! include 'DIMENSIONS'
1917 ! include 'COMMON.CONTROL'
1918 ! include 'COMMON.VAR'
1919 ! include 'COMMON.MD'
1921 ! include 'COMMON.LANGEVIN'
1923 ! include 'COMMON.LANGEVIN.lang0'
1925 ! include 'COMMON.CHAIN'
1926 ! include 'COMMON.DERIV'
1927 ! include 'COMMON.GEO'
1928 ! include 'COMMON.LOCAL'
1929 ! include 'COMMON.INTERACT'
1930 ! include 'COMMON.IOUNITS'
1931 ! include 'COMMON.NAMES'
1932 ! include 'COMMON.TIME1'
1933 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
1934 !el common /stochcalc/ stochforcvec
1935 real(kind=8) :: time00
1937 ! Compute friction and stochastic forces
1942 time_fric=time_fric+MPI_Wtime()-time00
1944 call stochastic_force(stochforcvec)
1945 time_stoch=time_stoch+MPI_Wtime()-time00
1948 ! Compute the acceleration due to friction forces (d_af_work) and stochastic
1949 ! forces (d_as_work)
1951 call ginv_mult(fric_work, d_af_work)
1952 call ginv_mult(stochforcvec, d_as_work)
1954 end subroutine sddir_precalc
1955 !-----------------------------------------------------------------------------
1956 subroutine sddir_verlet1
1957 ! Applying velocity Verlet algorithm - step 1 to velocities
1960 ! implicit real*8 (a-h,o-z)
1961 ! include 'DIMENSIONS'
1962 ! include 'COMMON.CONTROL'
1963 ! include 'COMMON.VAR'
1964 ! include 'COMMON.MD'
1966 ! include 'COMMON.LANGEVIN'
1968 ! include 'COMMON.LANGEVIN.lang0'
1970 ! include 'COMMON.CHAIN'
1971 ! include 'COMMON.DERIV'
1972 ! include 'COMMON.GEO'
1973 ! include 'COMMON.LOCAL'
1974 ! include 'COMMON.INTERACT'
1975 ! include 'COMMON.IOUNITS'
1976 ! include 'COMMON.NAMES'
1977 ! Revised 3/31/05 AL: correlation between random contributions to
1978 ! position and velocity increments included.
1979 real(kind=8) :: sqrt13 = 0.57735026918962576451d0 ! 1/sqrt(3)
1980 real(kind=8) :: adt,adt2
1981 integer :: i,j,ind,inres,mnum
1983 ! Add the contribution from BOTH friction and stochastic force to the
1984 ! coordinates, but ONLY the contribution from the friction forces to velocities
1987 adt=(d_a_old(j,0)+d_af_work(j))*d_time
1988 adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
1989 dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
1990 d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
1991 d_t(j,0)=d_t_old(j,0)+adt
1996 adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
1997 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
1998 dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
1999 d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
2000 d_t(j,i)=d_t_old(j,i)+adt
2006 ! iti=iabs(itype(i,mnum))
2007 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2008 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2009 .and.(mnum.ne.5)) then
2012 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
2013 adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
2014 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
2015 d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
2016 d_t(j,inres)=d_t_old(j,inres)+adt
2022 end subroutine sddir_verlet1
2023 !-----------------------------------------------------------------------------
2024 subroutine sddir_verlet2
2025 ! Calculating the adjusted velocities for accelerations
2028 ! implicit real*8 (a-h,o-z)
2029 ! include 'DIMENSIONS'
2030 ! include 'COMMON.CONTROL'
2031 ! include 'COMMON.VAR'
2032 ! include 'COMMON.MD'
2034 ! include 'COMMON.LANGEVIN'
2036 ! include 'COMMON.LANGEVIN.lang0'
2038 ! include 'COMMON.CHAIN'
2039 ! include 'COMMON.DERIV'
2040 ! include 'COMMON.GEO'
2041 ! include 'COMMON.LOCAL'
2042 ! include 'COMMON.INTERACT'
2043 ! include 'COMMON.IOUNITS'
2044 ! include 'COMMON.NAMES'
2045 real(kind=8),dimension(6*nres) :: stochforcvec,d_as_work1 !(MAXRES6) maxres6=6*maxres
2046 real(kind=8) :: cos60 = 0.5d0, sin60 = 0.86602540378443864676d0
2047 integer :: i,j,ind,inres,mnum
2048 ! Revised 3/31/05 AL: correlation between random contributions to
2049 ! position and velocity increments included.
2050 ! The correlation coefficients are calculated at low-friction limit.
2051 ! Also, friction forces are now not calculated with new velocities.
2053 ! call friction_force
2054 call stochastic_force(stochforcvec)
2056 ! Compute the acceleration due to friction forces (d_af_work) and stochastic
2057 ! forces (d_as_work)
2059 call ginv_mult(stochforcvec, d_as_work1)
2065 d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j)) &
2066 +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
2071 d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j)) &
2072 +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
2078 ! iti=iabs(itype(i,mnum))
2079 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2080 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2081 .and.(mnum.ne.5)) then
2084 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres) &
2085 +d_af_work(ind+j))+sin60*d_as_work(ind+j) &
2086 +cos60*d_as_work1(ind+j))*d_time
2092 end subroutine sddir_verlet2
2093 !-----------------------------------------------------------------------------
2094 subroutine max_accel
2096 ! Find the maximum difference in the accelerations of the the sites
2097 ! at the beginning and the end of the time step.
2100 ! implicit real*8 (a-h,o-z)
2101 ! include 'DIMENSIONS'
2102 ! include 'COMMON.CONTROL'
2103 ! include 'COMMON.VAR'
2104 ! include 'COMMON.MD'
2105 ! include 'COMMON.CHAIN'
2106 ! include 'COMMON.DERIV'
2107 ! include 'COMMON.GEO'
2108 ! include 'COMMON.LOCAL'
2109 ! include 'COMMON.INTERACT'
2110 ! include 'COMMON.IOUNITS'
2111 real(kind=8),dimension(3) :: aux,accel,accel_old
2112 real(kind=8) :: dacc
2116 ! aux(j)=d_a(j,0)-d_a_old(j,0)
2117 accel_old(j)=d_a_old(j,0)
2124 ! 7/3/08 changed to asymmetric difference
2126 ! accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
2127 accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
2128 accel(j)=accel(j)+0.5d0*d_a(j,i)
2129 ! if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
2130 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
2131 dacc=dabs(accel(j)-accel_old(j))
2132 ! write (iout,*) i,dacc
2133 if (dacc.gt.amax) amax=dacc
2141 accel_old(j)=d_a_old(j,0)
2146 accel_old(j)=accel_old(j)+d_a_old(j,1)
2147 accel(j)=accel(j)+d_a(j,1)
2152 ! iti=iabs(itype(i,mnum))
2153 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2154 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2155 .and.(mnum.ne.5)) then
2157 ! accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
2158 accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
2159 accel(j)=accel(j)+d_a(j,i+nres)
2163 ! if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
2164 if (dabs(accel(j)).gt.dabs(accel_old(j))) then
2165 dacc=dabs(accel(j)-accel_old(j))
2166 ! write (iout,*) "side-chain",i,dacc
2167 if (dacc.gt.amax) amax=dacc
2171 accel_old(j)=accel_old(j)+d_a_old(j,i)
2172 accel(j)=accel(j)+d_a(j,i)
2173 ! aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
2177 end subroutine max_accel
2178 !-----------------------------------------------------------------------------
2179 subroutine predict_edrift(epdrift)
2181 ! Predict the drift of the potential energy
2184 use control_data, only: lmuca
2185 ! implicit real*8 (a-h,o-z)
2186 ! include 'DIMENSIONS'
2187 ! include 'COMMON.CONTROL'
2188 ! include 'COMMON.VAR'
2189 ! include 'COMMON.MD'
2190 ! include 'COMMON.CHAIN'
2191 ! include 'COMMON.DERIV'
2192 ! include 'COMMON.GEO'
2193 ! include 'COMMON.LOCAL'
2194 ! include 'COMMON.INTERACT'
2195 ! include 'COMMON.IOUNITS'
2196 ! include 'COMMON.MUCA'
2197 real(kind=8) :: epdrift,epdriftij
2199 ! Drift of the potential energy
2205 epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
2206 if (lmuca) epdriftij=epdriftij*factor
2207 ! write (iout,*) "back",i,j,epdriftij
2208 if (epdriftij.gt.epdrift) epdrift=epdriftij
2212 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2215 dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
2216 if (lmuca) epdriftij=epdriftij*factor
2217 ! write (iout,*) "side",i,j,epdriftij
2218 if (epdriftij.gt.epdrift) epdrift=epdriftij
2222 epdrift=0.5d0*epdrift*d_time*d_time
2223 ! write (iout,*) "epdrift",epdrift
2225 end subroutine predict_edrift
2226 !-----------------------------------------------------------------------------
2227 subroutine verlet_bath
2229 ! Coupling to the thermostat by using the Berendsen algorithm
2232 ! implicit real*8 (a-h,o-z)
2233 ! include 'DIMENSIONS'
2234 ! include 'COMMON.CONTROL'
2235 ! include 'COMMON.VAR'
2236 ! include 'COMMON.MD'
2237 ! include 'COMMON.CHAIN'
2238 ! include 'COMMON.DERIV'
2239 ! include 'COMMON.GEO'
2240 ! include 'COMMON.LOCAL'
2241 ! include 'COMMON.INTERACT'
2242 ! include 'COMMON.IOUNITS'
2243 ! include 'COMMON.NAMES'
2244 real(kind=8) :: T_half,fact
2245 integer :: i,j,inres,mnum
2247 T_half=2.0d0/(dimen3*Rb)*EK
2248 fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
2249 ! write(iout,*) "T_half", T_half
2250 ! write(iout,*) "EK", EK
2251 ! write(iout,*) "fact", fact
2253 d_t(j,0)=fact*d_t(j,0)
2257 d_t(j,i)=fact*d_t(j,i)
2262 ! iti=iabs(itype(i,mnum))
2263 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2264 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2265 .and.(mnum.ne.5)) then
2268 d_t(j,inres)=fact*d_t(j,inres)
2273 end subroutine verlet_bath
2274 !-----------------------------------------------------------------------------
2276 ! Set up the initial conditions of a MD simulation
2279 use control, only:tcpu
2280 !el use io_basic, only:ilen
2283 use minimm, only:minim_dc,minimize,sc_move
2284 use io_config, only:readrst
2285 use io, only:statout
2286 ! implicit real*8 (a-h,o-z)
2287 ! include 'DIMENSIONS'
2290 character(len=16) :: form
2291 integer :: IERROR,ERRCODE
2293 ! include 'COMMON.SETUP'
2294 ! include 'COMMON.CONTROL'
2295 ! include 'COMMON.VAR'
2296 ! include 'COMMON.MD'
2298 ! include 'COMMON.LANGEVIN'
2300 ! include 'COMMON.LANGEVIN.lang0'
2302 ! include 'COMMON.CHAIN'
2303 ! include 'COMMON.DERIV'
2304 ! include 'COMMON.GEO'
2305 ! include 'COMMON.LOCAL'
2306 ! include 'COMMON.INTERACT'
2307 ! include 'COMMON.IOUNITS'
2308 ! include 'COMMON.NAMES'
2309 ! include 'COMMON.REMD'
2310 real(kind=8),dimension(0:n_ene) :: energia_long,energia_short
2311 real(kind=8),dimension(3) :: vcm,incr,L
2312 real(kind=8) :: xv,sigv,lowb,highb
2313 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
2314 character(len=256) :: qstr
2317 character(len=50) :: tytul
2318 logical :: file_exist
2319 !el common /gucio/ cm
2320 integer :: i,j,ipos,iq,iw,nft_sc,iretcode,nfun,itime,ierr,mnum
2321 real(kind=8) :: etot,tt0
2325 ! write(iout,*) "d_time", d_time
2326 ! Compute the standard deviations of stochastic forces for Langevin dynamics
2327 ! if the friction coefficients do not depend on surface area
2328 if (lang.gt.0 .and. .not.surfarea) then
2331 stdforcp(i)=stdfp(mnum)*dsqrt(gamp(mnum))
2335 stdforcsc(i)=stdfsc(iabs(itype(i,mnum)),mnum) &
2336 *dsqrt(gamsc(iabs(itype(i,mnum)),mnum))
2339 ! Open the pdb file for snapshotshots
2342 if (ilen(tmpdir).gt.0) &
2343 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2344 liczba(:ilen(liczba))//".pdb")
2346 file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2350 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) &
2351 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2352 liczba(:ilen(liczba))//".x")
2353 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2356 if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) &
2357 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// &
2358 liczba(:ilen(liczba))//".cx")
2359 cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) &
2365 if (ilen(tmpdir).gt.0) &
2366 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
2367 open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
2369 if (ilen(tmpdir).gt.0) &
2370 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
2371 cartname=prefix(:ilen(prefix))//"_MD.cx"
2375 write (qstr,'(256(1h ))')
2378 iq = qinfrag(i,iset)*10
2379 iw = wfrag(i,iset)/100
2381 if(me.eq.king.or..not.out1file) &
2382 write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
2383 write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
2388 iq = qinpair(i,iset)*10
2389 iw = wpair(i,iset)/100
2391 if(me.eq.king.or..not.out1file) &
2392 write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
2393 write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
2397 ! pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
2399 ! cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
2401 ! cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
2403 ! statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
2407 if (restart1file) then
2409 inquire(file=mremd_rst_name,exist=file_exist)
2410 write (*,*) me," Before broadcast: file_exist",file_exist
2412 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,&
2415 write (*,*) me," After broadcast: file_exist",file_exist
2416 ! inquire(file=mremd_rst_name,exist=file_exist)
2417 if(me.eq.king.or..not.out1file) &
2418 write(iout,*) "Initial state read by master and distributed"
2420 if (ilen(tmpdir).gt.0) &
2421 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_' &
2422 //liczba(:ilen(liczba))//'.rst')
2423 inquire(file=rest2name,exist=file_exist)
2426 if(.not.restart1file) then
2427 if(me.eq.king.or..not.out1file) &
2428 write(iout,*) "Initial state will be read from file ",&
2429 rest2name(:ilen(rest2name))
2432 call rescale_weights(t_bath)
2434 if(me.eq.king.or..not.out1file)then
2435 if (restart1file) then
2436 write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),&
2439 write(iout,*) "File ",rest2name(:ilen(rest2name)),&
2442 write(iout,*) "Initial velocities randomly generated"
2449 ! Generate initial velocities
2450 if(me.eq.king.or..not.out1file) &
2451 write(iout,*) "Initial velocities randomly generated"
2456 ! rest2name = prefix(:ilen(prefix))//'.rst'
2457 if(me.eq.king.or..not.out1file)then
2458 write (iout,*) "Initial velocities"
2460 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2461 (d_t(j,i+nres),j=1,3)
2463 ! Zeroing the total angular momentum of the system
2464 write(iout,*) "Calling the zero-angular momentum subroutine"
2467 ! Getting the potential energy and forces and velocities and accelerations
2469 ! write (iout,*) "velocity of the center of the mass:"
2470 ! write (iout,*) (vcm(j),j=1,3)
2472 d_t(j,0)=d_t(j,0)-vcm(j)
2474 ! Removing the velocity of the center of mass
2476 if(me.eq.king.or..not.out1file)then
2477 write (iout,*) "vcm right after adjustment:"
2478 write (iout,*) (vcm(j),j=1,3)
2480 if ((.not.rest).and.(indpdb.eq.0)) then
2482 if(iranconf.ne.0) then
2484 print *, 'Calling OVERLAP_SC'
2485 call overlap_sc(fail)
2488 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
2489 print *,'SC_move',nft_sc,etot
2490 if(me.eq.king.or..not.out1file) &
2491 write(iout,*) 'SC_move',nft_sc,etot
2495 print *, 'Calling MINIM_DC'
2496 call minim_dc(etot,iretcode,nfun)
2498 call geom_to_var(nvar,varia)
2499 print *,'Calling MINIMIZE.'
2500 call minimize(etot,varia,iretcode,nfun)
2501 call var_to_geom(nvar,varia)
2503 if(me.eq.king.or..not.out1file) &
2504 write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
2507 call chainbuild_cart
2512 kinetic_T=2.0d0/(dimen3*Rb)*EK
2513 if(me.eq.king.or..not.out1file)then
2523 call etotal(potEcomp)
2524 if (large) call enerprint(potEcomp)
2527 t_etotal=t_etotal+MPI_Wtime()-tt0
2529 t_etotal=t_etotal+tcpu()-tt0
2536 if (amax*d_time .gt. dvmax) then
2537 d_time=d_time*dvmax/amax
2538 if(me.eq.king.or..not.out1file) write (iout,*) &
2539 "Time step reduced to",d_time,&
2540 " because of too large initial acceleration."
2542 if(me.eq.king.or..not.out1file)then
2543 write(iout,*) "Potential energy and its components"
2544 call enerprint(potEcomp)
2545 ! write(iout,*) (potEcomp(i),i=0,n_ene)
2547 potE=potEcomp(0)-potEcomp(20)
2550 if (ntwe.ne.0) call statout(itime)
2551 if(me.eq.king.or..not.out1file) &
2552 write (iout,'(/a/3(a25,1pe14.5/))') "Initial:", &
2553 " Kinetic energy",EK," Potential energy",potE, &
2554 " Total energy",totE," Maximum acceleration ", &
2557 write (iout,*) "Initial coordinates"
2559 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),&
2562 write (iout,*) "Initial dC"
2564 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),&
2565 (dc(j,i+nres),j=1,3)
2567 write (iout,*) "Initial velocities"
2568 write (iout,"(13x,' backbone ',23x,' side chain')")
2570 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
2571 (d_t(j,i+nres),j=1,3)
2573 write (iout,*) "Initial accelerations"
2575 ! write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
2576 write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),&
2577 (d_a(j,i+nres),j=1,3)
2583 d_t_old(j,i)=d_t(j,i)
2584 d_a_old(j,i)=d_a(j,i)
2586 ! write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
2595 call etotal_short(energia_short)
2596 if (large) call enerprint(potEcomp)
2599 t_eshort=t_eshort+MPI_Wtime()-tt0
2601 t_eshort=t_eshort+tcpu()-tt0
2606 if(.not.out1file .and. large) then
2607 write (iout,*) "energia_long",energia_long(0),&
2608 " energia_short",energia_short(0),&
2609 " total",energia_long(0)+energia_short(0)
2610 write (iout,*) "Initial fast-force accelerations"
2612 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
2613 (d_a(j,i+nres),j=1,3)
2616 ! 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
2619 d_a_short(j,i)=d_a(j,i)
2628 call etotal_long(energia_long)
2629 if (large) call enerprint(potEcomp)
2632 t_elong=t_elong+MPI_Wtime()-tt0
2634 t_elong=t_elong+tcpu()-tt0
2639 if(.not.out1file .and. large) then
2640 write (iout,*) "energia_long",energia_long(0)
2641 write (iout,*) "Initial slow-force accelerations"
2643 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),&
2644 (d_a(j,i+nres),j=1,3)
2648 t_enegrad=t_enegrad+MPI_Wtime()-tt0
2650 t_enegrad=t_enegrad+tcpu()-tt0
2654 end subroutine init_MD
2655 !-----------------------------------------------------------------------------
2656 subroutine random_vel
2658 ! implicit real*8 (a-h,o-z)
2660 use random, only:anorm_distr
2662 ! include 'DIMENSIONS'
2663 ! include 'COMMON.CONTROL'
2664 ! include 'COMMON.VAR'
2665 ! include 'COMMON.MD'
2667 ! include 'COMMON.LANGEVIN'
2669 ! include 'COMMON.LANGEVIN.lang0'
2671 ! include 'COMMON.CHAIN'
2672 ! include 'COMMON.DERIV'
2673 ! include 'COMMON.GEO'
2674 ! include 'COMMON.LOCAL'
2675 ! include 'COMMON.INTERACT'
2676 ! include 'COMMON.IOUNITS'
2677 ! include 'COMMON.NAMES'
2678 ! include 'COMMON.TIME1'
2679 real(kind=8) :: xv,sigv,lowb,highb ,Ek1
2681 real(kind=8) ,allocatable, dimension(:) :: DDU1,DDU2,DL2,DL1,xsolv,DML,rs
2682 real(kind=8) :: sumx
2684 real(kind=8) ,allocatable, dimension(:) :: rsold
2685 real (kind=8),allocatable,dimension(:,:) :: matold
2688 integer :: i,j,ii,k,ind,mark,imark,mnum
2689 ! Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
2690 ! First generate velocities in the eigenspace of the G matrix
2691 ! write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
2698 sigv=dsqrt((Rb*t_bath)/geigen(i))
2701 d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
2703 write (iout,*) "i",i," ii",ii," geigen",geigen(i),&
2704 " d_t_work_new",d_t_work_new(ii)
2715 Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
2718 write (iout,*) "Ek from eigenvectors",Ek1
2719 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
2723 ! Transform velocities to UNRES coordinate space
2724 allocate (DL1(2*nres))
2725 allocate (DDU1(2*nres))
2726 allocate (DL2(2*nres))
2727 allocate (DDU2(2*nres))
2728 allocate (xsolv(2*nres))
2729 allocate (DML(2*nres))
2730 allocate (rs(2*nres))
2732 allocate (rsold(2*nres))
2733 allocate (matold(2*nres,2*nres))
2735 matold(1,1)=DMorig(1)
2736 matold(1,2)=DU1orig(1)
2737 matold(1,3)=DU2orig(1)
2738 write (*,*) DMorig(1),DU1orig(1),DU2orig(1)
2743 matold(i,j)=DMorig(i)
2744 matold(i,j-1)=DU1orig(i-1)
2746 matold(i,j-2)=DU2orig(i-2)
2754 matold(i,j+1)=DU1orig(i)
2760 matold(i,j+2)=DU2orig(i)
2764 matold(dimen,dimen)=DMorig(dimen)
2765 matold(dimen,dimen-1)=DU1orig(dimen-1)
2766 matold(dimen,dimen-2)=DU2orig(dimen-2)
2767 write(iout,*) "old gmatrix"
2768 call matout(dimen,dimen,2*nres,2*nres,matold)
2772 ! Find the ith eigenvector of the pentadiagonal inertiq matrix
2776 DML(j)=DMorig(j)-geigen(i)
2779 DML(j-1)=DMorig(j)-geigen(i)
2784 DDU1(imark-1)=DU2orig(imark-1)
2785 do j=imark+1,dimen-1
2786 DDU1(j-1)=DU1orig(j)
2794 DDU2(j)=DU2orig(j+1)
2803 write (iout,*) "DL2,DL1,DML,DDU1,DDU2"
2804 write(iout,'(10f10.5)') (DL2(k),k=3,dimen-1)
2805 write(iout,'(10f10.5)') (DL1(k),k=2,dimen-1)
2806 write(iout,'(10f10.5)') (DML(k),k=1,dimen-1)
2807 write(iout,'(10f10.5)') (DDU1(k),k=1,dimen-2)
2808 write(iout,'(10f10.5)') (DDU2(k),k=1,dimen-3)
2811 if (imark.gt.2) rs(imark-2)=-DU2orig(imark-2)
2812 if (imark.gt.1) rs(imark-1)=-DU1orig(imark-1)
2813 if (imark.lt.dimen) rs(imark)=-DU1orig(imark)
2814 if (imark.lt.dimen-1) rs(imark+1)=-DU2orig(imark)
2818 ! write (iout,*) "Vector rs"
2820 ! write (iout,*) j,rs(j)
2823 call FDIAG(dimen-1,DL2,DL1,DML,DDU1,DDU2,rs,xsolv,mark)
2830 sumx=-geigen(i)*xsolv(j)
2832 sumx=sumx+matold(j,k)*xsolv(k)
2835 sumx=sumx+matold(j,k)*xsolv(k-1)
2837 write(iout,'(i5,3f10.5)') j,sumx,rsold(j),sumx-rsold(j)
2840 sumx=-geigen(i)*xsolv(j-1)
2842 sumx=sumx+matold(j,k)*xsolv(k)
2845 sumx=sumx+matold(j,k)*xsolv(k-1)
2847 write(iout,'(i5,3f10.5)') j-1,sumx,rsold(j-1),sumx-rsold(j-1)
2851 "Solution of equations system",i," for eigenvalue",geigen(i)
2853 write(iout,'(i5,f10.5)') j,xsolv(j)
2856 do j=dimen-1,imark,-1
2861 write (iout,*) "Un-normalized eigenvector",i," for eigenvalue",geigen(i)
2863 write(iout,'(i5,f10.5)') j,xsolv(j)
2866 ! Normalize ith eigenvector
2869 sumx=sumx+xsolv(j)**2
2873 xsolv(j)=xsolv(j)/sumx
2876 write (iout,*) "Eigenvector",i," for eigenvalue",geigen(i)
2878 write(iout,'(i5,3f10.5)') j,xsolv(j)
2881 ! All done at this point for eigenvector i, exit loop
2889 write (iout,*) "Unable to find eigenvector",i
2892 ! write (iout,*) "i=",i
2894 ! write (iout,*) "k=",k
2897 ! write(iout,*) "ind",ind," ind1",3*(i-1)+k
2898 d_t_work(ind)=d_t_work(ind) &
2899 +xsolv(j)*d_t_work_new(3*(i-1)+k)
2902 enddo ! i (loop over the eigenvectors)
2905 write (iout,*) "d_t_work"
2907 write (iout,"(i5,f10.5)") i,d_t_work(i)
2912 ! if (itype(i,1).eq.10) then
2914 iti=iabs(itype(i,mnum))
2915 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
2916 if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
2917 .or.(mnum.eq.5)) then
2924 ! write (iout,*) "k",k," ii+k",ii+k," ii+j+k",ii+j+k,"EK1",Ek1
2925 Ek1=Ek1+0.5d0*mp(mnum)*((d_t_work(ii+j+k)+d_t_work_new(ii+k))/2)**2+&
2926 0.5d0*0.25d0*IP(mnum)*(d_t_work(ii+j+k)-d_t_work_new(ii+k))**2
2929 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2930 .and.(mnum.ne.5)) ii=ii+3
2931 write (iout,*) "i",i," itype",itype(i,mnum)," mass",msc(itype(i,mnum),mnum)
2932 write (iout,*) "ii",ii
2935 write (iout,*) "k",k," ii",ii,"EK1",EK1
2936 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2938 Ek1=Ek1+0.5d0*Isc(iabs(itype(i,mnum),mnum))*(d_t_work(ii)-d_t_work(ii-3))**2
2939 Ek1=Ek1+0.5d0*msc(iabs(itype(i,mnum)),mnum)*d_t_work(ii)**2
2941 write (iout,*) "i",i," ii",ii
2943 write (iout,*) "Ek from d_t_work",Ek1
2944 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
2960 d_t(k,j)=d_t_work(ind)
2964 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
2967 d_t(k,j+nres)=d_t_work(ind)
2973 write (iout,*) "Random velocities in the Calpha,SC space"
2975 write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3)
2978 write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3)
2985 ! if (itype(i,1).eq.10) then
2987 if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
2990 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
2994 d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
2995 d_t(j,i)=d_t(j,i+1)-d_t(j,i)
3000 write (iout,*)"Random velocities in the virtual-bond-vector space"
3002 write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3)
3005 write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3)
3008 write (iout,*) "Ek from d_t_work",Ek1
3009 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
3017 d_t_work(ind)=d_t_work(ind) &
3018 +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
3020 ! write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
3024 ! Transfer to the d_t vector
3026 d_t(j,0)=d_t_work(j)
3032 d_t(j,i)=d_t_work(ind)
3037 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3038 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3039 .and.(mnum.ne.5)) then
3042 d_t(j,i+nres)=d_t_work(ind)
3048 ! write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",&
3049 ! 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
3053 end subroutine random_vel
3054 !-----------------------------------------------------------------------------
3056 subroutine sd_verlet_p_setup
3057 ! Sets up the parameters of stochastic Verlet algorithm
3058 ! implicit real*8 (a-h,o-z)
3059 ! include 'DIMENSIONS'
3060 use control, only: tcpu
3065 ! include 'COMMON.CONTROL'
3066 ! include 'COMMON.VAR'
3067 ! include 'COMMON.MD'
3069 ! include 'COMMON.LANGEVIN'
3071 ! include 'COMMON.LANGEVIN.lang0'
3073 ! include 'COMMON.CHAIN'
3074 ! include 'COMMON.DERIV'
3075 ! include 'COMMON.GEO'
3076 ! include 'COMMON.LOCAL'
3077 ! include 'COMMON.INTERACT'
3078 ! include 'COMMON.IOUNITS'
3079 ! include 'COMMON.NAMES'
3080 ! include 'COMMON.TIME1'
3081 real(kind=8),dimension(6*nres) :: emgdt !(MAXRES6) maxres6=6*maxres
3082 real(kind=8) :: pterm,vterm,rho,rhoc,vsig
3083 real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
3084 prand_vec,vrand_vec1,vrand_vec2 !(MAXRES6) maxres6=6*maxres
3085 logical :: lprn = .false.
3086 real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
3087 real(kind=8) :: ktm,gdt,egdt,gdt2,gdt3,gdt4,gdt5,gdt6,gdt7,gdt8,&
3089 integer :: i,maxres2
3096 ! AL 8/17/04 Code adapted from tinker
3098 ! Get the frictional and random terms for stochastic dynamics in the
3099 ! eigenspace of mass-scaled UNRES friction matrix
3103 gdt = fricgam(i) * d_time
3105 ! Stochastic dynamics reduces to simple MD for zero friction
3107 if (gdt .le. zero) then
3108 pfric_vec(i) = 1.0d0
3109 vfric_vec(i) = d_time
3110 afric_vec(i) = 0.5d0 * d_time * d_time
3111 prand_vec(i) = 0.0d0
3112 vrand_vec1(i) = 0.0d0
3113 vrand_vec2(i) = 0.0d0
3115 ! Analytical expressions when friction coefficient is large
3118 if (gdt .ge. gdt_radius) then
3121 vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
3122 afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
3123 pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
3124 vterm = 1.0d0 - egdt**2
3125 rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
3127 ! Use series expansions when friction coefficient is small
3138 afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0 &
3139 - gdt5/120.0d0 + gdt6/720.0d0 &
3140 - gdt7/5040.0d0 + gdt8/40320.0d0 &
3141 - gdt9/362880.0d0) / fricgam(i)**2
3142 vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
3143 pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
3144 pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0 &
3145 + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0 &
3146 + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0 &
3147 + 127.0d0*gdt9/90720.0d0
3148 vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0 &
3149 - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0 &
3150 - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0 &
3151 - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
3152 rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0 &
3153 - 17.0d0*gdt2/1280.0d0 &
3154 + 17.0d0*gdt3/6144.0d0 &
3155 + 40967.0d0*gdt4/34406400.0d0 &
3156 - 57203.0d0*gdt5/275251200.0d0 &
3157 - 1429487.0d0*gdt6/13212057600.0d0)
3160 ! Compute the scaling factors of random terms for the nonzero friction case
3162 ktm = 0.5d0*d_time/fricgam(i)
3163 psig = dsqrt(ktm*pterm) / fricgam(i)
3164 vsig = dsqrt(ktm*vterm)
3165 rhoc = dsqrt(1.0d0 - rho*rho)
3167 vrand_vec1(i) = vsig * rho
3168 vrand_vec2(i) = vsig * rhoc
3173 "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",&
3176 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),&
3177 afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
3181 ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
3184 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
3185 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
3186 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
3187 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
3188 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
3189 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
3192 t_sdsetup=t_sdsetup+MPI_Wtime()
3194 t_sdsetup=t_sdsetup+tcpu()-tt0
3197 end subroutine sd_verlet_p_setup
3198 !-----------------------------------------------------------------------------
3199 subroutine eigtransf1(n,ndim,ab,d,c)
3203 real(kind=8) :: ab(ndim,ndim,n),c(ndim,n),d(ndim)
3209 c(i,j)=c(i,j)+ab(k,j,i)*d(k)
3214 end subroutine eigtransf1
3215 !-----------------------------------------------------------------------------
3216 subroutine eigtransf(n,ndim,a,b,d,c)
3220 real(kind=8) :: a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
3226 c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
3231 end subroutine eigtransf
3232 !-----------------------------------------------------------------------------
3233 subroutine sd_verlet1
3235 ! Applying stochastic velocity Verlet algorithm - step 1 to velocities
3237 ! implicit real*8 (a-h,o-z)
3238 ! include 'DIMENSIONS'
3239 ! include 'COMMON.CONTROL'
3240 ! include 'COMMON.VAR'
3241 ! include 'COMMON.MD'
3243 ! include 'COMMON.LANGEVIN'
3245 ! include 'COMMON.LANGEVIN.lang0'
3247 ! include 'COMMON.CHAIN'
3248 ! include 'COMMON.DERIV'
3249 ! include 'COMMON.GEO'
3250 ! include 'COMMON.LOCAL'
3251 ! include 'COMMON.INTERACT'
3252 ! include 'COMMON.IOUNITS'
3253 ! include 'COMMON.NAMES'
3254 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
3255 !el common /stochcalc/ stochforcvec
3256 logical :: lprn = .false.
3257 real(kind=8) :: ddt1,ddt2
3258 integer :: i,j,ind,inres
3260 ! write (iout,*) "dc_old"
3262 ! write (iout,'(i5,3f10.5,5x,3f10.5)')
3263 ! & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3266 dc_work(j)=dc_old(j,0)
3267 d_t_work(j)=d_t_old(j,0)
3268 d_a_work(j)=d_a_old(j,0)
3273 dc_work(ind+j)=dc_old(j,i)
3274 d_t_work(ind+j)=d_t_old(j,i)
3275 d_a_work(ind+j)=d_a_old(j,i)
3281 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3282 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3283 .and.(mnum.ne.5)) then
3285 dc_work(ind+j)=dc_old(j,i+nres)
3286 d_t_work(ind+j)=d_t_old(j,i+nres)
3287 d_a_work(ind+j)=d_a_old(j,i+nres)
3295 "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",&
3299 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
3300 vfric_mat(i,j),afric_mat(i,j),&
3301 prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3309 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) &
3310 +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3311 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3312 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3314 d_t_work_new(i)=ddt1+0.5d0*ddt2
3315 d_t_work(i)=ddt1+ddt2
3320 d_t(j,0)=d_t_work(j)
3325 dc(j,i)=dc_work(ind+j)
3326 d_t(j,i)=d_t_work(ind+j)
3332 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3333 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3334 .and.(mnum.ne.5)) then
3337 dc(j,inres)=dc_work(ind+j)
3338 d_t(j,inres)=d_t_work(ind+j)
3344 end subroutine sd_verlet1
3345 !-----------------------------------------------------------------------------
3346 subroutine sd_verlet2
3348 ! Calculating the adjusted velocities for accelerations
3350 ! implicit real*8 (a-h,o-z)
3351 ! include 'DIMENSIONS'
3352 ! include 'COMMON.CONTROL'
3353 ! include 'COMMON.VAR'
3354 ! include 'COMMON.MD'
3356 ! include 'COMMON.LANGEVIN'
3358 ! include 'COMMON.LANGEVIN.lang0'
3360 ! include 'COMMON.CHAIN'
3361 ! include 'COMMON.DERIV'
3362 ! include 'COMMON.GEO'
3363 ! include 'COMMON.LOCAL'
3364 ! include 'COMMON.INTERACT'
3365 ! include 'COMMON.IOUNITS'
3366 ! include 'COMMON.NAMES'
3367 !el real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV !(MAXRES6) maxres6=6*maxres
3368 real(kind=8),dimension(6*nres) :: stochforcvecV !(MAXRES6) maxres6=6*maxres
3369 !el common /stochcalc/ stochforcvec
3371 real(kind=8) :: ddt1,ddt2
3372 integer :: i,j,ind,inres
3373 ! Compute the stochastic forces which contribute to velocity change
3375 call stochastic_force(stochforcvecV)
3382 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3383 ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+ &
3384 vrand_mat2(i,j)*stochforcvecV(j)
3386 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
3390 d_t(j,0)=d_t_work(j)
3395 d_t(j,i)=d_t_work(ind+j)
3401 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3402 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3403 .and.(mnum.ne.5)) then
3406 d_t(j,inres)=d_t_work(ind+j)
3412 end subroutine sd_verlet2
3413 !-----------------------------------------------------------------------------
3414 subroutine sd_verlet_ciccotti_setup
3416 ! Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
3418 ! implicit real*8 (a-h,o-z)
3419 ! include 'DIMENSIONS'
3420 use control, only: tcpu
3425 ! include 'COMMON.CONTROL'
3426 ! include 'COMMON.VAR'
3427 ! include 'COMMON.MD'
3429 ! include 'COMMON.LANGEVIN'
3431 ! include 'COMMON.LANGEVIN.lang0'
3433 ! include 'COMMON.CHAIN'
3434 ! include 'COMMON.DERIV'
3435 ! include 'COMMON.GEO'
3436 ! include 'COMMON.LOCAL'
3437 ! include 'COMMON.INTERACT'
3438 ! include 'COMMON.IOUNITS'
3439 ! include 'COMMON.NAMES'
3440 ! include 'COMMON.TIME1'
3441 real(kind=8),dimension(6*nres) :: emgdt !(MAXRES6) maxres6=6*maxres
3442 real(kind=8) :: pterm,vterm,rho,rhoc,vsig
3443 real(kind=8),dimension(6*nres) :: pfric_vec,vfric_vec,afric_vec,&
3444 prand_vec,vrand_vec1,vrand_vec2 !(MAXRES6) maxres6=6*maxres
3445 logical :: lprn = .false.
3446 real(kind=8) :: zero = 1.0d-8, gdt_radius = 0.05d0
3447 real(kind=8) :: ktm,gdt,egdt,tt0
3448 integer :: i,maxres2
3455 ! AL 8/17/04 Code adapted from tinker
3457 ! Get the frictional and random terms for stochastic dynamics in the
3458 ! eigenspace of mass-scaled UNRES friction matrix
3462 write (iout,*) "i",i," fricgam",fricgam(i)
3463 gdt = fricgam(i) * d_time
3465 ! Stochastic dynamics reduces to simple MD for zero friction
3467 if (gdt .le. zero) then
3468 pfric_vec(i) = 1.0d0
3469 vfric_vec(i) = d_time
3470 afric_vec(i) = 0.5d0*d_time*d_time
3471 prand_vec(i) = afric_vec(i)
3472 vrand_vec2(i) = vfric_vec(i)
3474 ! Analytical expressions when friction coefficient is large
3479 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
3480 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
3481 prand_vec(i) = afric_vec(i)
3482 vrand_vec2(i) = vfric_vec(i)
3484 ! Compute the scaling factors of random terms for the nonzero friction case
3486 ! ktm = 0.5d0*d_time/fricgam(i)
3487 ! psig = dsqrt(ktm*pterm) / fricgam(i)
3488 ! vsig = dsqrt(ktm*vterm)
3489 ! prand_vec(i) = psig*afric_vec(i)
3490 ! vrand_vec2(i) = vsig*vfric_vec(i)
3495 "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",&
3498 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),&
3499 afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
3503 ! Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
3505 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
3506 call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
3507 call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
3508 call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
3509 call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
3511 t_sdsetup=t_sdsetup+MPI_Wtime()
3513 t_sdsetup=t_sdsetup+tcpu()-tt0
3516 end subroutine sd_verlet_ciccotti_setup
3517 !-----------------------------------------------------------------------------
3518 subroutine sd_verlet1_ciccotti
3520 ! Applying stochastic velocity Verlet algorithm - step 1 to velocities
3521 ! implicit real*8 (a-h,o-z)
3523 ! include 'DIMENSIONS'
3527 ! include 'COMMON.CONTROL'
3528 ! include 'COMMON.VAR'
3529 ! include 'COMMON.MD'
3531 ! include 'COMMON.LANGEVIN'
3533 ! include 'COMMON.LANGEVIN.lang0'
3535 ! include 'COMMON.CHAIN'
3536 ! include 'COMMON.DERIV'
3537 ! include 'COMMON.GEO'
3538 ! include 'COMMON.LOCAL'
3539 ! include 'COMMON.INTERACT'
3540 ! include 'COMMON.IOUNITS'
3541 ! include 'COMMON.NAMES'
3542 !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
3543 !el common /stochcalc/ stochforcvec
3544 logical :: lprn = .false.
3545 real(kind=8) :: ddt1,ddt2
3546 integer :: i,j,ind,inres
3547 ! write (iout,*) "dc_old"
3549 ! write (iout,'(i5,3f10.5,5x,3f10.5)')
3550 ! & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
3553 dc_work(j)=dc_old(j,0)
3554 d_t_work(j)=d_t_old(j,0)
3555 d_a_work(j)=d_a_old(j,0)
3560 dc_work(ind+j)=dc_old(j,i)
3561 d_t_work(ind+j)=d_t_old(j,i)
3562 d_a_work(ind+j)=d_a_old(j,i)
3567 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3569 dc_work(ind+j)=dc_old(j,i+nres)
3570 d_t_work(ind+j)=d_t_old(j,i+nres)
3571 d_a_work(ind+j)=d_a_old(j,i+nres)
3580 "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",&
3584 write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
3585 vfric_mat(i,j),afric_mat(i,j),&
3586 prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
3594 dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j) &
3595 +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
3596 ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
3597 ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
3599 d_t_work_new(i)=ddt1+0.5d0*ddt2
3600 d_t_work(i)=ddt1+ddt2
3605 d_t(j,0)=d_t_work(j)
3610 dc(j,i)=dc_work(ind+j)
3611 d_t(j,i)=d_t_work(ind+j)
3617 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3618 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3619 .and.(mnum.ne.5)) then
3622 dc(j,inres)=dc_work(ind+j)
3623 d_t(j,inres)=d_t_work(ind+j)
3629 end subroutine sd_verlet1_ciccotti
3630 !-----------------------------------------------------------------------------
3631 subroutine sd_verlet2_ciccotti
3633 ! Calculating the adjusted velocities for accelerations
3635 ! implicit real*8 (a-h,o-z)
3636 ! include 'DIMENSIONS'
3637 ! include 'COMMON.CONTROL'
3638 ! include 'COMMON.VAR'
3639 ! include 'COMMON.MD'
3641 ! include 'COMMON.LANGEVIN'
3643 ! include 'COMMON.LANGEVIN.lang0'
3645 ! include 'COMMON.CHAIN'
3646 ! include 'COMMON.DERIV'
3647 ! include 'COMMON.GEO'
3648 ! include 'COMMON.LOCAL'
3649 ! include 'COMMON.INTERACT'
3650 ! include 'COMMON.IOUNITS'
3651 ! include 'COMMON.NAMES'
3652 !el real(kind=8),dimension(6*nres) :: stochforcvec,stochforcvecV !(MAXRES6) maxres6=6*maxres
3653 real(kind=8),dimension(6*nres) :: stochforcvecV !(MAXRES6) maxres6=6*maxres
3654 !el common /stochcalc/ stochforcvec
3655 real(kind=8) :: ddt1,ddt2
3656 integer :: i,j,ind,inres
3658 ! Compute the stochastic forces which contribute to velocity change
3660 call stochastic_force(stochforcvecV)
3667 ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
3668 ! ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
3669 ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
3671 d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
3675 d_t(j,0)=d_t_work(j)
3680 d_t(j,i)=d_t_work(ind+j)
3686 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3688 ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3691 d_t(j,inres)=d_t_work(ind+j)
3697 end subroutine sd_verlet2_ciccotti
3699 !-----------------------------------------------------------------------------
3701 !-----------------------------------------------------------------------------
3702 subroutine inertia_tensor
3704 ! Calculating the intertia tensor for the entire protein in order to
3705 ! remove the perpendicular components of velocity matrix which cause
3706 ! the molecule to rotate.
3709 ! implicit real*8 (a-h,o-z)
3710 ! include 'DIMENSIONS'
3711 ! include 'COMMON.CONTROL'
3712 ! include 'COMMON.VAR'
3713 ! include 'COMMON.MD'
3714 ! include 'COMMON.CHAIN'
3715 ! include 'COMMON.DERIV'
3716 ! include 'COMMON.GEO'
3717 ! include 'COMMON.LOCAL'
3718 ! include 'COMMON.INTERACT'
3719 ! include 'COMMON.IOUNITS'
3720 ! include 'COMMON.NAMES'
3722 real(kind=8),dimension(3,3) :: Im,Imcp,eigvec,Id
3723 real(kind=8),dimension(3) :: pr,eigval,L,vp,vrot
3724 real(kind=8) :: M_SC,mag,mag2,M_PEP
3725 real(kind=8),dimension(3,0:nres) :: vpp !(3,0:MAXRES)
3726 real(kind=8),dimension(3) :: vs_p,pp,incr,v
3727 real(kind=8),dimension(3,3) :: pr1,pr2
3729 !el common /gucio/ cm
3730 integer :: iti,inres,i,j,k,mnum
3741 ! calculating the center of the mass of the protein
3745 if (itype(i,mnum).eq.ntyp1_molec(mnum)) cycle
3746 M_PEP=M_PEP+mp(mnum)
3748 cm(j)=cm(j)+(c(j,i)+0.5d0*dc(j,i))*mp(mnum)
3757 iti=iabs(itype(i,mnum))
3758 M_SC=M_SC+msc(iabs(iti),mnum)
3761 cm(j)=cm(j)+msc(iabs(iti),mnum)*c(j,inres)
3765 cm(j)=cm(j)/(M_SC+M_PEP)
3771 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
3773 Im(1,1)=Im(1,1)+mp(mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
3774 Im(1,2)=Im(1,2)-mp(mnum)*pr(1)*pr(2)
3775 Im(1,3)=Im(1,3)-mp(mnum)*pr(1)*pr(3)
3776 Im(2,3)=Im(2,3)-mp(mnum)*pr(2)*pr(3)
3777 Im(2,2)=Im(2,2)+mp(mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
3778 Im(3,3)=Im(3,3)+mp(mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
3783 iti=iabs(itype(i,mnum))
3786 pr(j)=c(j,inres)-cm(j)
3788 Im(1,1)=Im(1,1)+msc(iabs(iti),mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
3789 Im(1,2)=Im(1,2)-msc(iabs(iti),mnum)*pr(1)*pr(2)
3790 Im(1,3)=Im(1,3)-msc(iabs(iti),mnum)*pr(1)*pr(3)
3791 Im(2,3)=Im(2,3)-msc(iabs(iti),mnum)*pr(2)*pr(3)
3792 Im(2,2)=Im(2,2)+msc(iabs(iti),mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
3793 Im(3,3)=Im(3,3)+msc(iabs(iti),mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
3798 Im(1,1)=Im(1,1)+Ip(mnum)*(1-dc_norm(1,i)*dc_norm(1,i))* &
3799 vbld(i+1)*vbld(i+1)*0.25d0
3800 Im(1,2)=Im(1,2)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(2,i))* &
3801 vbld(i+1)*vbld(i+1)*0.25d0
3802 Im(1,3)=Im(1,3)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(3,i))* &
3803 vbld(i+1)*vbld(i+1)*0.25d0
3804 Im(2,3)=Im(2,3)+Ip(mnum)*(-dc_norm(2,i)*dc_norm(3,i))* &
3805 vbld(i+1)*vbld(i+1)*0.25d0
3806 Im(2,2)=Im(2,2)+Ip(mnum)*(1-dc_norm(2,i)*dc_norm(2,i))* &
3807 vbld(i+1)*vbld(i+1)*0.25d0
3808 Im(3,3)=Im(3,3)+Ip(mnum)*(1-dc_norm(3,i)*dc_norm(3,i))* &
3809 vbld(i+1)*vbld(i+1)*0.25d0
3815 ! if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
3816 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3817 .and.(mnum.ne.5)) then
3818 iti=iabs(itype(i,mnum))
3820 Im(1,1)=Im(1,1)+Isc(iti,mnum)*(1-dc_norm(1,inres)* &
3821 dc_norm(1,inres))*vbld(inres)*vbld(inres)
3822 Im(1,2)=Im(1,2)-Isc(iti,mnum)*(dc_norm(1,inres)* &
3823 dc_norm(2,inres))*vbld(inres)*vbld(inres)
3824 Im(1,3)=Im(1,3)-Isc(iti,mnum)*(dc_norm(1,inres)* &
3825 dc_norm(3,inres))*vbld(inres)*vbld(inres)
3826 Im(2,3)=Im(2,3)-Isc(iti,mnum)*(dc_norm(2,inres)* &
3827 dc_norm(3,inres))*vbld(inres)*vbld(inres)
3828 Im(2,2)=Im(2,2)+Isc(iti,mnum)*(1-dc_norm(2,inres)* &
3829 dc_norm(2,inres))*vbld(inres)*vbld(inres)
3830 Im(3,3)=Im(3,3)+Isc(iti,mnum)*(1-dc_norm(3,inres)* &
3831 dc_norm(3,inres))*vbld(inres)*vbld(inres)
3836 ! write(iout,*) "The angular momentum before adjustment:"
3837 ! write(iout,*) (L(j),j=1,3)
3843 ! Copying the Im matrix for the djacob subroutine
3851 ! Finding the eigenvectors and eignvalues of the inertia tensor
3852 call djacob(3,3,10000,1.0d-10,Imcp,eigvec,eigval)
3853 ! write (iout,*) "Eigenvalues & Eigenvectors"
3854 ! write (iout,'(5x,3f10.5)') (eigval(i),i=1,3)
3857 ! write (iout,'(i5,3f10.5)') i,(eigvec(i,j),j=1,3)
3859 ! Constructing the diagonalized matrix
3861 if (dabs(eigval(i)).gt.1.0d-15) then
3862 Id(i,i)=1.0d0/eigval(i)
3869 Imcp(i,j)=eigvec(j,i)
3875 pr1(i,j)=pr1(i,j)+Id(i,k)*Imcp(k,j)
3882 pr2(i,j)=pr2(i,j)+eigvec(i,k)*pr1(k,j)
3886 ! Calculating the total rotational velocity of the molecule
3889 vrot(i)=vrot(i)+pr2(i,j)*L(j)
3892 ! Resetting the velocities
3894 call vecpr(vrot(1),dc(1,i),vp)
3896 d_t(j,i)=d_t(j,i)-vp(j)
3901 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3902 .and.(mnum.ne.5)) then
3903 ! if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
3904 ! if(itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
3906 call vecpr(vrot(1),dc(1,inres),vp)
3908 d_t(j,inres)=d_t(j,inres)-vp(j)
3913 ! write(iout,*) "The angular momentum after adjustment:"
3914 ! write(iout,*) (L(j),j=1,3)
3917 end subroutine inertia_tensor
3918 !-----------------------------------------------------------------------------
3919 subroutine angmom(cm,L)
3922 ! implicit real*8 (a-h,o-z)
3923 ! include 'DIMENSIONS'
3924 ! include 'COMMON.CONTROL'
3925 ! include 'COMMON.VAR'
3926 ! include 'COMMON.MD'
3927 ! include 'COMMON.CHAIN'
3928 ! include 'COMMON.DERIV'
3929 ! include 'COMMON.GEO'
3930 ! include 'COMMON.LOCAL'
3931 ! include 'COMMON.INTERACT'
3932 ! include 'COMMON.IOUNITS'
3933 ! include 'COMMON.NAMES'
3935 real(kind=8),dimension(3) :: L,cm,pr,vp,vrot,incr,v,pp
3936 integer :: iti,inres,i,j,mnum
3937 ! Calculate the angular momentum
3947 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
3950 v(j)=incr(j)+0.5d0*d_t(j,i)
3953 incr(j)=incr(j)+d_t(j,i)
3955 call vecpr(pr(1),v(1),vp)
3957 L(j)=L(j)+mp(mnum)*vp(j)
3961 pp(j)=0.5d0*d_t(j,i)
3963 call vecpr(pr(1),pp(1),vp)
3965 L(j)=L(j)+Ip(mnum)*vp(j)
3973 iti=iabs(itype(i,mnum))
3976 pr(j)=c(j,inres)-cm(j)
3978 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3979 .and.(mnum.ne.5)) then
3981 v(j)=incr(j)+d_t(j,inres)
3988 call vecpr(pr(1),v(1),vp)
3989 ! write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),
3990 ! & " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
3992 L(j)=L(j)+msc(iabs(iti),mnum)*vp(j)
3994 ! write (iout,*) "L",(l(j),j=1,3)
3995 if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
3996 .and.(mnum.ne.5)) then
3998 v(j)=incr(j)+d_t(j,inres)
4000 call vecpr(dc(1,inres),d_t(1,inres),vp)
4002 L(j)=L(j)+Isc(iti,mnum)*vp(j)
4006 incr(j)=incr(j)+d_t(j,i)
4010 end subroutine angmom
4011 !-----------------------------------------------------------------------------
4012 subroutine vcm_vel(vcm)
4015 ! implicit real*8 (a-h,o-z)
4016 ! include 'DIMENSIONS'
4017 ! include 'COMMON.VAR'
4018 ! include 'COMMON.MD'
4019 ! include 'COMMON.CHAIN'
4020 ! include 'COMMON.DERIV'
4021 ! include 'COMMON.GEO'
4022 ! include 'COMMON.LOCAL'
4023 ! include 'COMMON.INTERACT'
4024 ! include 'COMMON.IOUNITS'
4025 real(kind=8),dimension(3) :: vcm,vv
4026 real(kind=8) :: summas,amas
4037 summas=summas+mp(mnum)
4039 vcm(j)=vcm(j)+mp(mnum)*(vv(j)+0.5d0*d_t(j,i))
4042 amas=msc(iabs(itype(i,mnum)),mnum)
4044 if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4045 .and.(mnum.ne.5)) then
4047 vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres))
4051 vcm(j)=vcm(j)+amas*vv(j)
4055 vv(j)=vv(j)+d_t(j,i)
4058 ! write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas
4060 vcm(j)=vcm(j)/summas
4063 end subroutine vcm_vel
4064 !-----------------------------------------------------------------------------
4066 !-----------------------------------------------------------------------------
4068 ! RATTLE algorithm for velocity Verlet - step 1, UNRES
4072 ! implicit real*8 (a-h,o-z)
4073 ! include 'DIMENSIONS'
4075 ! include 'COMMON.CONTROL'
4076 ! include 'COMMON.VAR'
4077 ! include 'COMMON.MD'
4079 ! include 'COMMON.LANGEVIN'
4081 ! include 'COMMON.LANGEVIN.lang0'
4083 ! include 'COMMON.CHAIN'
4084 ! include 'COMMON.DERIV'
4085 ! include 'COMMON.GEO'
4086 ! include 'COMMON.LOCAL'
4087 ! include 'COMMON.INTERACT'
4088 ! include 'COMMON.IOUNITS'
4089 ! include 'COMMON.NAMES'
4090 ! include 'COMMON.TIME1'
4091 !el real(kind=8) :: gginv(2*nres,2*nres),&
4092 !el gdc(3,2*nres,2*nres)
4093 real(kind=8) :: dC_uncor(3,2*nres) !,&
4094 !el real(kind=8) :: Cmat(2*nres,2*nres)
4095 real(kind=8) :: x(2*nres),xcorr(3,2*nres) !maxres2=2*maxres
4096 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
4097 !el common /przechowalnia/ nbond
4098 integer :: max_rattle = 5
4099 logical :: lprn = .false., lprn1 = .false., not_done
4100 real(kind=8) :: tol_rattle = 1.0d-5
4102 integer :: ii,i,j,jj,l,ind,ind1,nres2
4105 !el /common/ przechowalnia
4107 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4108 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4109 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4111 if (lprn) write (iout,*) "RATTLE1"
4115 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4116 .and.(mnum.ne.5)) nbond=nbond+1
4118 ! Make a folded form of the Ginv-matrix
4131 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4136 if (itype(k,1).ne.10 .and. itype(k,mnum).ne.ntyp1_molec(mnum)&
4137 .and.(mnum.ne.5)) then
4141 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4149 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4160 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4164 if (itype(k,1).ne.10) then
4168 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
4176 write (iout,*) "Matrix GGinv"
4177 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
4183 if (iter.gt.max_rattle) then
4184 write (iout,*) "Error - too many iterations in RATTLE."
4187 ! Calculate the matrix C = GG**(-1) dC_old o dC
4192 dC_uncor(j,ind1)=dC(j,i)
4196 if (itype(i,1).ne.10) then
4199 dC_uncor(j,ind1)=dC(j,i+nres)
4208 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
4212 if (itype(k,1).ne.10) then
4215 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
4220 ! Calculate deviations from standard virtual-bond lengths
4224 x(ind)=vbld(i+1)**2-vbl**2
4227 if (itype(i,1).ne.10) then
4229 x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4233 write (iout,*) "Coordinates and violations"
4235 write(iout,'(i5,3f10.5,5x,e15.5)') &
4236 i,(dC_uncor(j,i),j=1,3),x(i)
4238 write (iout,*) "Velocities and violations"
4242 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4243 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4247 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4248 .and.(mnum.ne.5)) then
4251 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4252 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4253 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4256 ! write (iout,*) "gdc"
4258 ! write (iout,*) "i",i
4260 ! write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4266 if (dabs(x(i)).gt.xmax) then
4270 if (xmax.lt.tol_rattle) then
4274 ! Calculate the matrix of the system of equations
4279 Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
4284 write (iout,*) "Matrix Cmat"
4285 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4287 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4288 ! Add constraint term to positions
4295 xx = xx+x(ii)*gdc(j,ind,ii)
4299 d_t_new(j,i)=d_t_new(j,i)-xx/d_time
4303 if (itype(i,1).ne.10) then
4308 xx = xx+x(ii)*gdc(j,ind,ii)
4311 dC(j,i+nres)=dC(j,i+nres)-xx
4312 d_t_new(j,i+nres)=d_t_new(j,i+nres)-xx/d_time
4316 ! Rebuild the chain using the new coordinates
4317 call chainbuild_cart
4319 write (iout,*) "New coordinates, Lagrange multipliers,",&
4320 " and differences between actual and standard bond lengths"
4324 xx=vbld(i+1)**2-vbl**2
4325 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4326 i,(dC(j,i),j=1,3),x(ind),xx
4330 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4333 xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4334 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4335 i,(dC(j,i+nres),j=1,3),x(ind),xx
4338 write (iout,*) "Velocities and violations"
4342 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4343 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4346 if (itype(i,1).ne.10) then
4348 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4349 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4350 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4357 10 write (iout,*) "Error - singularity in solving the system",&
4358 " of equations for Lagrange multipliers."
4362 "RATTLE inactive; use -DRATTLE switch at compile time."
4365 end subroutine rattle1
4366 !-----------------------------------------------------------------------------
4368 ! RATTLE algorithm for velocity Verlet - step 2, UNRES
4372 ! implicit real*8 (a-h,o-z)
4373 ! include 'DIMENSIONS'
4375 ! include 'COMMON.CONTROL'
4376 ! include 'COMMON.VAR'
4377 ! include 'COMMON.MD'
4379 ! include 'COMMON.LANGEVIN'
4381 ! include 'COMMON.LANGEVIN.lang0'
4383 ! include 'COMMON.CHAIN'
4384 ! include 'COMMON.DERIV'
4385 ! include 'COMMON.GEO'
4386 ! include 'COMMON.LOCAL'
4387 ! include 'COMMON.INTERACT'
4388 ! include 'COMMON.IOUNITS'
4389 ! include 'COMMON.NAMES'
4390 ! include 'COMMON.TIME1'
4391 !el real(kind=8) :: gginv(2*nres,2*nres),&
4392 !el gdc(3,2*nres,2*nres)
4393 real(kind=8) :: dC_uncor(3,2*nres) !,&
4394 !el Cmat(2*nres,2*nres)
4395 real(kind=8) :: x(2*nres) !maxres2=2*maxres
4396 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
4397 !el common /przechowalnia/ nbond
4398 integer :: max_rattle = 5
4399 logical :: lprn = .false., lprn1 = .false., not_done
4400 real(kind=8) :: tol_rattle = 1.0d-5
4404 !el /common/ przechowalnia
4405 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4406 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4407 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4409 if (lprn) write (iout,*) "RATTLE2"
4410 if (lprn) write (iout,*) "Velocity correction"
4411 ! Calculate the matrix G dC
4417 gdc(j,i,ind)=GGinv(i,ind)*dC(j,k)
4422 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4423 .and.(mnum.ne.5)) then
4426 gdc(j,i,ind)=GGinv(i,ind)*dC(j,k+nres)
4432 ! write (iout,*) "gdc"
4434 ! write (iout,*) "i",i
4436 ! write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4440 ! Calculate the matrix of the system of equations
4447 Cmat(ind,j)=Cmat(ind,j)+dC(k,i)*gdc(k,ind,j)
4453 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4454 .and.(mnum.ne.5)) then
4459 Cmat(ind,j)=Cmat(ind,j)+dC(k,i+nres)*gdc(k,ind,j)
4464 ! Calculate the scalar product dC o d_t_new
4468 x(ind)=scalar(d_t(1,i),dC(1,i))
4472 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4473 .and.(mnum.ne.5)) then
4475 x(ind)=scalar(d_t(1,i+nres),dC(1,i+nres))
4479 write (iout,*) "Velocities and violations"
4483 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4484 i,ind,(d_t(j,i),j=1,3),x(ind)
4488 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4489 .and.(mnum.ne.5)) then
4491 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4492 i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind)
4498 if (dabs(x(i)).gt.xmax) then
4502 if (xmax.lt.tol_rattle) then
4507 write (iout,*) "Matrix Cmat"
4508 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4510 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4511 ! Add constraint term to velocities
4518 xx = xx+x(ii)*gdc(j,ind,ii)
4520 d_t(j,i)=d_t(j,i)-xx
4525 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4526 .and.(mnum.ne.5)) then
4531 xx = xx+x(ii)*gdc(j,ind,ii)
4533 d_t(j,i+nres)=d_t(j,i+nres)-xx
4539 "New velocities, Lagrange multipliers violations"
4543 if (lprn) write (iout,'(2i5,3f10.5,5x,2e15.5)') &
4544 i,ind,(d_t(j,i),j=1,3),x(ind),scalar(d_t(1,i),dC(1,i))
4548 if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
4551 write (iout,'(2i5,3f10.5,5x,2e15.5)') &
4552 i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind),&
4553 scalar(d_t(1,i+nres),dC(1,i+nres))
4559 10 write (iout,*) "Error - singularity in solving the system",&
4560 " of equations for Lagrange multipliers."
4564 "RATTLE inactive; use -DRATTLE option at compile time."
4567 end subroutine rattle2
4568 !-----------------------------------------------------------------------------
4569 subroutine rattle_brown
4570 ! RATTLE/LINCS algorithm for Brownian dynamics, UNRES
4574 ! implicit real*8 (a-h,o-z)
4575 ! include 'DIMENSIONS'
4577 ! include 'COMMON.CONTROL'
4578 ! include 'COMMON.VAR'
4579 ! include 'COMMON.MD'
4581 ! include 'COMMON.LANGEVIN'
4583 ! include 'COMMON.LANGEVIN.lang0'
4585 ! include 'COMMON.CHAIN'
4586 ! include 'COMMON.DERIV'
4587 ! include 'COMMON.GEO'
4588 ! include 'COMMON.LOCAL'
4589 ! include 'COMMON.INTERACT'
4590 ! include 'COMMON.IOUNITS'
4591 ! include 'COMMON.NAMES'
4592 ! include 'COMMON.TIME1'
4593 !el real(kind=8) :: gginv(2*nres,2*nres),&
4594 !el gdc(3,2*nres,2*nres)
4595 real(kind=8) :: dC_uncor(3,2*nres) !,&
4596 !el real(kind=8) :: Cmat(2*nres,2*nres)
4597 real(kind=8) :: x(2*nres) !maxres2=2*maxres
4598 !el common /przechowalnia/ GGinv,gdc,Cmat,nbond
4599 !el common /przechowalnia/ nbond
4600 integer :: max_rattle = 5
4601 logical :: lprn = .false., lprn1 = .false., not_done
4602 real(kind=8) :: tol_rattle = 1.0d-5
4606 !el /common/ przechowalnia
4607 if(.not.allocated(GGinv)) allocate(GGinv(nres2,nres2))
4608 if(.not.allocated(gdc)) allocate(gdc(3,nres2,nres2))
4609 if(.not.allocated(Cmat)) allocate(Cmat(nres2,nres2))
4612 if (lprn) write (iout,*) "RATTLE_BROWN"
4615 if (itype(i,1).ne.10) nbond=nbond+1
4617 ! Make a folded form of the Ginv-matrix
4630 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4634 if (itype(k,1).ne.10) then
4638 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4645 if (itype(i,1).ne.10) then
4655 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
4659 if (itype(k,1).ne.10) then
4663 if (j.eq.1 .and. l.eq.1)GGinv(ii,jj)=fricmat(ind,ind1)
4671 write (iout,*) "Matrix GGinv"
4672 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
4678 if (iter.gt.max_rattle) then
4679 write (iout,*) "Error - too many iterations in RATTLE."
4682 ! Calculate the matrix C = GG**(-1) dC_old o dC
4687 dC_uncor(j,ind1)=dC(j,i)
4691 if (itype(i,1).ne.10) then
4694 dC_uncor(j,ind1)=dC(j,i+nres)
4703 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
4707 if (itype(k,1).ne.10) then
4710 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
4715 ! Calculate deviations from standard virtual-bond lengths
4719 x(ind)=vbld(i+1)**2-vbl**2
4722 if (itype(i,1).ne.10) then
4724 x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4728 write (iout,*) "Coordinates and violations"
4730 write(iout,'(i5,3f10.5,5x,e15.5)') &
4731 i,(dC_uncor(j,i),j=1,3),x(i)
4733 write (iout,*) "Velocities and violations"
4737 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4738 i,ind,(d_t(j,i),j=1,3),scalar(d_t(1,i),dC_old(1,i))
4741 if (itype(i,1).ne.10) then
4743 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4744 i+nres,ind,(d_t(j,i+nres),j=1,3),&
4745 scalar(d_t(1,i+nres),dC_old(1,i+nres))
4748 write (iout,*) "gdc"
4750 write (iout,*) "i",i
4752 write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
4758 if (dabs(x(i)).gt.xmax) then
4762 if (xmax.lt.tol_rattle) then
4766 ! Calculate the matrix of the system of equations
4771 Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
4776 write (iout,*) "Matrix Cmat"
4777 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
4779 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
4780 ! Add constraint term to positions
4787 xx = xx+x(ii)*gdc(j,ind,ii)
4790 d_t(j,i)=d_t(j,i)+xx/d_time
4795 if (itype(i,1).ne.10) then
4800 xx = xx+x(ii)*gdc(j,ind,ii)
4803 d_t(j,i+nres)=d_t(j,i+nres)+xx/d_time
4804 dC(j,i+nres)=dC(j,i+nres)+xx
4808 ! Rebuild the chain using the new coordinates
4809 call chainbuild_cart
4811 write (iout,*) "New coordinates, Lagrange multipliers,",&
4812 " and differences between actual and standard bond lengths"
4816 xx=vbld(i+1)**2-vbl**2
4817 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4818 i,(dC(j,i),j=1,3),x(ind),xx
4821 if (itype(i,1).ne.10) then
4823 xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
4824 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
4825 i,(dC(j,i+nres),j=1,3),x(ind),xx
4828 write (iout,*) "Velocities and violations"
4832 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4833 i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
4836 if (itype(i,1).ne.10) then
4838 write (iout,'(2i5,3f10.5,5x,e15.5)') &
4839 i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
4840 scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
4847 10 write (iout,*) "Error - singularity in solving the system",&
4848 " of equations for Lagrange multipliers."
4852 "RATTLE inactive; use -DRATTLE option at compile time"
4855 end subroutine rattle_brown
4856 !-----------------------------------------------------------------------------
4858 !-----------------------------------------------------------------------------
4859 subroutine friction_force
4864 ! implicit real*8 (a-h,o-z)
4865 ! include 'DIMENSIONS'
4866 ! include 'COMMON.VAR'
4867 ! include 'COMMON.CHAIN'
4868 ! include 'COMMON.DERIV'
4869 ! include 'COMMON.GEO'
4870 ! include 'COMMON.LOCAL'
4871 ! include 'COMMON.INTERACT'
4872 ! include 'COMMON.MD'
4874 ! include 'COMMON.LANGEVIN'
4876 ! include 'COMMON.LANGEVIN.lang0'
4878 ! include 'COMMON.IOUNITS'
4879 !el real(kind=8),dimension(6*nres) :: gamvec !(MAXRES6) maxres6=6*maxres
4880 !el common /syfek/ gamvec
4881 real(kind=8) :: vv(3),vvtot(3,nres),v_work(6*nres) !,&
4882 !el ginvfric(2*nres,2*nres) !maxres2=2*maxres
4883 !el common /przechowalnia/ ginvfric
4885 logical :: lprn = .false., checkmode = .false.
4886 integer :: i,j,ind,k,nres2,nres6,mnum
4890 if(.not.allocated(gamvec)) allocate(gamvec(nres6)) !(MAXRES6)
4891 if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
4899 d_t_work(j)=d_t(j,0)
4904 d_t_work(ind+j)=d_t(j,i)
4910 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
4911 .and.(mnum.ne.5)) then
4913 d_t_work(ind+j)=d_t(j,i+nres)
4919 call fricmat_mult(d_t_work,fric_work)
4921 if (.not.checkmode) return
4924 write (iout,*) "d_t_work and fric_work"
4926 write (iout,'(i3,2e15.5)') i,d_t_work(i),fric_work(i)
4930 friction(j,0)=fric_work(j)
4935 friction(j,i)=fric_work(ind+j)
4941 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
4942 .and.(mnum.ne.5)) then
4943 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
4945 friction(j,i+nres)=fric_work(ind+j)
4951 write(iout,*) "Friction backbone"
4953 write(iout,'(i5,3e15.5,5x,3e15.5)') &
4954 i,(friction(j,i),j=1,3),(d_t(j,i),j=1,3)
4956 write(iout,*) "Friction side chain"
4958 write(iout,'(i5,3e15.5,5x,3e15.5)') &
4959 i,(friction(j,i+nres),j=1,3),(d_t(j,i+nres),j=1,3)
4968 vvtot(j,i)=vv(j)+0.5d0*d_t(j,i)
4969 vvtot(j,i+nres)=vv(j)+d_t(j,i+nres)
4970 vv(j)=vv(j)+d_t(j,i)
4973 write (iout,*) "vvtot backbone and sidechain"
4975 write (iout,'(i5,3e15.5,5x,3e15.5)') i,(vvtot(j,i),j=1,3),&
4976 (vvtot(j,i+nres),j=1,3)
4981 v_work(ind+j)=vvtot(j,i)
4987 v_work(ind+j)=vvtot(j,i+nres)
4991 write (iout,*) "v_work gamvec and site-based friction forces"
4993 write (iout,'(i5,3e15.5)') i,v_work(i),gamvec(i),&
4997 ! fric_work1(i)=0.0d0
4999 ! fric_work1(i)=fric_work1(i)-A(j,i)*gamvec(j)*v_work(j)
5002 ! write (iout,*) "fric_work and fric_work1"
5004 ! write (iout,'(i5,2e15.5)') i,fric_work(i),fric_work1(i)
5010 ginvfric(i,j)=ginvfric(i,j)+ginv(i,k)*fricmat(k,j)
5014 write (iout,*) "ginvfric"
5016 write (iout,'(i5,100f8.3)') i,(ginvfric(i,j),j=1,dimen)
5018 write (iout,*) "symmetry check"
5021 write (iout,*) i,j,ginvfric(i,j)-ginvfric(j,i)
5026 end subroutine friction_force
5027 !-----------------------------------------------------------------------------
5028 subroutine setup_fricmat
5032 use control_data, only:time_Bcast
5033 use control, only:tcpu
5035 ! implicit real*8 (a-h,o-z)
5039 real(kind=8) :: time00
5041 ! include 'DIMENSIONS'
5042 ! include 'COMMON.VAR'
5043 ! include 'COMMON.CHAIN'
5044 ! include 'COMMON.DERIV'
5045 ! include 'COMMON.GEO'
5046 ! include 'COMMON.LOCAL'
5047 ! include 'COMMON.INTERACT'
5048 ! include 'COMMON.MD'
5049 ! include 'COMMON.SETUP'
5050 ! include 'COMMON.TIME1'
5051 ! integer licznik /0/
5054 ! include 'COMMON.LANGEVIN'
5056 ! include 'COMMON.LANGEVIN.lang0'
5058 ! include 'COMMON.IOUNITS'
5060 integer :: i,j,ind,ind1,m
5061 logical :: lprn = .false.
5062 real(kind=8) :: dtdi !el ,gamvec(2*nres)
5063 !el real(kind=8),dimension(2*nres,2*nres) :: ginvfric,fcopy
5064 real(kind=8),dimension(2*nres,2*nres) :: fcopy
5065 !el real(kind=8),dimension(2*nres*(2*nres+1)/2) :: Ghalf !(mmaxres2) (mmaxres2=(maxres2*(maxres2+1)/2))
5066 !el common /syfek/ gamvec
5067 real(kind=8) :: work(8*2*nres)
5068 integer :: iwork(2*nres)
5069 !el common /przechowalnia/ ginvfric,Ghalf,fcopy
5070 integer :: ii,iti,k,l,nzero,nres2,nres6,ierr,mnum
5072 if (fg_rank.ne.king) goto 10
5077 if(.not.allocated(gamvec)) allocate(gamvec(nres2)) !(MAXRES2)
5078 if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
5079 !el if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
5080 !el allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
5081 if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !maxres2=2*maxres
5083 !el if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
5084 ! Zeroing out fricmat
5090 ! Load the friction coefficients corresponding to peptide groups
5095 gamvec(ind1)=gamp(mnum)
5097 ! Load the friction coefficients corresponding to side chains
5101 gamsc(ntyp1_molec(j),j)=1.0d0
5108 gamvec(ii)=gamsc(iabs(iti),mnum)
5110 if (surfarea) call sdarea(gamvec)
5112 ! write (iout,*) "Matrix A and vector gamma"
5114 ! write (iout,'(i2,$)') i
5116 ! write (iout,'(f4.1,$)') A(i,j)
5118 ! write (iout,'(f8.3)') gamvec(i)
5122 write (iout,*) "Vector gamvec"
5124 write (iout,'(i5,f10.5)') i, gamvec(i)
5128 ! The friction matrix
5133 dtdi=dtdi+A(j,k)*A(j,i)*gamvec(j)
5140 write (iout,'(//a)') "Matrix fricmat"
5141 call matout2(dimen,dimen,nres2,nres2,fricmat)
5143 if (lang.eq.2 .or. lang.eq.3) then
5144 ! Mass-scale the friction matrix if non-direct integration will be performed
5150 Ginvfric(i,j)=Ginvfric(i,j)+ &
5151 Gsqrm(i,k)*Gsqrm(l,j)*fricmat(k,l)
5156 ! Diagonalize the friction matrix
5161 Ghalf(ind)=Ginvfric(i,j)
5164 call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
5167 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
5168 " mass-scaled friction matrix"
5169 call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
5171 ! Precompute matrices for tinker stochastic integrator
5178 mt1(i,j)=mt1(i,j)+fricvec(k,i)*gsqrm(k,j)
5179 mt2(i,j)=mt2(i,j)+fricvec(k,i)*gsqrp(k,j)
5185 else if (lang.eq.4) then
5186 ! Diagonalize the friction matrix
5191 Ghalf(ind)=fricmat(i,j)
5194 call gldiag(nres2,dimen,dimen,Ghalf,work,fricgam,fricvec,&
5197 write (iout,'(//2a)') "Eigenvectors and eigenvalues of the",&
5199 call eigout(dimen,dimen,nres2,nres2,fricvec,fricgam)
5201 ! Determine the number of zero eigenvalues of the friction matrix
5202 nzero=max0(dimen-dimen1,0)
5203 ! do while (fricgam(nzero+1).le.1.0d-5 .and. nzero.lt.dimen)
5206 write (iout,*) "Number of zero eigenvalues:",nzero
5211 fricmat(i,j)=fricmat(i,j) &
5212 +fricvec(i,k)*fricvec(j,k)/fricgam(k)
5217 write (iout,'(//a)') "Generalized inverse of fricmat"
5218 call matout(dimen,dimen,nres6,nres6,fricmat)
5223 if (nfgtasks.gt.1) then
5224 if (fg_rank.eq.0) then
5225 ! The matching BROADCAST for fg processors is called in ERGASTULUM
5231 call MPI_Bcast(10,1,MPI_INTEGER,king,FG_COMM,IERROR)
5233 time_Bcast=time_Bcast+MPI_Wtime()-time00
5235 time_Bcast=time_Bcast+tcpu()-time00
5237 ! print *,"Processor",myrank,
5238 ! & " BROADCAST iorder in SETUP_FRICMAT"
5241 write (iout,*) "setup_fricmat licznik"!,licznik !sp
5247 ! Scatter the friction matrix
5248 call MPI_Scatterv(fricmat(1,1),nginv_counts(0),&
5249 nginv_start(0),MPI_DOUBLE_PRECISION,fcopy(1,1),&
5250 myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
5253 time_scatter=time_scatter+MPI_Wtime()-time00
5254 time_scatter_fmat=time_scatter_fmat+MPI_Wtime()-time00
5256 time_scatter=time_scatter+tcpu()-time00
5257 time_scatter_fmat=time_scatter_fmat+tcpu()-time00
5261 do j=1,2*my_ng_count
5262 fricmat(j,i)=fcopy(i,j)
5265 ! write (iout,*) "My chunk of fricmat"
5266 ! call MATOUT2(my_ng_count,dimen,maxres2,maxres2,fcopy)
5270 end subroutine setup_fricmat
5271 !-----------------------------------------------------------------------------
5272 subroutine stochastic_force(stochforcvec)
5275 use random, only:anorm_distr
5276 ! implicit real*8 (a-h,o-z)
5277 ! include 'DIMENSIONS'
5278 use control, only: tcpu
5283 ! include 'COMMON.VAR'
5284 ! include 'COMMON.CHAIN'
5285 ! include 'COMMON.DERIV'
5286 ! include 'COMMON.GEO'
5287 ! include 'COMMON.LOCAL'
5288 ! include 'COMMON.INTERACT'
5289 ! include 'COMMON.MD'
5290 ! include 'COMMON.TIME1'
5292 ! include 'COMMON.LANGEVIN'
5294 ! include 'COMMON.LANGEVIN.lang0'
5296 ! include 'COMMON.IOUNITS'
5298 real(kind=8) :: x,sig,lowb,highb
5299 real(kind=8) :: ff(3),force(3,0:2*nres),zeta2,lowb2
5300 real(kind=8) :: highb2,sig2,forcvec(6*nres),stochforcvec(6*nres)
5301 real(kind=8) :: time00
5302 logical :: lprn = .false.
5303 integer :: i,j,ind,mnum
5307 stochforc(j,i)=0.0d0
5317 ! Compute the stochastic forces acting on bodies. Store in force.
5323 force(j,i)=anorm_distr(x,sig,lowb,highb)
5331 force(j,i+nres)=anorm_distr(x,sig2,lowb2,highb2)
5335 time_fsample=time_fsample+MPI_Wtime()-time00
5337 time_fsample=time_fsample+tcpu()-time00
5339 ! Compute the stochastic forces acting on virtual-bond vectors.
5345 stochforc(j,i)=ff(j)+0.5d0*force(j,i)
5348 ff(j)=ff(j)+force(j,i)
5350 ! if (itype(i+1,1).ne.ntyp1) then
5352 if (itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
5354 stochforc(j,i)=stochforc(j,i)+force(j,i+nres+1)
5355 ff(j)=ff(j)+force(j,i+nres+1)
5360 stochforc(j,0)=ff(j)+force(j,nnt+nres)
5364 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5365 .and.(mnum.ne.5)) then
5366 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
5368 stochforc(j,i+nres)=force(j,i+nres)
5374 stochforcvec(j)=stochforc(j,0)
5379 stochforcvec(ind+j)=stochforc(j,i)
5385 if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
5386 .and.(mnum.ne.5)) then
5387 ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
5389 stochforcvec(ind+j)=stochforc(j,i+nres)
5395 write (iout,*) "stochforcvec"
5397 write(iout,'(i5,e15.5)') i,stochforcvec(i)
5399 write(iout,*) "Stochastic forces backbone"
5401 write(iout,'(i5,3e15.5)') i,(stochforc(j,i),j=1,3)
5403 write(iout,*) "Stochastic forces side chain"
5405 write(iout,'(i5,3e15.5)') &
5406 i,(stochforc(j,i+nres),j=1,3)
5414 write (iout,*) i,ind
5416 forcvec(ind+j)=force(j,i)
5421 write (iout,*) i,ind
5423 forcvec(j+ind)=force(j,i+nres)
5428 write (iout,*) "forcvec"
5432 write (iout,'(2i3,2f10.5)') i,j,force(j,i),&
5439 write (iout,'(2i3,2f10.5)') i,j,force(j,i+nres),&
5448 end subroutine stochastic_force
5449 !-----------------------------------------------------------------------------
5450 subroutine sdarea(gamvec)
5452 ! Scale the friction coefficients according to solvent accessible surface areas
5453 ! Code adapted from TINKER
5457 ! implicit real*8 (a-h,o-z)
5458 ! include 'DIMENSIONS'
5459 ! include 'COMMON.CONTROL'
5460 ! include 'COMMON.VAR'
5461 ! include 'COMMON.MD'
5463 ! include 'COMMON.LANGEVIN'
5465 ! include 'COMMON.LANGEVIN.lang0'
5467 ! include 'COMMON.CHAIN'
5468 ! include 'COMMON.DERIV'
5469 ! include 'COMMON.GEO'
5470 ! include 'COMMON.LOCAL'
5471 ! include 'COMMON.INTERACT'
5472 ! include 'COMMON.IOUNITS'
5473 ! include 'COMMON.NAMES'
5474 real(kind=8),dimension(2*nres) :: radius,gamvec !(maxres2)
5475 real(kind=8),parameter :: twosix = 1.122462048309372981d0
5476 logical :: lprn = .false.
5477 real(kind=8) :: probe,area,ratio
5478 integer :: i,j,ind,iti,mnum
5480 ! determine new friction coefficients every few SD steps
5482 ! set the atomic radii to estimates of sigma values
5484 ! print *,"Entered sdarea"
5490 ! Load peptide group radii
5493 radius(i)=pstok(mnum)
5495 ! Load side chain radii
5499 radius(i+nres)=restok(iti,mnum)
5502 ! write (iout,*) "i",i," radius",radius(i)
5505 radius(i) = radius(i) / twosix
5506 if (radius(i) .ne. 0.0d0) radius(i) = radius(i) + probe
5509 ! scale atomic friction coefficients by accessible area
5511 if (lprn) write (iout,*) &
5512 "Original gammas, surface areas, scaling factors, new gammas, ",&
5513 "std's of stochastic forces"
5516 if (radius(i).gt.0.0d0) then
5517 call surfatom (i,area,radius)
5518 ratio = dmax1(area/(4.0d0*pi*radius(i)**2),1.0d-1)
5519 if (lprn) write (iout,'(i5,3f10.5,$)') &
5520 i,gamvec(ind+1),area,ratio
5523 gamvec(ind) = ratio * gamvec(ind)
5526 stdforcp(i)=stdfp(mnum)*dsqrt(gamvec(ind))
5527 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcp(i)
5531 if (radius(i+nres).gt.0.0d0) then
5532 call surfatom (i+nres,area,radius)
5533 ratio = dmax1(area/(4.0d0*pi*radius(i+nres)**2),1.0d-1)
5534 if (lprn) write (iout,'(i5,3f10.5,$)') &
5535 i,gamvec(ind+1),area,ratio
5538 gamvec(ind) = ratio * gamvec(ind)
5541 stdforcsc(i)=stdfsc(itype(i,mnum),mnum)*dsqrt(gamvec(ind))
5542 if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcsc(i)
5547 end subroutine sdarea
5548 !-----------------------------------------------------------------------------
5550 !-----------------------------------------------------------------------------
5553 ! ###################################################
5554 ! ## COPYRIGHT (C) 1996 by Jay William Ponder ##
5555 ! ## All Rights Reserved ##
5556 ! ###################################################
5558 ! ################################################################
5560 ! ## subroutine surfatom -- exposed surface area of an atom ##
5562 ! ################################################################
5565 ! "surfatom" performs an analytical computation of the surface
5566 ! area of a specified atom; a simplified version of "surface"
5568 ! literature references:
5570 ! T. J. Richmond, "Solvent Accessible Surface Area and
5571 ! Excluded Volume in Proteins", Journal of Molecular Biology,
5574 ! L. Wesson and D. Eisenberg, "Atomic Solvation Parameters
5575 ! Applied to Molecular Dynamics of Proteins in Solution",
5576 ! Protein Science, 1, 227-235 (1992)
5578 ! variables and parameters:
5580 ! ir number of atom for which area is desired
5581 ! area accessible surface area of the atom
5582 ! radius radii of each of the individual atoms
5585 subroutine surfatom(ir,area,radius)
5587 ! implicit real*8 (a-h,o-z)
5588 ! include 'DIMENSIONS'
5590 ! include 'COMMON.GEO'
5591 ! include 'COMMON.IOUNITS'
5593 integer :: nsup,nstart_sup
5594 ! double precision c,dc,dc_old,d_c_work,xloc,xrot,dc_norm
5595 ! common /chain/ c(3,maxres2+2),dc(3,0:maxres2),dc_old(3,0:maxres2),
5596 ! & xloc(3,maxres),xrot(3,maxres),dc_norm(3,0:maxres2),
5597 ! & dc_work(MAXRES6),nres,nres0
5598 integer,parameter :: maxarc=300
5602 integer :: mi,ni,narc
5603 integer :: key(maxarc)
5604 integer :: intag(maxarc)
5605 integer :: intag1(maxarc)
5606 real(kind=8) :: area,arcsum
5607 real(kind=8) :: arclen,exang
5608 real(kind=8) :: delta,delta2
5609 real(kind=8) :: eps,rmove
5610 real(kind=8) :: xr,yr,zr
5611 real(kind=8) :: rr,rrsq
5612 real(kind=8) :: rplus,rminus
5613 real(kind=8) :: axx,axy,axz
5614 real(kind=8) :: ayx,ayy
5615 real(kind=8) :: azx,azy,azz
5616 real(kind=8) :: uxj,uyj,uzj
5617 real(kind=8) :: tx,ty,tz
5618 real(kind=8) :: txb,tyb,td
5619 real(kind=8) :: tr2,tr,txr,tyr
5620 real(kind=8) :: tk1,tk2
5621 real(kind=8) :: thec,the,t,tb
5622 real(kind=8) :: txk,tyk,tzk
5623 real(kind=8) :: t1,ti,tf,tt
5624 real(kind=8) :: txj,tyj,tzj
5625 real(kind=8) :: ccsq,cc,xysq
5626 real(kind=8) :: bsqk,bk,cosine
5627 real(kind=8) :: dsqj,gi,pix2
5628 real(kind=8) :: therk,dk,gk
5629 real(kind=8) :: risqk,rik
5630 real(kind=8) :: radius(2*nres) !(maxatm) (maxatm=maxres2)
5631 real(kind=8) :: ri(maxarc),risq(maxarc)
5632 real(kind=8) :: ux(maxarc),uy(maxarc),uz(maxarc)
5633 real(kind=8) :: xc(maxarc),yc(maxarc),zc(maxarc)
5634 real(kind=8) :: xc1(maxarc),yc1(maxarc),zc1(maxarc)
5635 real(kind=8) :: dsq(maxarc),bsq(maxarc)
5636 real(kind=8) :: dsq1(maxarc),bsq1(maxarc)
5637 real(kind=8) :: arci(maxarc),arcf(maxarc)
5638 real(kind=8) :: ex(maxarc),lt(maxarc),gr(maxarc)
5639 real(kind=8) :: b(maxarc),b1(maxarc),bg(maxarc)
5640 real(kind=8) :: kent(maxarc),kout(maxarc)
5641 real(kind=8) :: ther(maxarc)
5642 logical :: moved,top
5643 logical :: omit(maxarc)
5646 ! maxatm = 2*nres !maxres2 maxres2=2*maxres
5647 ! maxlight = 8*maxatm
5650 ! maxtors = 4*maxatm
5652 ! zero out the surface area for the sphere of interest
5655 ! write (2,*) "ir",ir," radius",radius(ir)
5656 if (radius(ir) .eq. 0.0d0) return
5658 ! set the overlap significance and connectivity shift
5662 delta2 = delta * delta
5667 ! store coordinates and radius of the sphere of interest
5675 ! initialize values of some counters and summations
5684 ! test each sphere to see if it overlaps the sphere of interest
5687 if (i.eq.ir .or. radius(i).eq.0.0d0) goto 30
5688 rplus = rr + radius(i)
5690 if (abs(tx) .ge. rplus) goto 30
5692 if (abs(ty) .ge. rplus) goto 30
5694 if (abs(tz) .ge. rplus) goto 30
5696 ! check for sphere overlap by testing distance against radii
5698 xysq = tx*tx + ty*ty
5699 if (xysq .lt. delta2) then
5706 if (rplus-cc .le. delta) goto 30
5707 rminus = rr - radius(i)
5709 ! check to see if sphere of interest is completely buried
5711 if (cc-abs(rminus) .le. delta) then
5712 if (rminus .le. 0.0d0) goto 170
5716 ! check for too many overlaps with sphere of interest
5718 if (io .ge. maxarc) then
5720 20 format (/,' SURFATOM -- Increase the Value of MAXARC')
5724 ! get overlap between current sphere and sphere of interest
5733 gr(io) = (ccsq+rplus*rminus) / (2.0d0*rr*b1(io))
5739 ! case where no other spheres overlap the sphere of interest
5742 area = 4.0d0 * pi * rrsq
5746 ! case where only one sphere overlaps the sphere of interest
5749 area = pix2 * (1.0d0 + gr(1))
5750 area = mod(area,4.0d0*pi) * rrsq
5754 ! case where many spheres intersect the sphere of interest;
5755 ! sort the intersecting spheres by their degree of overlap
5757 call sort2 (io,gr,key)
5760 intag(i) = intag1(k)
5769 ! get radius of each overlap circle on surface of the sphere
5774 risq(i) = rrsq - gi*gi
5775 ri(i) = sqrt(risq(i))
5776 ther(i) = 0.5d0*pi - asin(min(1.0d0,max(-1.0d0,gr(i))))
5779 ! find boundary of inaccessible area on sphere of interest
5782 if (.not. omit(k)) then
5789 ! check to see if J circle is intersecting K circle;
5790 ! get distance between circle centers and sum of radii
5793 if (omit(j)) goto 60
5794 cc = (txk*xc(j)+tyk*yc(j)+tzk*zc(j))/(bk*b(j))
5795 cc = acos(min(1.0d0,max(-1.0d0,cc)))
5796 td = therk + ther(j)
5798 ! check to see if circles enclose separate regions
5800 if (cc .ge. td) goto 60
5802 ! check for circle J completely inside circle K
5804 if (cc+ther(j) .lt. therk) goto 40
5806 ! check for circles that are essentially parallel
5808 if (cc .gt. delta) goto 50
5813 ! check to see if sphere of interest is completely buried
5816 if (pix2-cc .le. td) goto 170
5822 ! find T value of circle intersections
5825 if (omit(k)) goto 110
5840 ! rotation matrix elements
5852 if (.not. omit(j)) then
5857 ! rotate spheres so K vector colinear with z-axis
5859 uxj = txj*axx + tyj*axy - tzj*axz
5860 uyj = tyj*ayy - txj*ayx
5861 uzj = txj*azx + tyj*azy + tzj*azz
5862 cosine = min(1.0d0,max(-1.0d0,uzj/b(j)))
5863 if (acos(cosine) .lt. therk+ther(j)) then
5864 dsqj = uxj*uxj + uyj*uyj
5869 tr2 = risqk*dsqj - tb*tb
5875 ! get T values of intersection for K circle
5878 tb = min(1.0d0,max(-1.0d0,tb))
5880 if (tyb-txr .lt. 0.0d0) tk1 = pix2 - tk1
5882 tb = min(1.0d0,max(-1.0d0,tb))
5884 if (tyb+txr .lt. 0.0d0) tk2 = pix2 - tk2
5885 thec = (rrsq*uzj-gk*bg(j)) / (rik*ri(j)*b(j))
5886 if (abs(thec) .lt. 1.0d0) then
5888 else if (thec .ge. 1.0d0) then
5890 else if (thec .le. -1.0d0) then
5894 ! see if "tk1" is entry or exit point; check t=0 point;
5895 ! "ti" is exit point, "tf" is entry point
5897 cosine = min(1.0d0,max(-1.0d0, &
5898 (uzj*gk-uxj*rik)/(b(j)*rr)))
5899 if ((acos(cosine)-ther(j))*(tk2-tk1) .le. 0.0d0) then
5907 if (narc .ge. maxarc) then
5909 70 format (/,' SURFATOM -- Increase the Value',&
5913 if (tf .le. ti) then
5934 ! special case; K circle without intersections
5936 if (narc .le. 0) goto 90
5938 ! general case; sum up arclength and set connectivity code
5940 call sort2 (narc,arci,key)
5945 if (narc .gt. 1) then
5948 if (t .lt. arci(j)) then
5949 arcsum = arcsum + arci(j) - t
5950 exang = exang + ex(ni)
5952 if (jb .ge. maxarc) then
5954 80 format (/,' SURFATOM -- Increase the Value',&
5959 kent(jb) = maxarc*i + k
5961 kout(jb) = maxarc*k + i
5970 arcsum = arcsum + pix2 - t
5972 exang = exang + ex(ni)
5975 kent(jb) = maxarc*i + k
5977 kout(jb) = maxarc*k + i
5984 arclen = arclen + gr(k)*arcsum
5987 if (arclen .eq. 0.0d0) goto 170
5988 if (jb .eq. 0) goto 150
5990 ! find number of independent boundaries and check connectivity
5994 if (kout(k) .ne. 0) then
6001 if (m .eq. kent(ii)) then
6004 if (j .eq. jb) goto 150
6016 ! attempt to fix connectivity error by moving atom slightly
6020 140 format (/,' SURFATOM -- Connectivity Error at Atom',i6)
6029 ! compute the exposed surface area for the sphere of interest
6032 area = ib*pix2 + exang + arclen
6033 area = mod(area,4.0d0*pi) * rrsq
6035 ! attempt to fix negative area by moving atom slightly
6037 if (area .lt. 0.0d0) then
6040 160 format (/,' SURFATOM -- Negative Area at Atom',i6)
6051 end subroutine surfatom
6052 !----------------------------------------------------------------
6053 !----------------------------------------------------------------
6054 subroutine alloc_MD_arrays
6055 !EL Allocation of arrays used by MD module
6057 integer :: nres2,nres6
6060 !----------------------
6064 allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
6065 allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
6066 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
6067 allocate(fricvec(nres2,nres2))
6068 allocate(pfric_mat(nres2,nres2),vfric_mat(nres2,nres2))
6069 allocate(afric_mat(nres2,nres2),prand_mat(nres2,nres2))
6070 allocate(vrand_mat1(nres2,nres2),vrand_mat2(nres2,nres2)) !(MAXRES2,MAXRES2)
6071 allocate(pfric0_mat(nres2,nres2,0:maxflag_stoch))
6072 allocate(afric0_mat(nres2,nres2,0:maxflag_stoch))
6073 allocate(vfric0_mat(nres2,nres2,0:maxflag_stoch))
6074 allocate(prand0_mat(nres2,nres2,0:maxflag_stoch))
6075 allocate(vrand0_mat1(nres2,nres2,0:maxflag_stoch))
6076 allocate(vrand0_mat2(nres2,nres2,0:maxflag_stoch)) !(MAXRES2,MAXRES2,0:maxflag_stoch)
6077 allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
6079 allocate(mt1(nres2,nres2),mt2(nres2,nres2),mt3(nres2,nres2)) !(maxres2,maxres2)
6080 !----------------------
6082 ! commom.langevin.lang0
6084 allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
6085 if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
6086 allocate(fricvec(nres2,nres2)) !(MAXRES2,MAXRES2)
6087 allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
6088 allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
6091 !el if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2))
6092 !----------------------
6093 ! commom.hairpin in CSA module
6094 !----------------------
6095 ! common.mce in MCM_MD module
6096 !----------------------
6098 ! common /mdgrad/ in module.energy
6099 ! common /back_constr/ in module.energy
6100 ! common /qmeas/ in module.energy
6103 allocate(potEcomp(0:n_ene+4)) !(0:n_ene+4)
6105 allocate(d_t(3,0:nres2),d_a(3,0:nres2),d_t_old(3,0:nres2)) !(3,0:MAXRES2)
6106 allocate(d_a_work(nres6)) !(6*MAXRES)
6108 allocate(DM(nres2),DU1(nres2),DU2(nres2))
6109 allocate(DMorig(nres2),DU1orig(nres2),DU2orig(nres2))
6111 allocate(Gmat(nres2,nres2),A(nres2,nres2))
6112 if(.not.allocated(Ginv)) allocate(Ginv(nres2,nres2)) !in control: ergastulum
6113 allocate(Gsqrp(nres2,nres2),Gsqrm(nres2,nres2),Gvec(nres2,nres2)) !(maxres2,maxres2)
6115 allocate(Geigen(nres2)) !(maxres2)
6116 if(.not.allocated(vtot)) allocate(vtot(nres2)) !(maxres2)
6117 ! common /inertia/ in io_conf: parmread
6118 ! real(kind=8),dimension(:),allocatable :: ISC,msc !(ntyp+1)
6119 ! common /langevin/in io read_MDpar
6120 ! real(kind=8),dimension(:),allocatable :: gamsc !(ntyp1)
6121 ! real(kind=8),dimension(:),allocatable :: stdfsc !(ntyp)
6122 ! in io_conf: parmread
6123 ! real(kind=8),dimension(:),allocatable :: restok !(ntyp+1)
6124 ! common /mdpmpi/ in control: ergastulum
6125 if(.not.allocated(ng_start)) allocate(ng_start(0:nfgtasks-1))
6126 if(.not.allocated(ng_counts)) allocate(ng_counts(0:nfgtasks-1))
6127 if(.not.allocated(nginv_counts)) allocate(nginv_counts(0:nfgtasks-1)) !(0:MaxProcs-1)
6128 if(.not.allocated(nginv_start)) allocate(nginv_start(0:nfgtasks)) !(0:MaxProcs)
6129 !----------------------
6130 ! common.muca in read_muca
6131 ! common /double_muca/
6132 ! real(kind=8) :: elow,ehigh,factor,hbin,factor_min
6133 ! real(kind=8),dimension(:),allocatable :: emuca,nemuca,&
6134 ! nemuca2,hist !(4*maxres)
6135 ! common /integer_muca/
6136 ! integer :: nmuca,imtime,muca_smooth
6138 ! real(kind=8),dimension(:),allocatable :: elowi,ehighi !(maxprocs)
6139 !----------------------
6141 ! common /mdgrad/ in module.energy
6142 ! common /back_constr/ in module.energy
6143 ! common /qmeas/ in module.energy
6147 allocate(d_t_work(nres6),d_t_work_new(nres6),d_af_work(nres6))
6148 allocate(d_as_work(nres6),kinetic_force(nres6)) !(MAXRES6)
6149 allocate(d_t_new(3,0:nres2),d_a_old(3,0:nres2),d_a_short(3,0:nres2)) !,d_a !(3,0:MAXRES2)
6150 allocate(stdforcp(nres),stdforcsc(nres)) !(MAXRES)
6151 !----------------------
6153 allocate(D_ban(nres6)) !(MAXRES6) maxres6=6*maxres
6154 ! common /stochcalc/ stochforcvec
6155 allocate(stochforcvec(nres6)) !(MAXRES6) maxres6=6*maxres
6156 !----------------------
6158 end subroutine alloc_MD_arrays
6159 !-----------------------------------------------------------------------------
6160 !-----------------------------------------------------------------------------