! Calculate energy and forces
call zerograd
call etotal(potEcomp)
- potE=potEcomp(0)-potEcomp(20)
+ potE=potEcomp(0)-potEcomp(51)
call cartgrad
totT=totT+d_time
totTafm=totT
end subroutine gauss
!-----------------------------------------------------------------------------
! kinetic_lesyng.f
+#ifdef FIVEDIAG
+ subroutine kinetic(KE_total)
+!c----------------------------------------------------------------
+!c This subroutine calculates the total kinetic energy of the chain
+!c-----------------------------------------------------------------
+!c 3/5/2020 AL Corrected for multichain systems, no fake peptide groups
+!c inside, implemented with five-diagonal inertia matrix
+ use energy_data
+ implicit none
+ real(kind=8):: KE_total,KEt_p,KEt_sc,KEr_p,KEr_sc,mag1,mag2
+ integer i,j,k,iti,mnum
+ real(kind=8),dimension(3) :: incr,v
+
+ KEt_p=0.0d0
+ KEt_sc=0.0d0
+ KEr_p=0.0D0
+ KEr_sc=0.0D0
+!c write (iout,*) "ISC",(isc(itype(i)),i=1,nres)
+!c The translational part for peptide virtual bonds
+ do j=1,3
+ incr(j)=d_t(j,0)
+ enddo
+ do i=nnt,nct-1 !czy na pewno nct-1??
+ mnum=molnum(i)
+!c write (iout,*) "Kinetic trp:",i,(incr(j),j=1,3
+!c Skip dummy peptide groups
+ if (itype(i,mnum).ne.ntyp1_molec(mnum)&
+ .and. itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
+ do j=1,3
+ v(j)=incr(j)+0.5d0*d_t(j,i)
+ enddo
+ if (mnum.eq.5) mp(mnum)=0.0d0
+! if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
+!c write (iout,*) "Kinetic trp:",i,(v(j),j=1,3)
+ vtot(i)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
+ KEt_p=KEt_p+mp(mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
+ endif
+ do j=1,3
+ incr(j)=incr(j)+d_t(j,i)
+ enddo
+ enddo
+!c write(iout,*) 'KEt_p', KEt_p
+!c The translational part for the side chain virtual bond
+!c Only now we can initialize incr with zeros. It must be equal
+!c to the velocities of the first Calpha.
+ do j=1,3
+ incr(j)=d_t(j,0)
+ enddo
+ do i=nnt,nct
+ mnum=molnum(i)
+ iti=iabs(itype(i,mnum))
+ if (mnum.eq.5) iti=itype(i,mnum)
+ if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
+ .or.mnum.ge.3) then
+ do j=1,3
+ v(j)=incr(j)
+ enddo
+ else
+ do j=1,3
+ v(j)=incr(j)+d_t(j,nres+i)
+ enddo
+ endif
+! if (mnum.ne.5) then
+! write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3)
+! write (iout,*) "i",i," msc",msc(iti,mnum)," v",(v(j),j=1,3)
+ KEt_sc=KEt_sc+msc(iti,mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
+ vtot(i+nres)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
+! endif
+ do j=1,3
+ incr(j)=incr(j)+d_t(j,i)
+ enddo
+ enddo
+! goto 111
+! write(iout,*) 'KEt_sc', KEt_sc
+! The part due to stretching and rotation of the peptide groups
+ do i=nnt,nct-1
+ mnum=molnum(i)
+ if (itype(i,mnum).ne.ntyp1_molec(mnum)&
+ .and.itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
+ if (mnum.eq.5) Ip(mnum)=Isc(itype(i,mnum),mnum)
+! write (iout,*) "i",i
+! write (iout,*) "i",i," mag1",mag1," mag2",mag2
+ do j=1,3
+ incr(j)=d_t(j,i)
+ enddo
+!c write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3)
+ KEr_p=KEr_p+Ip(mnum)*(incr(1)*incr(1)+incr(2)*incr(2) &
+ +incr(3)*incr(3))
+ endif
+ enddo
+!c goto 111
+!c write(iout,*) 'KEr_p', KEr_p
+!c The rotational part of the side chain virtual bond
+ do i=nnt,nct
+ mnum=molnum(i)
+ iti=iabs(itype(i,mnum))
+ if (itype(i,1).ne.10.and.itype(i,mnum).ne.ntyp1_molec(mnum)&
+ .and.mnum.lt.3) then
+ do j=1,3
+ incr(j)=d_t(j,nres+i)
+ enddo
+! write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3)
+ KEr_sc=KEr_sc+Isc(iti,mnum)*(incr(1)*incr(1)+incr(2)*incr(2)+&
+ incr(3)*incr(3))
+ endif
+ enddo
+!c The total kinetic energy
+ 111 continue
+! write(iout,*) ' KEt_p',KEt_p,' KEt_sc',KEt_sc,' KEr_p',KEr_p, &
+! ' KEr_sc', KEr_sc
+ KE_total=0.5d0*(KEt_p+KEt_sc+0.25d0*KEr_p+KEr_sc)
+!c write (iout,*) "KE_total",KE_total
+ return
+ end subroutine kinetic
+#else
+
!-----------------------------------------------------------------------------
subroutine kinetic(KE_total)
!----------------------------------------------------------------
! include 'COMMON.IOUNITS'
real(kind=8) :: KE_total,mscab
- integer :: i,j,k,iti,mnum
+ integer :: i,j,k,iti,mnum,term
real(kind=8) :: KEt_p,KEt_sc,KEr_p,KEr_sc,incr(3),&
mag1,mag2,v(3)
#ifdef DEBUG
do j=1,3
incr(j)=d_t(j,0)
enddo
- do i=nnt,nct-1
+ term=nct-1
+! if (molnum(nct).gt.3) term=nct
+ do i=nnt,term
mnum=molnum(i)
- if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
-! write (iout,*) "Kinetic trp:",i,(incr(j),j=1,3)
+ if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
+! write (iout,*) "Kinetic trp:",i,(incr(j),j=1,3),mp(mnum)
+ if (mnum.gt.4) then
do j=1,3
v(j)=incr(j)+0.5d0*d_t(j,i)
- enddo
+ enddo
+ else
+ do j=1,3
+ v(j)=incr(j)+0.5d0*d_t(j,i)
+ enddo
+ endif
vtot(i)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
KEt_p=KEt_p+mp(mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
do j=1,3
do i=nnt,nct
mnum=molnum(i)
iti=iabs(itype(i,mnum))
- if (mnum.eq.5) then
- mscab=0.0d0
- else
+! if (mnum.ge.4) then
+! mscab=0.0d0
+! else
mscab=msc(iti,mnum)
- endif
+! endif
! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
- .or.(mnum.eq.5)) then
+ .or.(mnum.ge.4)) then
do j=1,3
v(j)=incr(j)
enddo
enddo
endif
! write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3)
-! write (iout,*) "i",i," msc",msc(iti)," v",(v(j),j=1,3)
+! write (iout,*) "i",i," msc",msc(iti,mnum)," v",(v(j),j=1,3)
KEt_sc=KEt_sc+mscab*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
vtot(i+nres)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
do j=1,3
iti=iabs(itype(i,mnum))
! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
- .and.(mnum.ne.5)) then
+ .and.(mnum.lt.4)) then
do j=1,3
incr(j)=d_t(j,nres+i)
enddo
return
end subroutine kinetic
!-----------------------------------------------------------------------------
+#endif
+ subroutine kinetic_CASC(KE_total)
+!c----------------------------------------------------------------
+!c Compute the kinetic energy of the system using the Calpha-SC
+!c coordinate system
+!c-----------------------------------------------------------------
+ implicit none
+ double precision KE_total
+
+ integer i,j,k,iti,ichain,innt,inct,mnum
+ double precision KEt_p,KEt_sc,KEr_p,KEr_sc,incr(3),&
+ mag1,mag2,v(3)
+#ifdef FIVEDIAG
+ KEt_p=0.0d0
+ KEt_sc=0.0d0
+ KEr_p=0.0D0
+ KEr_sc=0.0D0
+!c write (iout,*) "ISC",(isc(itype(i)),i=1,nres)
+!c The translational part for peptide virtual bonds
+ do ichain=1,nchain
+
+ innt=chain_border(1,ichain)
+ inct=chain_border(2,ichain)
+!c write (iout,*) "Kinetic_CASC chain",ichain," innt",innt,
+!c & " inct",inct
+
+ do i=innt,inct-1
+ mnum=molnum(i)
+ if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
+ if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
+!c write (iout,*) i,(d_t(j,i),j=1,3),(d_t(j,i+1),j=1,3)
+ do j=1,3
+ v(j)=0.5d0*(d_t(j,i)+d_t(j,i+1))
+ enddo
+!c write (iout,*) "Kinetic trp i",i," v",(v(j),j=1,3)
+ KEt_p=KEt_p+mp(mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
+ enddo
+!c write(iout,*) 'KEt_p', KEt_p
+!c The translational part for the side chain virtual bond
+!c Only now we can initialize incr with zeros. It must be equal
+!c to the velocities of the first Calpha.
+ do i=innt,inct
+ mnum=molnum(i)
+ iti=iabs(itype(i,mnum))
+ if (itype(i,1).eq.10.or.mnum.ge.3.or. itype(i,mnum).eq.ntyp1_molec(mnum)) then
+!c write (iout,*) i,iti,(d_t(j,i),j=1,3)
+ do j=1,3
+ v(j)=d_t(j,i)
+ enddo
+ else
+!c write (iout,*) i,iti,(d_t(j,nres+i),j=1,3)
+ do j=1,3
+ v(j)=d_t(j,nres+i)
+ enddo
+ endif
+!c write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3)
+!c write (iout,*) "i",i," msc",msc(iti)," v",(v(j),j=1,3)
+ KEt_sc=KEt_sc+msc(iti,mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
+ enddo
+!c goto 111
+!c write(iout,*) 'KEt_sc', KEt_sc
+!c The part due to stretching and rotation of the peptide groups
+ do i=innt,inct-1
+ mnum=molnum(i)
+ do j=1,3
+ incr(j)=d_t(j,i+1)-d_t(j,i)
+ enddo
+ if (mnum.eq.5) Ip(mnum)=Isc(itype(i,mnum),mnum)
+!c write (iout,*) i,(incr(j),j=1,3)
+!c write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3)
+ KEr_p=KEr_p+Ip(mnum)*(incr(1)*incr(1)+incr(2)*incr(2)&
+ +incr(3)*incr(3))
+ enddo
+!c goto 111
+!c write(iout,*) 'KEr_p', KEr_p
+!c The rotational part of the side chain virtual bond
+ do i=innt,inct
+ mnum=molnum(i)
+ iti=iabs(itype(i,mnum))
+! if (iti.ne.10.and.mnum.lt.3) then
+ if (itype(i,1).ne.10.and.mnum.lt.3.and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
+ do j=1,3
+ incr(j)=d_t(j,nres+i)-d_t(j,i)
+ enddo
+!c write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3)
+ KEr_sc=KEr_sc+Isc(iti,mnum)*(incr(1)*incr(1)+incr(2)*incr(2)+&
+ incr(3)*incr(3))
+ endif
+ enddo
+
+ enddo ! ichain
+!c The total kinetic energy
+ 111 continue
+!c write(iout,*) ' KEt_p',KEt_p,' KEt_sc',KEt_sc,' KEr_p',KEr_p,
+!c & ' KEr_sc', KEr_sc
+ KE_total=0.5d0*(KEt_p+KEt_sc+0.25d0*KEr_p+KEr_sc)
+!c write (iout,*) "KE_total",KE_tota
+#else
+ write (iout,*) "Need to compile with -DFIVEDIAG to use this sub!"
+ stop
+#endif
+ return
+ end subroutine kinetic_CASC
+
! MD_A-MTS.F
!-----------------------------------------------------------------------------
subroutine MD
character(len=50) :: tytul
!el common /gucio/ cm
integer :: i,j,nharp
- integer,dimension(4,nres/3) :: iharp !(4,nres/3)(4,maxres/3)
+ integer,dimension(4,nres) :: iharp !(4,nres/3)(4,maxres/3)
! logical :: ovrtim
real(kind=8) :: tt0,scalfac
integer :: nres2,itime
tt0=tcpu()
#endif
do itime=1,n_timestep
+ if (large) print *,itime,ntwe
if (ovrtim()) exit
if (large.and. mod(itime,ntwe).eq.0) &
write (iout,*) "itime",itime
enddo
#endif
else if (lang.eq.1 .or. lang.eq.4) then
+ print *,"before setup_fricmat"
call setup_fricmat
+ print *,"after setup_fricmat"
endif
write (iout,'(a,i10)') &
"Friction matrix reset based on surface area, itime",itime
call RESPA_step(itime)
else
! Variable time step algorithm.
+ if (large) print *,"before verlet_step"
call velverlet_step(itime)
+ if (large) print *,"after verlet_step"
endif
else
#ifdef BROWN
itime_mat=itime
if (ntwe.ne.0) then
if (mod(itime,ntwe).eq.0) then
+! call returnbox
call statout(itime)
- call returnbox
+! call returnbox
! call check_ecartint
endif
#ifdef VOUT
enddo
do i=nnt,nct
mnum=molnum(i)
- if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1.and.mnum.ne.5) then
+ if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum).and.mnum.lt.4) then
do j=1,3
ind=ind+1
v_work(ind)=d_t(j,i+nres)
endif
if (mod(itime,ntwx).eq.0) then
call returnbox
+ call enerprint(potEcomp)
+
write (tytul,'("time",f8.2)') totT
if(mdpdb) then
call hairpin(.true.,nharp,iharp)
call secondary2(.true.)
call pdbout(potE,tytul,ipdb)
+! call enerprint(potEcomp)
else
call cartout(totT)
endif
+ if (fodson) then
+ write(iout,*) "starting fodstep"
+ call fodstep(nfodstep)
+ write(iout,*) "after fodstep"
+ call statout(itime)
+ if(mdpdb) then
+ call hairpin(.true.,nharp,iharp)
+ call secondary2(.true.)
+ call pdbout(potE,tytul,ipdb)
+ else
+ call cartout(totT)
+ endif
+ endif
+
endif
if (rstcount.eq.1000.or.itime.eq.n_timestep) then
open(irest2,file=rest2name,status='unknown')
icount_scale=0
if (lang.eq.1) then
call sddir_precalc
+ if (large) print *,"after sddir_precalc"
else if (lang.eq.2 .or. lang.eq.3) then
#ifndef LANG0
call stochastic_force(stochforcvec)
call chainbuild_cart
if (rattle) call rattle1
if (ntwe.ne.0) then
- if (large.and. mod(itime,ntwe).eq.0) then
+ if (large) then !.and. mod(itime,ntwe).eq.0) then
write (iout,*) "Cartesian and internal coordinates: step 1"
call cartprint
call intout
t_etotal=t_etotal+tcpu()-tt0
#endif
#endif
- potE=potEcomp(0)-potEcomp(20)
+ potE=potEcomp(0)-potEcomp(51)
call cartgrad
! Get the new accelerations
call lagrangian
t_elong=t_elong+tcpu()-tt0
#endif
#endif
+ potE=potEcomp(0)-potEcomp(51)
call cartgrad
call lagrangian
#ifdef MPI
do i=0,n_ene
potEcomp(i)=energia_short(i)+energia_long(i)
enddo
- potE=potEcomp(0)-potEcomp(20)
+ potE=potEcomp(0)-potEcomp(51)
! potE=energia_short(0)+energia_long(0)
totT=totT+d_time
totTafm=totT
mnum=molnum(i)
! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
- .and.(mnum.ne.5)) then
+ .and.(mnum.lt.4)) then
inres=i+nres
do j=1,3
d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
mnum=molnum(i)
! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
- .and.(mnum.ne.5)) then
+ .and.(mnum.lt.4)) then
inres=i+nres
do j=1,3
adt=d_a_old(j,inres)*d_time
! iti=iabs(itype(i,mnum))
! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
- .and.(mnum.ne.5)) then
+ .and.(mnum.lt.4)) then
inres=i+nres
do j=1,3
d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
!el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres
!el common /stochcalc/ stochforcvec
real(kind=8) :: time00
+ integer :: i
!
! Compute friction and stochastic forces
!
#ifdef MPI
time00=MPI_Wtime()
+ if (large) print *,"before friction_force"
call friction_force
+ if (large) print *,"after friction_force"
time_fric=time_fric+MPI_Wtime()-time00
time00=MPI_Wtime()
call stochastic_force(stochforcvec)
! Compute the acceleration due to friction forces (d_af_work) and stochastic
! forces (d_as_work)
!
+! call ginv_mult(fric_work, d_af_work)
+! call ginv_mult(stochforcvec, d_as_work)
+#ifdef FIVEDIAG
+ call fivediaginv_mult(dimen,fric_work, d_af_work)
+ call fivediaginv_mult(dimen,stochforcvec, d_as_work)
+ if (large) then
+ write(iout,*),"dimen",dimen
+ do i=1,dimen
+ write(iout,*) "fricwork",fric_work(i), d_af_work(i)
+ write(iout,*) "stochforcevec", stochforcvec(i), d_as_work(i)
+ enddo
+ endif
+#else
call ginv_mult(fric_work, d_af_work)
call ginv_mult(stochforcvec, d_as_work)
+#endif
+
return
end subroutine sddir_precalc
!-----------------------------------------------------------------------------
do j=1,3
adt=(d_a_old(j,0)+d_af_work(j))*d_time
adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
+! write(iout,*) i,"adt",adt,"ads",adt2,d_a_old(j,0),d_af_work(j),d_time
dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
d_t(j,0)=d_t_old(j,0)+adt
do j=1,3
adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
+! write(iout,*) i,"adt",adt,"ads",adt2,d_a_old(j,i),d_af_work(ind+j)
dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
d_t(j,i)=d_t_old(j,i)+adt
! iti=iabs(itype(i,mnum))
! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
- .and.(mnum.ne.5)) then
+ .and.(mnum.lt.4)) then
inres=i+nres
do j=1,3
adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
+! write(iout,*) i,"adt",adt,"ads",adt2,d_a_old(j,inres),d_af_work(ind+j)
dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
d_t(j,inres)=d_t_old(j,inres)+adt
ind=ind+3
endif
enddo
+
return
end subroutine sddir_verlet1
!-----------------------------------------------------------------------------
! Compute the acceleration due to friction forces (d_af_work) and stochastic
! forces (d_as_work)
!
+#ifdef FIVEDIAG
+ call fivediaginv_mult(6*nres,stochforcvec, d_as_work1)
+#else
call ginv_mult(stochforcvec, d_as_work1)
+#endif
!
! Update velocities
! iti=iabs(itype(i,mnum))
! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
- .and.(mnum.ne.5)) then
+ .and.(mnum.lt.4)) then
inres=i+nres
do j=1,3
d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres) &
! iti=iabs(itype(i,mnum))
! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
- .and.(mnum.ne.5)) then
+ .and.(mnum.lt.4)) then
do j=1,3
! accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
endif
! Side chains
if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1.and.&
- molnum(i).ne.5) then
+ molnum(i).lt.4) then
do j=1,3
epdriftij= &
dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
! iti=iabs(itype(i,mnum))
! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
- .and.(mnum.ne.5)) then
+ .and.(mnum.lt.4)) then
inres=i+nres
do j=1,3
d_t(j,inres)=fact*d_t(j,inres)
use minimm, only:minim_dc,minimize,sc_move
use io_config, only:readrst
use io, only:statout
+ use random, only: iran_num
! implicit real*8 (a-h,o-z)
! include 'DIMENSIONS'
#ifdef MP
character(len=16) :: form
integer :: IERROR,ERRCODE
#endif
- integer :: iranmin,itrial,itmp
+ integer :: iranmin,itrial,itmp,n_model_try,k, &
+ i_model
+ integer, dimension(:),allocatable :: list_model_try
+ integer, dimension(0:nodes-1) :: i_start_models
! include 'COMMON.SETUP'
! include 'COMMON.CONTROL'
! include 'COMMON.VAR'
! include 'COMMON.IOUNITS'
! include 'COMMON.NAMES'
! include 'COMMON.REMD'
- real(kind=8),dimension(0:n_ene) :: energia_long,energia_short
+ real(kind=8),dimension(0:n_ene) :: energia_long,energia_short,energia
real(kind=8),dimension(3) :: vcm,incr,L
real(kind=8) :: xv,sigv,lowb,highb
real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
if(me.eq.king.or..not.out1file)then
write (iout,*) "Initial velocities"
do i=0,nres
- write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
+ write (iout,'(i6,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
(d_t(j,i+nres),j=1,3)
enddo
! Zeroing the total angular momentum of the system
write (iout,*) "vcm right after adjustment:"
write (iout,*) (vcm(j),j=1,3)
endif
- if (.not.rest) then
+
+
+
! call chainbuild
+
+ if ((.not.rest).or.(forceminim)) then
+ if (forceminim) call chainbuild_cart
+ 122 continue
if(iranconf.ne.0 .or.indpdb.gt.0.and..not.unres_pdb .or.preminim) then
if (overlapsc) then
print *, 'Calling OVERLAP_SC'
call minimize(etot,varia,iretcode,nfun)
call var_to_geom(nvar,varia)
endif
+ write(iout,*) "just before minimin"
+ call cartprint
if(me.eq.king.or..not.out1file) &
write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
endif
+ write(iout,*) "just after minimin"
+ call cartprint
if(iranconf.ne.0) then
!c 8/22/17 AL Loop to produce a low-energy random conformation
DO iranmin=1,40
endif
if(me.eq.king.or..not.out1file) &
write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
-
+ write(iout,*) "just after minimin"
+ call cartprint
if (isnan(etot) .or. etot.gt.4.0d6) then
write (iout,*) "Energy too large",etot, &
" trying another random conformation"
stop
#endif
44 continue
+ else if (preminim) then
+ if (start_from_model) then
+ n_model_try=0
+ fail=.true.
+ list_model_try=0
+ do while (fail .and. n_model_try.lt.nmodel_start)
+ write (iout,*) "n_model_try",n_model_try
+ do
+ i_model=iran_num(1,nmodel_start)
+ do k=1,n_model_try
+ if (i_model.eq.list_model_try(k)) exit
+ enddo
+ if (k.gt.n_model_try) exit
+ enddo
+ n_model_try=n_model_try+1
+ list_model_try(n_model_try)=i_model
+ if (me.eq.king .or. .not. out1file) &
+ write (iout,*) 'Trying to start from model ',&
+ pdbfiles_chomo(i_model)(:ilen(pdbfiles_chomo(i_model)))
+ do i=1,2*nres
+ do j=1,3
+ c(j,i)=chomo(j,i,i_model)
+ enddo
+ enddo
+ call int_from_cart(.true.,.false.)
+ call sc_loc_geom(.false.)
+ dc(:,0)=c(:,1)
+ do i=1,nres-1
+ do j=1,3
+ dc(j,i)=c(j,i+1)-c(j,i)
+ dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
+ enddo
+ enddo
+ do i=2,nres-1
+ do j=1,3
+ dc(j,i+nres)=c(j,i+nres)-c(j,i)
+ dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
+ enddo
+ enddo
+ if (me.eq.king.or..not.out1file) then
+ write (iout,*) "Energies before removing overlaps"
+ call etotal(energia(0))
+ call enerprint(energia(0))
+ endif
+! Remove SC overlaps if requested
+ if (overlapsc) then
+ write (iout,*) 'Calling OVERLAP_SC'
+ call overlap_sc(fail)
+ if (fail) then
+ write (iout,*)&
+ "Failed to remove overlap from model",i_model
+ cycle
+ endif
+ endif
+ if (me.eq.king.or..not.out1file) then
+ write (iout,*) "Energies after removing overlaps"
+ call etotal(energia(0))
+ call enerprint(energia(0))
+ endif
+#ifdef SEARCHSC
+! Search for better SC rotamers if requested
+ if (searchsc) then
+ call sc_move(2,nres-1,10,1d10,nft_sc,etot)
+ print *,'SC_move',nft_sc,etot
+ if (me.eq.king.or..not.out1file)&
+ write(iout,*) 'SC_move',nft_sc,etot
+ endif
+ call etotal(energia(0))
+#endif
+ enddo
+ call MPI_Gather(i_model,1,MPI_INTEGER,i_start_models(0),&
+ 1,MPI_INTEGER,king,CG_COMM,IERROR)
+ if (n_model_try.gt.nmodel_start .and.&
+ (me.eq.king .or. out1file)) then
+ write (iout,*)&
+ "All models have irreparable overlaps. Trying randoms starts."
+ iranconf=1
+ i_model=nmodel_start+1
+ goto 122
+ endif
+ else
+! Remove SC overlaps if requested
+ if (overlapsc) then
+ write (iout,*) 'Calling OVERLAP_SC'
+ call overlap_sc(fail)
+ if (fail) then
+ write (iout,*)&
+ "Failed to remove overlap"
+ endif
+ endif
+ if (me.eq.king.or..not.out1file) then
+ write (iout,*) "Energies after removing overlaps"
+ call etotal(energia(0))
+ call enerprint(energia(0))
+ endif
+ endif
+! 8/22/17 AL Minimize initial structure
+ if (dccart) then
+ if (me.eq.king.or..not.out1file) write(iout,*)&
+ 'Minimizing initial PDB structure: Calling MINIM_DC'
+ call minim_dc(etot,iretcode,nfun)
+ else
+ call geom_to_var(nvar,varia)
+ if(me.eq.king.or..not.out1file) write (iout,*)&
+ 'Minimizing initial PDB structure: Calling MINIMIZE.'
+ call minimize(etot,varia,iretcode,nfun)
+ call var_to_geom(nvar,varia)
+#ifdef LBFGS
+ if (me.eq.king.or..not.out1file)&
+ write(iout,*) 'LBFGS return code is ',status,' eval ',nfun
+ if(me.eq.king.or..not.out1file)&
+ write(iout,*) 'LBFGS return code is ',status,' eval ',nfun
+#else
+ if (me.eq.king.or..not.out1file)&
+ write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
+ if(me.eq.king.or..not.out1file)&
+ write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
+#endif
+ endif
+ endif
+ if (nmodel_start.gt.0 .and. me.eq.king) then
+ write (iout,'(a)') "Task Starting model"
+ do i=0,nodes-1
+ if (i_start_models(i).gt.nmodel_start) then
+ write (iout,'(i4,2x,a)') i,"RANDOM STRUCTURE"
+ else
+ write(iout,'(i4,2x,a)')i,pdbfiles_chomo(i_start_models(i)) &
+ (:ilen(pdbfiles_chomo(i_start_models(i))))
+ endif
+ enddo
endif
endif
call chainbuild_cart
call kinetic(EK)
+ write(iout,*) "just after kinetic"
+ call cartprint
if (tbf) then
call verlet_bath
endif
kinetic_T=2.0d0/(dimen3*Rb)*EK
if(me.eq.king.or..not.out1file)then
+ write(iout,*) "just after verlet_bath"
call cartprint
call intout
endif
#endif
potE=potEcomp(0)
call cartgrad
+ write(iout,*) "before lagrangian"
call lagrangian
+ write(iout,*) "before max_accel"
call max_accel
if (amax*d_time .gt. dvmax) then
d_time=d_time*dvmax/amax
call enerprint(potEcomp)
! write(iout,*) (potEcomp(i),i=0,n_ene)
endif
- potE=potEcomp(0)-potEcomp(20)
+ potE=potEcomp(0)-potEcomp(51)
totE=EK+potE
itime=0
itime_mat=itime
write (iout,*) "Initial velocities"
write (iout,"(13x,' backbone ',23x,' side chain')")
do i=0,nres
- write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
+ write (iout,'(i6,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
(d_t(j,i+nres),j=1,3)
enddo
write (iout,*) "Initial accelerations"
! include 'COMMON.NAMES'
! include 'COMMON.TIME1'
real(kind=8) :: xv,sigv,lowb,highb ,Ek1
+#ifdef FIVEDIAG
+ integer ichain,n,innt,inct,ibeg,ierr
+ real(kind=8) ,allocatable, dimension(:):: work
+ integer,allocatable,dimension(:) :: iwork
+! double precision Ghalf(mmaxres2_chain),Geigen(maxres2_chain),&
+! Gvec(maxres2_chain,maxres2_chain)
+! common /przechowalnia/Ghalf,Geigen,Gvec
+!#ifdef DEBUG
+! double precision inertia(maxres2_chain,maxres2_chain)
+!#endif
+#endif
!#define DEBUG
#ifdef FIVEDIAG
- real(kind=8) ,allocatable, dimension(:) :: DDU1,DDU2,DL2,DL1,xsolv,DML,rs
- real(kind=8) :: sumx
+ real(kind=8) ,allocatable, dimension(:) :: xsolv,DML,rs
+ real(kind=8) :: sumx,Ek2,Ek3,aux
#ifdef DEBUG
real(kind=8) ,allocatable, dimension(:) :: rsold
- real (kind=8),allocatable,dimension(:,:) :: matold
+ real (kind=8),allocatable,dimension(:,:) :: matold,inertia
integer :: iti
#endif
#endif
- integer :: i,j,ii,k,ind,mark,imark,mnum
+ integer :: i,j,ii,k,mark,imark,mnum,nres2
+ integer(kind=8) :: ind
! Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
+!#undef DEBUG
! First generate velocities in the eigenspace of the G matrix
! write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
! call flush(iout)
- xv=0.0d0
- ii=0
- do i=1,dimen
- do k=1,3
- ii=ii+1
- sigv=dsqrt((Rb*t_bath)/geigen(i))
- lowb=-5*sigv
- highb=5*sigv
- d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
+#ifdef FIVEDIAG
+ if(.not.allocated(work)) then
+ allocate(work(48*nres))
+ allocate(iwork(6*nres))
+ endif
+ print *,"IN RANDOM VEL"
+ nres2=2*nres
+! print *,size(ghalf)
+#undef DEBUG
#ifdef DEBUG
- write (iout,*) "i",i," ii",ii," geigen",geigen(i),&
- " d_t_work_new",d_t_work_new(ii)
+ write (iout,*) "Random_vel, fivediag"
+ flush(iout)
+ allocate(inertia(2*nres,2*nres))
#endif
- enddo
- enddo
+ d_t=0.0d0
+ Ek2=0.0d0
+ EK=0.0d0
+ Ek3=0.0d0
#ifdef DEBUG
-! diagnostics
- Ek1=0.0d0
- ii=0
- do i=1,dimen
- do k=1,3
- ii=ii+1
- Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
- enddo
- enddo
- write (iout,*) "Ek from eigenvectors",Ek1
- write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
-! end diagnostics
+ write(iout,*), nchain
#endif
-#ifdef FIVEDIAG
-! Transform velocities to UNRES coordinate space
- allocate (DL1(2*nres))
- allocate (DDU1(2*nres))
- allocate (DL2(2*nres))
- allocate (DDU2(2*nres))
- allocate (xsolv(2*nres))
- allocate (DML(2*nres))
- allocate (rs(2*nres))
+ do ichain=1,nchain
+ ind=0
+ if(.not.allocated(ghalf)) print *,"COCO"
+ if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2))
+ ghalf=0.0d0
+ n=dimen_chain(ichain)
+ innt=iposd_chain(ichain)
+ inct=innt+n-1
#ifdef DEBUG
- allocate (rsold(2*nres))
- allocate (matold(2*nres,2*nres))
- matold=0
- matold(1,1)=DMorig(1)
- matold(1,2)=DU1orig(1)
- matold(1,3)=DU2orig(1)
- write (*,*) DMorig(1),DU1orig(1),DU2orig(1)
-
- do i=2,dimen-1
- do j=1,dimen
- if (i.eq.j) then
- matold(i,j)=DMorig(i)
- matold(i,j-1)=DU1orig(i-1)
- if (j.ge.3) then
- matold(i,j-2)=DU2orig(i-2)
- endif
-
- endif
-
- enddo
- do j=1,dimen-1
- if (i.eq.j) then
- matold(i,j+1)=DU1orig(i)
-
- end if
- enddo
- do j=1,dimen-2
- if(i.eq.j) then
- matold(i,j+2)=DU2orig(i)
- endif
- enddo
- enddo
- matold(dimen,dimen)=DMorig(dimen)
- matold(dimen,dimen-1)=DU1orig(dimen-1)
- matold(dimen,dimen-2)=DU2orig(dimen-2)
- write(iout,*) "old gmatrix"
- call matout(dimen,dimen,2*nres,2*nres,matold)
+ write (iout,*) "Chain",ichain," n",n," start",innt
+ do i=innt,inct
+ if (i.lt.inct-1) then
+ write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i),DU1orig(i),&
+ DU2orig(i)
+ else if (i.eq.inct-1) then
+ write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i),DU1orig(i)
+ else
+ write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i)
+ endif
+ enddo
#endif
- d_t_work=0.0d0
- do i=1,dimen
-! Find the ith eigenvector of the pentadiagonal inertiq matrix
- do imark=dimen,1,-1
- do j=1,imark-1
- DML(j)=DMorig(j)-geigen(i)
- enddo
- do j=imark+1,dimen
- DML(j-1)=DMorig(j)-geigen(i)
- enddo
- do j=1,imark-2
- DDU1(j)=DU1orig(j)
- enddo
- DDU1(imark-1)=DU2orig(imark-1)
- do j=imark+1,dimen-1
- DDU1(j-1)=DU1orig(j)
- enddo
- do j=1,imark-3
- DDU2(j)=DU2orig(j)
- enddo
- DDU2(imark-2)=0.0d0
- DDU2(imark-1)=0.0d0
- do j=imark,dimen-3
- DDU2(j)=DU2orig(j+1)
- enddo
- do j=1,dimen-3
- DL2(j+2)=DDU2(j)
- enddo
- do j=1,dimen-2
- DL1(j+1)=DDU1(j)
- enddo
-#ifdef DEBUG
- write (iout,*) "DL2,DL1,DML,DDU1,DDU2"
- write(iout,'(10f10.5)') (DL2(k),k=3,dimen-1)
- write(iout,'(10f10.5)') (DL1(k),k=2,dimen-1)
- write(iout,'(10f10.5)') (DML(k),k=1,dimen-1)
- write(iout,'(10f10.5)') (DDU1(k),k=1,dimen-2)
- write(iout,'(10f10.5)') (DDU2(k),k=1,dimen-3)
-#endif
- rs=0.0d0
- if (imark.gt.2) rs(imark-2)=-DU2orig(imark-2)
- if (imark.gt.1) rs(imark-1)=-DU1orig(imark-1)
- if (imark.lt.dimen) rs(imark)=-DU1orig(imark)
- if (imark.lt.dimen-1) rs(imark+1)=-DU2orig(imark)
-#ifdef DEBUG
- rsold=rs
-#endif
-! write (iout,*) "Vector rs"
-! do j=1,dimen-1
-! write (iout,*) j,rs(j)
-! enddo
- xsolv=0.0d0
- call FDIAG(dimen-1,DL2,DL1,DML,DDU1,DDU2,rs,xsolv,mark)
-
- if (mark.eq.1) then
-
-#ifdef DEBUG
-! Check solution
- do j=1,imark-1
- sumx=-geigen(i)*xsolv(j)
- do k=1,imark-1
- sumx=sumx+matold(j,k)*xsolv(k)
- enddo
- do k=imark+1,dimen
- sumx=sumx+matold(j,k)*xsolv(k-1)
- enddo
- write(iout,'(i5,3f10.5)') j,sumx,rsold(j),sumx-rsold(j)
- enddo
- do j=imark+1,dimen
- sumx=-geigen(i)*xsolv(j-1)
- do k=1,imark-1
- sumx=sumx+matold(j,k)*xsolv(k)
- enddo
- do k=imark+1,dimen
- sumx=sumx+matold(j,k)*xsolv(k-1)
- enddo
- write(iout,'(i5,3f10.5)') j-1,sumx,rsold(j-1),sumx-rsold(j-1)
- enddo
-! end check
- write (iout,*)&
- "Solution of equations system",i," for eigenvalue",geigen(i)
- do j=1,dimen-1
- write(iout,'(i5,f10.5)') j,xsolv(j)
- enddo
-#endif
- do j=dimen-1,imark,-1
- xsolv(j+1)=xsolv(j)
- enddo
- xsolv(imark)=1.0d0
+ ghalf(ind+1)=dmorig(innt)
+ ghalf(ind+2)=du1orig(innt)
+ ghalf(ind+3)=dmorig(innt+1)
+ ind=ind+3
+ do i=3,n
+ ind=ind+i-3
+ write (iout,*) "i",i," ind",ind," indu2",innt+i-2,&
+ " indu1",innt+i-1," indm",innt+i
+ ghalf(ind+1)=du2orig(innt-1+i-2)
+ ghalf(ind+2)=du1orig(innt-1+i-1)
+ ghalf(ind+3)=dmorig(innt-1+i)
+!c write (iout,'(3(a,i2,1x))') "DU2",innt-1+i-2,
+!c & "DU1",innt-1+i-1,"DM ",innt-1+i
+ ind=ind+3
+ enddo
#ifdef DEBUG
- write (iout,*) "Un-normalized eigenvector",i," for eigenvalue",geigen(i)
- do j=1,dimen
- write(iout,'(i5,f10.5)') j,xsolv(j)
- enddo
+ ind=0
+ do i=1,n
+ do j=1,i
+ ind=ind+1
+ inertia(i,j)=ghalf(ind)
+ inertia(j,i)=ghalf(ind)
+ enddo
+ enddo
#endif
-! Normalize ith eigenvector
- sumx=0.0d0
- do j=1,dimen
- sumx=sumx+xsolv(j)**2
- enddo
- sumx=dsqrt(sumx)
- do j=1,dimen
- xsolv(j)=xsolv(j)/sumx
- enddo
#ifdef DEBUG
- write (iout,*) "Eigenvector",i," for eigenvalue",geigen(i)
- do j=1,dimen
- write(iout,'(i5,3f10.5)') j,xsolv(j)
- enddo
+ write (iout,*) "Chain ",ichain," ind",ind," dim",n*(n+1)/2
+ write (iout,*) "Five-diagonal inertia matrix, lower triangle"
+! call matoutr(n,ghalf)
#endif
-! All done at this point for eigenvector i, exit loop
- exit
-
- endif ! mark.eq.1
-
- enddo ! imark
-
- if (mark.ne.1) then
- write (iout,*) "Unable to find eigenvector",i
- endif
-
-! write (iout,*) "i=",i
- do k=1,3
-! write (iout,*) "k=",k
- do j=1,dimen
- ind=(j-1)*3+k
-! write(iout,*) "ind",ind," ind1",3*(i-1)+k
- d_t_work(ind)=d_t_work(ind) &
- +xsolv(j)*d_t_work_new(3*(i-1)+k)
- enddo
- enddo
- enddo ! i (loop over the eigenvectors)
-
-#ifdef DEBUG
- write (iout,*) "d_t_work"
- do i=1,3*dimen
- write (iout,"(i5,f10.5)") i,d_t_work(i)
- enddo
- Ek1=0.0d0
- ii=0
- do i=nnt,nct
-! if (itype(i,1).eq.10) then
- mnum=molnum(i)
- if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
- iti=iabs(itype(i,mnum))
-! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
- if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
- .or.(mnum.eq.5)) then
- j=ii+3
- else
- j=ii+6
+ call gldiag(nres*2,n,n,Ghalf,work,Geigen,Gvec,ierr,iwork)
+ if (large) then
+ write (iout,'(//a,i3)')&
+ "Eigenvectors and eigenvalues of the G matrix chain",ichain
+ call eigout(n,n,nres*2,nres*2,Gvec,Geigen)
endif
- if (i.lt.nct) then
+#ifdef DIAGCHECK
+!c check diagonalization
+ do i=1,n
+ do j=1,n
+ aux=0.0d0
+ do k=1,n
+ do l=1,n
+ aux=aux+gvec(k,i)*gvec(l,j)*inertia(k,l)
+ enddo
+ enddo
+ if (i.eq.j) then
+ write (iout,*) i,j,aux,geigen(i)
+ else
+ write (iout,*) i,j,aux
+ endif
+ enddo
+ enddo
+#endif
+ xv=0.0d0
+ ii=0
+ do i=1,n
do k=1,3
-! write (iout,*) "k",k," ii+k",ii+k," ii+j+k",ii+j+k,"EK1",Ek1
- Ek1=Ek1+0.5d0*mp(mnum)*((d_t_work(ii+j+k)+d_t_work_new(ii+k))/2)**2+&
- 0.5d0*0.25d0*IP(mnum)*(d_t_work(ii+j+k)-d_t_work_new(ii+k))**2
+ ii=ii+1
+ sigv=dsqrt((Rb*t_bath)/geigen(i))
+ lowb=-5*sigv
+ highb=5*sigv
+ d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
+ EK=EK+0.5d0*geigen(i)*d_t_work_new(ii)**2
+!c write (iout,*) "i",i," ii",ii," geigen",geigen(i),
+!c & " d_t_work_new",d_t_work_new(ii)
enddo
- endif
- if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
- .and.(mnum.ne.5)) ii=ii+3
- write (iout,*) "i",i," itype",itype(i,mnum)," mass",msc(itype(i,mnum),mnum)
- write (iout,*) "ii",ii
+ enddo
do k=1,3
- ii=ii+1
- write (iout,*) "k",k," ii",ii,"EK1",EK1
- if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
- .and.(mnum.ne.5))&
- Ek1=Ek1+0.5d0*Isc(iabs(itype(i,mnum)),mnum)*(d_t_work(ii)-d_t_work(ii-3))**2
- Ek1=Ek1+0.5d0*msc(iabs(itype(i,mnum)),mnum)*d_t_work(ii)**2
+ do i=1,n
+ ind=(i-1)*3+k
+ d_t_work(ind)=0.0d0
+ do j=1,n
+ d_t_work(ind)=d_t_work(ind)&
+ +Gvec(i,j)*d_t_work_new((j-1)*3+k)
+ enddo
+!c write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
+!c call flush(iout)
+ enddo
enddo
- write (iout,*) "i",i," ii",ii
- enddo
- write (iout,*) "Ek from d_t_work",Ek1
- write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
-#endif
- if(allocated(DDU1)) deallocate(DDU1)
- if(allocated(DDU2)) deallocate(DDU2)
- if(allocated(DL2)) deallocate(DL2)
- if(allocated(DL1)) deallocate(DL1)
- if(allocated(xsolv)) deallocate(xsolv)
- if(allocated(DML)) deallocate(DML)
- if(allocated(rs)) deallocate(rs)
#ifdef DEBUG
- if(allocated(matold)) deallocate(matold)
- if(allocated(rsold)) deallocate(rsold)
-#endif
- ind=1
- do j=nnt,nct
+ aux=0.0d0
do k=1,3
- d_t(k,j)=d_t_work(ind)
- ind=ind+1
+ do i=1,n
+ do j=1,n
+ aux=aux+inertia(i,j)*d_t_work(3*(i-1)+k)*d_t_work(3*(j-1)+k)
+ enddo
+ enddo
enddo
- mnum=molnum(j)
- if (itype(j,1).ne.10 .and. itype(j,mnum).ne.ntyp1_molec(mnum)&
- .and.(mnum.ne.5)) then
- do k=1,3
- d_t(k,j+nres)=d_t_work(ind)
+ Ek3=Ek3+aux/2
+#endif
+!c Transfer to the d_t vector
+ innt=chain_border(1,ichain)
+ inct=chain_border(2,ichain)
+ ind=0
+!c write (iout,*) "ichain",ichain," innt",innt," inct",inct
+ do i=innt,inct
+ do j=1,3
ind=ind+1
+ d_t(j,i)=d_t_work(ind)
enddo
- end if
+ mnum=molnum(i)
+ if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum).and.mnum.le.2) then
+ do j=1,3
+ ind=ind+1
+ d_t(j,i+nres)=d_t_work(ind)
+ enddo
+ endif
+ enddo
enddo
-#ifdef DEBUG
- write (iout,*) "Random velocities in the Calpha,SC space"
- do i=nnt,nct
- write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3)
+ if (large) then
+ write (iout,*)
+ write (iout,*) "Random velocities in the Calpha,SC space"
+ do i=1,nres
+ mnum=molnum(i)
+ write (iout,'(a3,1h(,i5,1h),3f10.5,3x,3f10.5)')&
+ restyp(itype(i,mnum),mnum),i,(d_t(j,i),j=1,3),(d_t(j,i+nres),j=1,3)
+ enddo
+ endif
+ call kinetic_CASC(Ek1)
+!
+! Transform the velocities to virtual-bond space
+!
+#define WLOS
+#ifdef WLOS
+ if (nnt.eq.1) then
+ d_t(:,0)=d_t(:,1)
+ endif
+ do i=1,nres
+ mnum=molnum(i)
+ if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum).or.mnum.ge.3) then
+ do j=1,3
+ d_t(j,i)=d_t(j,i+1)-d_t(j,i)
+ enddo
+ else
+ do j=1,3
+ d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
+ d_t(j,i)=d_t(j,i+1)-d_t(j,i)
+ enddo
+ end if
enddo
- do i=nnt,nct
- write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3)
+ d_t(:,nres)=0.0d0
+ d_t(:,nct)=0.0d0
+ d_t(:,2*nres)=0.0d0
+ if (nnt.gt.1) then
+ d_t(:,0)=d_t(:,1)
+ d_t(:,1)=0.0d0
+ endif
+!c d_a(:,0)=d_a(:,1)
+!c d_a(:,1)=0.0d0
+!c write (iout,*) "Shifting accelerations"
+ do ichain=2,nchain
+!c write (iout,*) "ichain",chain_border1(1,ichain)-1,
+!c & chain_border1(1,ichain)
+ d_t(:,chain_border1(1,ichain)-1)=d_t(:,chain_border1(1,ichain))
+ d_t(:,chain_border1(1,ichain))=0.0d0
+ enddo
+!c write (iout,*) "Adding accelerations"
+ do ichain=2,nchain
+!c write (iout,*) "chain",ichain,chain_border1(1,ichain)-1,
+!c & chain_border(2,ichain-1)
+ d_t(:,chain_border1(1,ichain)-1)=&
+ d_t(:,chain_border1(1,ichain)-1)+d_t(:,chain_border(2,ichain-1))
+ d_t(:,chain_border(2,ichain-1))=0.0d0
+ enddo
+ do ichain=2,nchain
+ write (iout,*) "chain",ichain,chain_border1(1,ichain)-1,&
+ chain_border(2,ichain-1)
+ d_t(:,chain_border1(1,ichain)-1)=&
+ d_t(:,chain_border1(1,ichain)-1)+d_t(:,chain_border(2,ichain-1))
+ d_t(:,chain_border(2,ichain-1))=0.0d0
enddo
-#endif
+#else
+ ibeg=0
+!c do j=1,3
+!c d_t(j,0)=d_t(j,nnt)
+!c enddo
+ do ichain=1,nchain
+ innt=chain_border(1,ichain)
+ inct=chain_border(2,ichain)
+!c write (iout,*) "ichain",ichain," innt",innt," inct",inct
+!c write (iout,*) "ibeg",ibeg
do j=1,3
- d_t(j,0)=d_t(j,nnt)
+ d_t(j,ibeg)=d_t(j,innt)
enddo
- do i=nnt,nct
-! if (itype(i,1).eq.10) then
- mnum=molnum(i)
- if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
- .or.(mnum.eq.5)) then
+ ibeg=inct+1
+ do i=innt,inct
+ mnum=molnum(i)
+ if (iabs(itype(i,1).eq.10).or.mnum.ge.3) then
+!c write (iout,*) "i",i,(d_t(j,i),j=1,3),(d_t(j,i+1),j=1,3)
do j=1,3
d_t(j,i)=d_t(j,i+1)-d_t(j,i)
enddo
enddo
end if
enddo
+ enddo
+#endif
+ if (large) then
+ write (iout,*)
+ write (iout,*)&
+ "Random velocities in the virtual-bond-vector space"
+ write (iout,'(3hORG,1h(,i5,1h),3f10.5)') 0,(d_t(j,0),j=1,3)
+ do i=1,nres
+ write (iout,'(a3,1h(,i5,1h),3f10.5,3x,3f10.5)')&
+ restyp(itype(i,mnum),mnum),i,(d_t(j,i),j=1,3),(d_t(j,i+nres),j=1,3)
+ enddo
+ write (iout,*)
+ write (iout,*) "Kinetic energy from inertia matrix eigenvalues",&
+ Ek
+ write (iout,*)&
+ "Kinetic temperatures from inertia matrix eigenvalues",&
+ 2*Ek/(3*dimen*Rb)
#ifdef DEBUG
- write (iout,*)"Random velocities in the virtual-bond-vector space"
- do i=nnt,nct-1
- write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3)
+ write (iout,*) "Kinetic energy from inertia matrix",Ek3
+ write (iout,*) "Kinetic temperatures from inertia",&
+ 2*Ek3/(3*dimen*Rb)
+#endif
+ write (iout,*) "Kinetic energy from velocities in CA-SC space",&
+ Ek1
+ write (iout,*)&
+ "Kinetic temperatures from velovities in CA-SC space",&
+ 2*Ek1/(3*dimen*Rb)
+ call kinetic(Ek1)
+ write (iout,*)&
+ "Kinetic energy from virtual-bond-vector velocities",Ek1
+ write (iout,*)&
+ "Kinetic temperature from virtual-bond-vector velocities ",&
+ 2*Ek1/(dimen3*Rb)
+ endif
+#else
+ xv=0.0d0
+ ii=0
+ do i=1,dimen
+ do k=1,3
+ ii=ii+1
+ sigv=dsqrt((Rb*t_bath)/geigen(i))
+ lowb=-5*sigv
+ highb=5*sigv
+ d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
+#ifdef DEBUG
+ write (iout,*) "i",i," ii",ii," geigen",geigen(i),&
+ " d_t_work_new",d_t_work_new(ii)
+#endif
+ enddo
enddo
- do i=nnt,nct
- write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3)
+#ifdef DEBUG
+! diagnostics
+ Ek1=0.0d0
+ ii=0
+ do i=1,dimen
+ do k=1,3
+ ii=ii+1
+ Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
+ enddo
enddo
- call kinetic(Ek1)
- write (iout,*) "Ek from d_t_work",Ek1
+ write (iout,*) "Ek from eigenvectors",Ek1
write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
+! end diagnostics
#endif
-#else
+
do k=0,2
do i=1,dimen
ind=(i-1)*3+k+1
mnum=molnum(i)
! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
- .and.(mnum.ne.5)) then
+ .and.(mnum.lt.4)) then
do j=1,3
ind=ind+1
d_t(j,i+nres)=d_t_work(ind)
! 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
! call flush(iout)
! write(iout,*) "end init MD"
+#undef DEBUG
return
end subroutine random_vel
!-----------------------------------------------------------------------------
mnum=molnum(i)
! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
- .and.(mnum.ne.5)) then
+ .and.(mnum.lt.4)) then
do j=1,3
dc_work(ind+j)=dc_old(j,i+nres)
d_t_work(ind+j)=d_t_old(j,i+nres)
mnum=molnum(i)
! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
- .and.(mnum.ne.5)) then
+ .and.(mnum.lt.4)) then
inres=i+nres
do j=1,3
dc(j,inres)=dc_work(ind+j)
mnum=molnum(i)
! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
- .and.(mnum.ne.5)) then
+ .and.(mnum.lt.4)) then
inres=i+nres
do j=1,3
d_t(j,inres)=d_t_work(ind+j)
" vrand_mat2"
do i=1,dimen
do j=1,dimen
- write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
- vfric_mat(i,j),afric_mat(i,j),&
+ write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),&
+ vfric_mat(i,j),afric_mat(i,j),&
prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
enddo
enddo
mnum=molnum(i)
! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
- .and.(mnum.ne.5)) then
+ .and.(mnum.lt.4)) then
inres=i+nres
do j=1,3
dc(j,inres)=dc_work(ind+j)
enddo
d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
enddo
-#endif
+#endif
+ do j=1,3
+ d_t(j,0)=d_t_work(j)
+ enddo
+ ind=3
+ do i=nnt,nct-1
+ do j=1,3
+ d_t(j,i)=d_t_work(ind+j)
+ enddo
+ ind=ind+3
+ enddo
+ do i=nnt,nct
+ mnum=molnum(i)
+ if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
+ .and.(mnum.lt.4))
+! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
+ inres=i+nres
+ do j=1,3
+ d_t(j,inres)=d_t_work(ind+j)
+ enddo
+ ind=ind+3
+ endif
+ enddo
+ return
+ end subroutine sd_verlet2_ciccotti
+#endif
+!-----------------------------------------------------------------------------
+! moments.f
+!-----------------------------------------------------------------------------
+#ifdef FIVEDIAG
+ subroutine inertia_tensor
+ use comm_gucio
+ use energy_data
+ real(kind=8) Im(3,3),Imcp(3,3),pr(3),M_SC,&
+ eigvec(3,3),Id(3,3),eigval(3),L(3),vp(3),vrot(3),&
+ vpp(3,0:MAXRES),vs_p(3),pr1(3,3),&
+ pr2(3,3),pp(3),incr(3),v(3),mag,mag2,M_PEP
+ integer iti,inres,i,j,k,mnum,mnum1
+ do i=1,3
+ do j=1,3
+ Im(i,j)=0.0d0
+ pr1(i,j)=0.0d0
+ pr2(i,j)=0.0d0
+ enddo
+ L(i)=0.0d0
+ cm(i)=0.0d0
+ vrot(i)=0.0d0
+ enddo
+ M_PEP=0.0d0
+
+!c caulating the center of the mass of the protein
+ do i=nnt,nct-1
+ mnum=molnum(i)
+ mnum1=molnum(i+1)
+ if (itype(i,mnum).eq.ntyp1_molec(mnum)&
+ .or. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) cycle
+! if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
+ if (mnum.ge.5) mp(mnum)=0.0d0
+ M_PEP=M_PEP+mp(mnum)
+
+ do j=1,3
+ cm(j)=cm(j)+(c(j,i)+0.5d0*dc(j,i))*mp(mnum)
+ enddo
+ enddo
+! do j=1,3
+! cm(j)=mp*cm(j)
+! enddo
+ M_SC=0.0d0
+ do i=nnt,nct
+ mnum=molnum(i)
+ mnum1=molnum(i+1)
+ iti=iabs(itype(i,mnum))
+ if (iti.eq.ntyp1_molec(mnum)) cycle
+ M_SC=M_SC+msc(iabs(iti),mnum)
+ inres=i+nres
+ do j=1,3
+ cm(j)=cm(j)+msc(iabs(iti),mnum)*c(j,inres)
+ enddo
+ enddo
+ do j=1,3
+ cm(j)=cm(j)/(M_SC+M_PEP)
+ enddo
+ do i=nnt,nct-1
+ mnum=molnum(i)
+ mnum1=molnum(i+1)
+! if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
+ if (mnum.ge.5) mp(mnum)=0.0d0
+ if (itype(i,mnum).eq.ntyp1_molec(mnum)&
+ .or. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) cycle
+ do j=1,3
+ pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
+ enddo
+ Im(1,1)=Im(1,1)+mp(mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
+ Im(1,2)=Im(1,2)-mp(mnum)*pr(1)*pr(2)
+ Im(1,3)=Im(1,3)-mp(mnum)*pr(1)*pr(3)
+ Im(2,3)=Im(2,3)-mp(mnum)*pr(2)*pr(3)
+ Im(2,2)=Im(2,2)+mp(mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
+ Im(3,3)=Im(3,3)+mp(mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
+ enddo
+
+ do i=nnt,nct
+ mnum=molnum(i)
+ iti=iabs(itype(i,mnum))
+ if (iti.eq.ntyp1_molec(mnum)) cycle
+ inres=i+nres
+ do j=1,3
+ pr(j)=c(j,inres)-cm(j)
+ enddo
+ Im(1,1)=Im(1,1)+msc(iabs(iti),mnum)*(pr(2)*pr(2)+pr(3)*pr(3))
+ Im(1,2)=Im(1,2)-msc(iabs(iti),mnum)*pr(1)*pr(2)
+ Im(1,3)=Im(1,3)-msc(iabs(iti),mnum)*pr(1)*pr(3)
+ Im(2,3)=Im(2,3)-msc(iabs(iti),mnum)*pr(2)*pr(3)
+ Im(2,2)=Im(2,2)+msc(iabs(iti),mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
+ Im(3,3)=Im(3,3)+msc(iabs(iti),mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
+ enddo
+ do i=nnt,nct-1
+ mnum=molnum(i)
+ mnum1=molnum(i+1)
+ if (itype(i,mnum).eq.ntyp1_molec(mnum)&
+ .or. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) cycle
+ Im(1,1)=Im(1,1)+Ip(mnum)*(1-dc_norm(1,i)*dc_norm(1,i))*&
+ vbld(i+1)*vbld(i+1)*0.25d0
+ Im(1,2)=Im(1,2)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(2,i))*&
+ vbld(i+1)*vbld(i+1)*0.25d0
+ Im(1,3)=Im(1,3)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(3,i))*&
+ vbld(i+1)*vbld(i+1)*0.25d0
+ Im(2,3)=Im(2,3)+Ip(mnum)*(-dc_norm(2,i)*dc_norm(3,i))*&
+ vbld(i+1)*vbld(i+1)*0.25d0
+ Im(2,2)=Im(2,2)+Ip(mnum)*(1-dc_norm(2,i)*dc_norm(2,i))*&
+ vbld(i+1)*vbld(i+1)*0.25d0
+ Im(3,3)=Im(3,3)+Ip(mnum)*(1-dc_norm(3,i)*dc_norm(3,i))*&
+ vbld(i+1)*vbld(i+1)*0.25d0
+ enddo
+ do i=nnt,nct
+ mnum=molnum(i)
+ mnum1=molnum(i+1)
+ iti=iabs(itype(i,mnum))
+ if (iti.ne.10 .and. iti.ne.ntyp1_molec(mnum).and.mnum.le.2) then
+ inres=i+nres
+ Im(1,1)=Im(1,1)+Isc(iti,mnum)*(1-dc_norm(1,inres)*&
+ dc_norm(1,inres))*vbld(inres)*vbld(inres)
+ Im(1,2)=Im(1,2)-Isc(iti,mnum)*(dc_norm(1,inres)*&
+ dc_norm(2,inres))*vbld(inres)*vbld(inres)
+ Im(1,3)=Im(1,3)-Isc(iti,mnum)*(dc_norm(1,inres)*&
+ dc_norm(3,inres))*vbld(inres)*vbld(inres)
+ Im(2,3)=Im(2,3)-Isc(iti,mnum)*(dc_norm(2,inres)*&
+ dc_norm(3,inres))*vbld(inres)*vbld(inres)
+ Im(2,2)=Im(2,2)+Isc(iti,mnum)*(1-dc_norm(2,inres)*&
+ dc_norm(2,inres))*vbld(inres)*vbld(inres)
+ Im(3,3)=Im(3,3)+Isc(iti,mnum)*(1-dc_norm(3,inres)*&
+ dc_norm(3,inres))*vbld(inres)*vbld(inres)
+ endif
+ enddo
+
+ call angmom(cm,L)
+ Im(2,1)=Im(1,2)
+ Im(3,1)=Im(1,3)
+ Im(3,2)=Im(2,3)
+
+!c Copng the Im matrix for the djacob subroutine
+ do i=1,3
+ do j=1,3
+ Imcp(i,j)=Im(i,j)
+ Id(i,j)=0.0d0
+ enddo
+ enddo
+!c Finding the eigenvectors and eignvalues of the inertia tensor
+ call djacob(3,3,10000,1.0d-10,Imcp,eigvec,eigval)
+ do i=1,3
+ if (dabs(eigval(i)).gt.1.0d-15) then
+ Id(i,i)=1.0d0/eigval(i)
+ else
+ Id(i,i)=0.0d0
+ endif
+ enddo
+ do i=1,3
+ do j=1,3
+ Imcp(i,j)=eigvec(j,i)
+ enddo
+ enddo
+ do i=1,3
+ do j=1,3
+ do k=1,3
+ pr1(i,j)=pr1(i,j)+Id(i,k)*Imcp(k,j)
+ enddo
+ enddo
+ enddo
+ do i=1,3
+ do j=1,3
+ do k=1,3
+ pr2(i,j)=pr2(i,j)+eigvec(i,k)*pr1(k,j)
+ enddo
+ enddo
+ enddo
+!c Calculating the total rotational velocity of the molecule
+ do i=1,3
+ do j=1,3
+ vrot(i)=vrot(i)+pr2(i,j)*L(j)
+ enddo
+ enddo
+ do i=nnt,nct-1
+ mnum=molnum(i)
+ mnum1=molnum(i+1)
+! write (iout,*) itype(i+1,mnum1),itype(i,mnum)
+ if (itype(i+1,mnum1).ne.ntyp1_molec(mnum1) &
+ .and. itype(i,mnum).eq.ntyp1_molec(mnum) .or.&
+ itype(i,mnum).ne.ntyp1_molec(mnum) &
+ .and. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) cycle
+ call vecpr(vrot(1),dc(1,i),vp)
+ do j=1,3
+ d_t(j,i)=d_t(j,i)-vp(j)
+ enddo
+ enddo
+ do i=nnt,nct
+ mnum=molnum(i)
+ if(itype(i,1).ne.10 .and.itype(i,mnum).ne.ntyp1_molec(mnum)&
+ .and.mnum.le.2) then
+ inres=i+nres
+ call vecpr(vrot(1),dc(1,inres),vp)
+ do j=1,3
+ d_t(j,inres)=d_t(j,inres)-vp(j)
+ enddo
+ endif
+ enddo
+ call angmom(cm,L)
+ return
+ end subroutine inertia_tensor
+!------------------------------------------------------------
+ subroutine angmom(cm,L)
+ implicit none
+ double precision L(3),cm(3),pr(3),vp(3),vrot(3),incr(3),v(3),&
+ pp(3),mscab
+ integer iti,inres,i,j,mnum,mnum1
+!c Calculate the angular momentum
+ do j=1,3
+ L(j)=0.0d0
+ enddo
+ do j=1,3
+ incr(j)=d_t(j,0)
+ enddo
+ do i=nnt,nct-1
+ mnum=molnum(i)
+ mnum1=molnum(i+1)
+! if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
+ if (mnum.ge.5) mp(mnum)=0.0d0
+ if (itype(i,mnum).eq.ntyp1_molec(mnum)&
+ .or. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) cycle
+ do j=1,3
+ pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
+ enddo
+ do j=1,3
+ v(j)=incr(j)+0.5d0*d_t(j,i)
+ enddo
+ do j=1,3
+ incr(j)=incr(j)+d_t(j,i)
+ enddo
+ call vecpr(pr(1),v(1),vp)
+ do j=1,3
+ L(j)=L(j)+mp(mnum)*vp(j)
+ enddo
+ do j=1,3
+ pr(j)=0.5d0*dc(j,i)
+ pp(j)=0.5d0*d_t(j,i)
+ enddo
+ call vecpr(pr(1),pp(1),vp)
+ do j=1,3
+ L(j)=L(j)+Ip(mnum)*vp(j)
+ enddo
+ enddo
do j=1,3
- d_t(j,0)=d_t_work(j)
+ incr(j)=d_t(j,0)
enddo
- ind=3
- do i=nnt,nct-1
+ do i=nnt,nct
+ mnum=molnum(i)
+ iti=iabs(itype(i,mnum))
+ inres=i+nres
do j=1,3
- d_t(j,i)=d_t_work(ind+j)
+ pr(j)=c(j,inres)-cm(j)
enddo
- ind=ind+3
- enddo
- do i=nnt,nct
- mnum=molnum(i)
- if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
- .and.(mnum.ne.5))
-! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
- inres=i+nres
+ if (itype(i,1).ne.10 .and.itype(i,mnum).ne.ntyp1_molec(mnum)&
+ .and.mnum.le.2) then
do j=1,3
- d_t(j,inres)=d_t_work(ind+j)
+ v(j)=incr(j)+d_t(j,inres)
+ enddo
+ else
+ do j=1,3
+ v(j)=incr(j)
enddo
- ind=ind+3
endif
- enddo
+ call vecpr(pr(1),v(1),vp)
+!c write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),
+!c & " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
+! if (mnum.gt.4) then
+! mscab=0.0d0
+! else
+ mscab=msc(iti,mnum)
+! endif
+ do j=1,3
+ L(j)=L(j)+mscab*vp(j)
+ enddo
+!c write (iout,*) "L",(l(j),j=1,3)
+ if (itype(i,1).ne.10 .and.itype(i,mnum).ne.ntyp1_molec(mnum)&
+ .and.mnum.le.2) then
+ do j=1,3
+ v(j)=incr(j)+d_t(j,inres)
+ enddo
+ call vecpr(dc(1,inres),d_t(1,inres),vp)
+ do j=1,3
+ L(j)=L(j)+Isc(iti,mnum)*vp(j)
+ enddo
+ endif
+ do j=1,3
+ incr(j)=incr(j)+d_t(j,i)
+ enddo
+ enddo
return
- end subroutine sd_verlet2_ciccotti
-#endif
-!-----------------------------------------------------------------------------
-! moments.f
-!-----------------------------------------------------------------------------
+ end subroutine angmom
+!---------------------------------------------------------------
+ subroutine vcm_vel(vcm)
+ double precision vcm(3),vv(3),summas,amas
+ integer i,j,mnum,mnum1
+ do j=1,3
+ vcm(j)=0.0d0
+ vv(j)=d_t(j,0)
+ enddo
+ summas=0.0d0
+ do i=nnt,nct
+ mnum=molnum(i)
+ if ((mnum.ge.5).or.(mnum.eq.3))&
+ mp(mnum)=msc(itype(i,mnum),mnum)
+ if (i.lt.nct) then
+ summas=summas+mp(mnum)
+ do j=1,3
+ vcm(j)=vcm(j)+mp(mnum)*(vv(j)+0.5d0*d_t(j,i))
+ enddo
+ endif
+ if (mnum.ne.4) then
+ amas=msc(iabs(itype(i,mnum)),mnum)
+ else
+ amas=0.0d0
+ endif
+! amas=msc(iabs(itype(i)))
+ summas=summas+amas
+! if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
+ .and.(mnum.lt.4)) then
+ do j=1,3
+ vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres))
+ enddo
+ else
+ do j=1,3
+ vcm(j)=vcm(j)+amas*vv(j)
+ enddo
+ endif
+ do j=1,3
+ vv(j)=vv(j)+d_t(j,i)
+ enddo
+ enddo
+!c write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas
+ do j=1,3
+ vcm(j)=vcm(j)/summas
+ enddo
+ return
+ end subroutine vcm_vel
+#else
subroutine inertia_tensor
! Calculating the intertia tensor for the entire protein in order to
M_PEP=0.0d0
do i=nnt,nct-1
mnum=molnum(i)
- if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
- if (itype(i,mnum).eq.ntyp1_molec(mnum)) cycle
+ if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
+ write(iout,*) "WTF",itype(i,mnum),i,mnum,mp(mnum)
+! if (itype(i,mnum).eq.ntyp1_molec(mnum)) cycle
M_PEP=M_PEP+mp(mnum)
+
do j=1,3
cm(j)=cm(j)+(c(j,i)+0.5d0*dc(j,i))*mp(mnum)
enddo
M_SC=0.0d0
do i=nnt,nct
mnum=molnum(i)
- if (mnum.eq.5) cycle
+ if (mnum.ge.5) cycle
iti=iabs(itype(i,mnum))
M_SC=M_SC+msc(iabs(iti),mnum)
inres=i+nres
do j=1,3
cm(j)=cm(j)/(M_SC+M_PEP)
enddo
-
+! write(iout,*) "Center of mass:",cm
do i=nnt,nct-1
mnum=molnum(i)
- if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
+ if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
do j=1,3
pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
enddo
Im(2,2)=Im(2,2)+mp(mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
Im(3,3)=Im(3,3)+mp(mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
enddo
+
+! write(iout,*) "The angular momentum before msc add"
+! do i=1,3
+! write (iout,*) (Im(i,j),j=1,3)
+! enddo
do i=nnt,nct
mnum=molnum(i)
iti=iabs(itype(i,mnum))
- if (mnum.eq.5) cycle
+ if (mnum.ge.5) cycle
inres=i+nres
do j=1,3
pr(j)=c(j,inres)-cm(j)
Im(2,2)=Im(2,2)+msc(iabs(iti),mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
Im(3,3)=Im(3,3)+msc(iabs(iti),mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
enddo
+! write(iout,*) "The angular momentum before Ip add"
+! do i=1,3
+! write (iout,*) (Im(i,j),j=1,3)
+! enddo
do i=nnt,nct-1
mnum=molnum(i)
Im(3,3)=Im(3,3)+Ip(mnum)*(1-dc_norm(3,i)*dc_norm(3,i))* &
vbld(i+1)*vbld(i+1)*0.25d0
enddo
+! write(iout,*) "The angular momentum before Isc add"
+! do i=1,3
+! write (iout,*) (Im(i,j),j=1,3)
+! enddo
do i=nnt,nct
mnum=molnum(i)
! if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
- .and.(mnum.ne.5)) then
+ .and.(mnum.lt.4)) then
iti=iabs(itype(i,mnum))
inres=i+nres
Im(1,1)=Im(1,1)+Isc(iti,mnum)*(1-dc_norm(1,inres)* &
endif
enddo
+! write(iout,*) "The angular momentum before agnom:"
+! do i=1,3
+! write (iout,*) (Im(i,j),j=1,3)
+! enddo
+
call angmom(cm,L)
! write(iout,*) "The angular momentum before adjustment:"
! write(iout,*) (L(j),j=1,3)
-
+! do i=1,3
+! write (iout,*) (Im(i,j),j=1,3)
+! enddo
Im(2,1)=Im(1,2)
Im(3,1)=Im(1,3)
Im(3,2)=Im(2,3)
do i=nnt,nct-1
call vecpr(vrot(1),dc(1,i),vp)
do j=1,3
+! print *,"HERE2",d_t(j,i),vp(j)
d_t(j,i)=d_t(j,i)-vp(j)
+! print *,"HERE2",d_t(j,i),vp(j)
enddo
enddo
do i=nnt,nct
mnum=molnum(i)
if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
- .and.(mnum.ne.5)) then
+ .and.(mnum.lt.4)) then
! if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
! if(itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
inres=i+nres
enddo
do i=nnt,nct-1
mnum=molnum(i)
- if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
+ if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
do j=1,3
pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
enddo
call vecpr(pr(1),v(1),vp)
do j=1,3
L(j)=L(j)+mp(mnum)*vp(j)
+! print *,"HERE3",J,i,L(j),mp(mnum),Ip(mnum),mnum
enddo
do j=1,3
pr(j)=0.5d0*dc(j,i)
pp(j)=0.5d0*d_t(j,i)
enddo
call vecpr(pr(1),pp(1),vp)
+! print *,vp,"vp"
do j=1,3
L(j)=L(j)+Ip(mnum)*vp(j)
enddo
mnum=molnum(i)
iti=iabs(itype(i,mnum))
inres=i+nres
- if (mnum.eq.5) then
+ if (mnum.gt.4) then
mscab=0.0d0
else
mscab=msc(iti,mnum)
do j=1,3
pr(j)=c(j,inres)-cm(j)
enddo
+ !endif
if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
- .and.(mnum.ne.5)) then
+ .and.(mnum.lt.4)) then
do j=1,3
v(j)=incr(j)+d_t(j,inres)
enddo
v(j)=incr(j)
enddo
endif
+! print *,i,pr,"pr",v
call vecpr(pr(1),v(1),vp)
! write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),&
! " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
L(j)=L(j)+mscab*vp(j)
enddo
! write (iout,*) "L",(l(j),j=1,3)
+! print *,"L",(l(j),j=1,3),i,vp(1)
+
if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
- .and.(mnum.ne.5)) then
+ .and.(mnum.lt.4)) then
do j=1,3
v(j)=incr(j)+d_t(j,inres)
enddo
summas=0.0d0
do i=nnt,nct
mnum=molnum(i)
- if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
+ if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
if (i.lt.nct) then
summas=summas+mp(mnum)
do j=1,3
vcm(j)=vcm(j)+mp(mnum)*(vv(j)+0.5d0*d_t(j,i))
+! print *,i,j,vv(j),d_t(j,i)
enddo
endif
- if (mnum.ne.5) then
+ if (mnum.ne.4) then
amas=msc(iabs(itype(i,mnum)),mnum)
else
amas=0.0d0
endif
summas=summas+amas
if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
- .and.(mnum.ne.5)) then
+ .and.(mnum.lt.4)) then
do j=1,3
vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres))
enddo
enddo
return
end subroutine vcm_vel
+#endif
!-----------------------------------------------------------------------------
! rattle.F
!-----------------------------------------------------------------------------
do i=nnt,nct
mnum=molnum(i)
if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
- .and.(mnum.ne.5)) nbond=nbond+1
+ .and.(mnum.lt.4)) nbond=nbond+1
enddo
! Make a folded form of the Ginv-matrix
ind=0
do k=nnt,nct
mnum=molnum(k)
if (itype(k,1).ne.10 .and. itype(k,mnum).ne.ntyp1_molec(mnum)&
- .and.(mnum.ne.5)) then
+ .and.(mnum.lt.4)) then
jj=jj+1
do l=1,3
ind1=ind1+1
do i=nnt,nct
mnum=molnum(i)
if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
- .and.(mnum.ne.5))
+ .and.(mnum.lt.4))
ii=ii+1
do j=1,3
ind=ind+1
do i=nnt,nct
mnum=molnum(i)
if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
- .and.(mnum.ne.5)) then
+ .and.(mnum.lt.4)) then
ind=ind+1
write (iout,'(2i5,3f10.5,5x,e15.5)') &
do i=nnt,nct
mnum=molnum(i)
if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
- .and.(mnum.ne.5))
+ .and.(mnum.lt.4))
ind=ind+1
xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
do k=nnt,nct
mnum=molnum(i)
if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
- .and.(mnum.ne.5)) then
+ .and.(mnum.lt.4)) then
ind=ind+1
do j=1,3
gdc(j,i,ind)=GGinv(i,ind)*dC(j,k+nres)
do i=nnt,nct
mnum=molnum(i)
if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
- .and.(mnum.ne.5)) then
+ .and.(mnum.lt.4)) then
ind=ind+1
do j=1,nbond
Cmat(ind,j)=0.0d0
do i=nnt,nct
mnum=molnum(i)
if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
- .and.(mnum.ne.5)) then
+ .and.(mnum.lt.4)) then
ind=ind+1
x(ind)=scalar(d_t(1,i+nres),dC(1,i+nres))
endif
do i=nnt,nct
mnum=molnum(i)
if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
- .and.(mnum.ne.5)) then
+ .and.(mnum.lt.4)) then
ind=ind+1
write (iout,'(2i5,3f10.5,5x,e15.5)') &
i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind)
do i=nnt,nct
mnum=molnum(i)
if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
- .and.(mnum.ne.5)) then
+ .and.(mnum.lt.4)) then
ind=ind+1
do j=1,3
xx=0.0d0
do i=nnt,nct
mnum=molnum(i)
if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
- .and.(mnum.ne.5))
+ .and.(mnum.lt.4))
ind=ind+1
write (iout,'(2i5,3f10.5,5x,2e15.5)') &
i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind),&
! include 'COMMON.IOUNITS'
!el real(kind=8),dimension(6*nres) :: gamvec !(MAXRES6) maxres6=6*maxres
!el common /syfek/ gamvec
- real(kind=8) :: vv(3),vvtot(3,nres),v_work(6*nres) !,&
+#ifdef FIVEDIAG
+ integer iposc,ichain,n,innt,inct
+ double precision rs(nres*2)
+ real(kind=8) ::v_work(3,6*nres),vvec(2*nres)
+#else
+ real(kind=8) :: v_work(6*nres)
+#endif
+
+ real(kind=8) :: vv(3),vvtot(3,nres)!,v_work(6*nres) !,&
!el ginvfric(2*nres,2*nres) !maxres2=2*maxres
!el common /przechowalnia/ ginvfric
- logical :: lprn = .false., checkmode = .false.
+ logical :: lprn, checkmode
integer :: i,j,ind,k,nres2,nres6,mnum
nres2=2*nres
nres6=6*nres
+ lprn=.false.
+ checkmode=.false.
+! if (large) lprn=.true.
+! if (large) checkmode=.true.
+#ifdef FIVEDIAG
+!c Here accelerations due to friction forces are computed right after forces.
+ d_t_work(:6*nres)=0.0d0
+ do j=1,3
+ v_work(j,1)=d_t(j,0)
+ v_work(j,nnt)=d_t(j,0)
+ enddo
+ do i=nnt+1,nct
+ do j=1,3
+ v_work(j,i)=v_work(j,i-1)+d_t(j,i-1)
+ enddo
+ enddo
+ do i=nnt,nct
+ mnum=molnum(i)
+ if (iabs(itype(i,1)).ne.10 .and. iabs(itype(i,mnum)).ne.ntyp1_molec(mnum).and.mnum.lt.3) then
+ do j=1,3
+ v_work(j,i+nres)=v_work(j,i)+d_t(j,i+nres)
+ enddo
+ endif
+ enddo
+#ifdef DEBUG
+ write (iout,*) "v_work"
+ do i=1,2*nres
+ write (iout,'(i5,3f10.5)') i,(v_work(j,i),j=1,3)
+ enddo
+#endif
+ do j=1,3
+ ind=0
+ do ichain=1,nchain
+ n=dimen_chain(ichain)
+ iposc=iposd_chain(ichain)
+!c write (iout,*) "friction_force j",j," ichain",ichain,
+!c & " n",n," iposc",iposc,iposc+n-1
+ innt=chain_border(1,ichain)
+ inct=chain_border(2,ichain)
+!c diagnostics
+!c innt=chain_border(1,1)
+!c inct=chain_border(2,1)
+ do i=innt,inct
+ vvec(ind+1)=v_work(j,i)
+ ind=ind+1
+! if (iabs(itype(i)).ne.10) then
+ if (iabs(itype(i,1)).ne.10 .and. iabs(itype(i,mnum)).ne.ntyp1_molec(mnum).and.mnum.lt.3) then
+ vvec(ind+1)=v_work(j,i+nres)
+ ind=ind+1
+ endif
+ enddo
+#ifdef DEBUG
+ write (iout,*) "vvec ind",ind," n",n
+ write (iout,'(f10.5)') (vvec(i),i=iposc,ind)
+#endif
+!c write (iout,*) "chain",i," ind",ind," n",n
+#ifdef TIMING
+#ifdef MPI
+ time01=MPI_Wtime()
+#else
+ time01=tcpu()
+#endif
+#endif
+! if (large) print *,"before fivediagmult"
+ call fivediagmult(n,DMfric(iposc),DU1fric(iposc),&
+ DU2fric(iposc),vvec(iposc),rs)
+! if (large) print *,"after fivediagmult"
+#ifdef TIMING
+#ifdef MPI
+ time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
+#else
+ time_fricmatmult=time_fricmatmult+tcpu()-time01
+#endif
+#endif
+#ifdef DEBUG
+ write (iout,*) "rs"
+ write (iout,'(f10.5)') (rs(i),i=1,n)
+#endif
+ do i=iposc,iposc+n-1
+! write (iout,*) "ichain",ichain," i",i," j",j,&
+! "index",3*(i-1)+j,"rs",rs(i-iposc+1)
+ fric_work(3*(i-1)+j)=-rs(i-iposc+1)
+ enddo
+ enddo
+ enddo
+#ifdef DEBUG
+ write (iout,*) "Vector fric_work dimen3",dimen3
+ write (iout,'(3f10.5)') (fric_work(j),j=1,dimen3)
+#endif
+#else
if(.not.allocated(gamvec)) allocate(gamvec(nres6)) !(MAXRES6)
if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
do i=0,nres2
do i=nnt,nct
mnum=molnum(i)
if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
- .and.(mnum.ne.5)) then
+ .and.(mnum.lt.4)) then
do j=1,3
d_t_work(ind+j)=d_t(j,i+nres)
enddo
do i=nnt,nct
mnum=molnum(i)
if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
- .and.(mnum.ne.5)) then
+ .and.(mnum.lt.4)) then
! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
do j=1,3
friction(j,i+nres)=fric_work(ind+j)
enddo
enddo
endif
+#endif
return
end subroutine friction_force
!-----------------------------------------------------------------------------
!#endif
! include 'COMMON.IOUNITS'
integer :: IERROR
- integer :: i,j,ind,ind1,m
+ integer :: i,j,ind,ind1,m,ichain,innt,inct
logical :: lprn = .false.
real(kind=8) :: dtdi !el ,gamvec(2*nres)
!el real(kind=8),dimension(2*nres,2*nres) :: ginvfric,fcopy
nres2=2*nres
nres6=6*nres
#ifdef MPI
+#ifndef FIVEDIAG
if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
if (fg_rank.ne.king) goto 10
#endif
+#endif
! nres2=2*nres
! nres6=6*nres
if(.not.allocated(gamvec)) allocate(gamvec(nres2)) !(MAXRES2)
+#ifndef FIVEDIAG
if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
!el allocate(fcopy(nres2,nres2)) !maxres2=2*maxres
if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !maxres2=2*maxres
if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
+#endif
+#ifdef FIVEDIAG
+ if (.not.allocated(DMfric)) allocate(DMfric(nres2))
+ if (.not.allocated(DU1fric)) allocate(DU1fric(nres2))
+ if (.not.allocated(DU2fric)) allocate(DU2fric(nres2))
+! Load the friction coefficients corresponding to peptide groups
+ ind1=0
+ do i=nnt,nct-1
+ mnum=molnum(i)
+ ind1=ind1+1
+ gamvec(ind1)=gamp(mnum)
+ enddo
+!HERE TEST
+ if (molnum(nct).eq.5) then
+ mnum=molnum(i)
+ ind1=ind1+1
+ gamvec(ind1)=gamp(mnum)
+ endif
+! Load the friction coefficients corresponding to side chains
+ m=nct-nnt
+ ind=0
+ do j=1,2
+ gamsc(ntyp1_molec(j),j)=1.0d0
+ enddo
+ do i=nnt,nct
+ mnum=molnum(i)
+ ind=ind+1
+ ii = ind+m
+ iti=itype(i,mnum)
+ gamvec(ii)=gamsc(iabs(iti),mnum)
+ enddo
+ if (surfarea) call sdarea(gamvec)
+ DMfric=0.0d0
+ DU1fric=0.0d0
+ DU2fric=0.0d0
+ ind=1
+ do ichain=1,nchain
+ innt=chain_border(1,ichain)
+ inct=chain_border(2,ichain)
+!c write (iout,*) "ichain",ichain," innt",innt," inct",inct
+!c DMfric part
+ mnum=molnum(innt)
+ DMfric(ind)=gamvec(innt-nnt+1)/4
+ if (iabs(itype(innt,1)).eq.10.or.mnum.gt.2) then
+ DMfric(ind)=DMfric(ind)+gamvec(m+innt-nnt+1)
+ ind=ind+1
+ else
+ DMfric(ind+1)=gamvec(m+innt-nnt+1)
+ ind=ind+2
+ endif
+!c write (iout,*) "DMfric init ind",ind
+!c DMfric
+ do i=innt+1,inct-1
+ mnum=molnum(i)
+ DMfric(ind)=gamvec(i-nnt+1)/2
+ if (iabs(itype(i,1)).eq.10.or.mnum.gt.2) then
+ DMfric(ind)=DMfric(ind)+gamvec(m+i-nnt+1)
+ ind=ind+1
+ else
+ DMfric(ind+1)=gamvec(m+i-nnt+1)
+ ind=ind+2
+ endif
+ enddo
+!c write (iout,*) "DMfric endloop ind",ind
+ if (inct.gt.innt) then
+ DMfric(ind)=gamvec(inct-1-nnt+1)/4
+ mnum=molnum(inct)
+ if (iabs(itype(inct,1)).eq.10.or.mnum.gt.2) then
+ DMfric(ind)=DMfric(ind)+gamvec(inct+m-nnt+1)
+ ind=ind+1
+ else
+ DMfric(ind+1)=gamvec(inct+m-nnt+1)
+ ind=ind+2
+ endif
+ endif
+!c write (iout,*) "DMfric end ind",ind
+ enddo
+!c DU1fric part
+ do ichain=1,nchain
+ mnum=molnum(i)
+
+ ind=iposd_chain(ichain)
+ innt=chain_border(1,ichain)
+ inct=chain_border(2,ichain)
+ do i=innt,inct
+ if (iabs(itype(i,1)).ne.10.and.mnum.le.2) then
+ ind=ind+2
+ else
+ DU1fric(ind)=gamvec(i-nnt+1)/4
+ ind=ind+1
+ endif
+ enddo
+ enddo
+!c DU2fric part
+ do ichain=1,nchain
+ mnum=molnum(i)
+ ind=iposd_chain(ichain)
+ innt=chain_border(1,ichain)
+ inct=chain_border(2,ichain)
+ do i=innt,inct-1
+ if (iabs(itype(i,1)).ne.10.and.mnum.le.2) then
+ DU2fric(ind)=gamvec(i-nnt+1)/4
+ DU2fric(ind+1)=0.0d0
+ ind=ind+2
+ else
+ DU2fric(ind)=0.0d0
+ ind=ind+1
+ endif
+ enddo
+ enddo
+ if (lprn) then
+ write(iout,*)"The upper part of the five-diagonal friction matrix"
+ do ichain=1,nchain
+ write (iout,'(a,i5)') 'Chain',ichain
+ innt=iposd_chain(ichain)
+ inct=iposd_chain(ichain)+dimen_chain(ichain)-1
+ do i=innt,inct
+ if (i.lt.inct-1) then
+ write (iout,'(2i3,3f10.5)') i,i-innt+1,DMfric(i),DU1fric(i),&
+ DU2fric(i)
+ else if (i.eq.inct-1) then
+ write (iout,'(2i3,3f10.5)') i,i-innt+1,DMfric(i),DU1fric(i)
+ else
+ write (iout,'(2i3,3f10.5)') i,i-innt+1,DMfric(i)
+ endif
+ enddo
+ enddo
+ endif
+ 10 continue
+#else
+
+
! Zeroing out fricmat
do i=1,dimen
do j=1,dimen
ind1=ind1+1
gamvec(ind1)=gamp(mnum)
enddo
+!HERE TEST
+ if (molnum(nct).eq.5) then
+ mnum=molnum(i)
+ ind1=ind1+1
+ gamvec(ind1)=gamp(mnum)
+ endif
! Load the friction coefficients corresponding to side chains
m=nct-nnt
ind=0
! call MATOUT2(my_ng_count,dimen,maxres2,maxres2,fcopy)
endif
#endif
+#endif
return
end subroutine setup_fricmat
!-----------------------------------------------------------------------------
real(kind=8) :: time00
logical :: lprn = .false.
integer :: i,j,ind,mnum
+#ifdef FIVEDIAG
+ integer ichain,innt,inct,iposc
+#endif
do i=0,2*nres
do j=1,3
#else
time_fsample=time_fsample+tcpu()-time00
#endif
+#ifdef FIVEDIAG
+ ind=0
+ do ichain=1,nchain
+ innt=chain_border(1,ichain)
+ inct=chain_border(2,ichain)
+ iposc=iposd_chain(ichain)
+!c for debugging only
+!c innt=chain_border(1,1)
+!c inct=chain_border(2,1)
+!c iposc=iposd_chain(1)
+!c write (iout,*)"stochastic_force ichain=",ichain," innt",innt,
+!c & " inct",inct," iposc",iposc
+ do j=1,3
+ stochforcvec(ind+j)=0.5d0*force(j,innt)
+ enddo
+ if (iabs(itype(innt,molnum(innt))).eq.10.or.molnum(innt).ge.3) then
+ do j=1,3
+ stochforcvec(ind+j)=stochforcvec(ind+j)+force(j,innt+nres)
+ enddo
+ ind=ind+3
+ else
+ ind=ind+3
+ do j=1,3
+ stochforcvec(ind+j)=force(j,innt+nres)
+ enddo
+ ind=ind+3
+ endif
+ do i=innt+1,inct-1
+ do j=1,3
+ stochforcvec(ind+j)=0.5d0*(force(j,i)+force(j,i-1))
+ enddo
+ if (iabs(itype(i,1).eq.10).or.molnum(i).ge.3) then
+ do j=1,3
+ stochforcvec(ind+j)=stochforcvec(ind+j)+force(j,i+nres)
+ enddo
+ ind=ind+3
+ else
+ ind=ind+3
+ do j=1,3
+ stochforcvec(ind+j)=force(j,i+nres)
+ enddo
+ ind=ind+3
+ endif
+ enddo
+ do j=1,3
+ stochforcvec(ind+j)=0.5d0*force(j,inct-1)
+ enddo
+ if (iabs(itype(inct,molnum(inct))).eq.10.or.molnum(inct).ge.3) then
+ do j=1,3
+ stochforcvec(ind+j)=stochforcvec(ind+j)+force(j,inct+nres)
+ enddo
+ ind=ind+3
+ else
+ ind=ind+3
+ do j=1,3
+ stochforcvec(ind+j)=force(j,inct+nres)
+ enddo
+ ind=ind+3
+ endif
+!c write (iout,*) "chain",ichain," ind",ind
+ enddo
+#ifdef DEBUG
+ write (iout,*) "stochforcvec"
+ write (iout,'(3f10.5)') (stochforcvec(j),j=1,ind)
+#endif
+#else
! Compute the stochastic forces acting on virtual-bond vectors.
do j=1,3
ff(j)=0.0d0
do i=nnt,nct
mnum=molnum(i)
if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
- .and.(mnum.ne.5)) then
+ .and.(mnum.lt.4)) then
! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
do j=1,3
stochforc(j,i+nres)=force(j,i+nres)
do i=nnt,nct
mnum=molnum(i)
if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
- .and.(mnum.ne.5)) then
+ .and.(mnum.lt.4)) then
! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
do j=1,3
stochforcvec(ind+j)=stochforc(j,i+nres)
enddo
endif
-
+#endif
return
end subroutine stochastic_force
!-----------------------------------------------------------------------------
! commom.langevin.lang0
! common /langforc/
allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
+#ifndef FIVEDIAG
if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
allocate(fricvec(nres2,nres2)) !(MAXRES2,MAXRES2)
+#endif
allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
#endif
-
+#ifndef FIVEDIAG
if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2))
+#endif
!----------------------
! commom.hairpin in CSA module
!----------------------
#ifdef FIVEDIAG
allocate(DM(nres2),DU1(nres2),DU2(nres2))
allocate(DMorig(nres2),DU1orig(nres2),DU2orig(nres2))
+!#ifdef DEBUG
+ allocate(Gvec(nres2,nres2))
+!#endif
#else
write (iout,*) "Before A Allocation",nfgtasks-1
call flush(iout)
projects / unres4.git / blobdiff