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
+ implicit none
+ include 'DIMENSIONS'
+ include 'COMMON.VAR'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.GEO'
+ include 'COMMON.LOCAL'
+ include 'COMMON.INTERACT'
+ include 'COMMON.MD'
+ include 'COMMON.LAGRANGE.5diag'
+ include 'COMMON.IOUNITS'
+ double precision KE_total
+ integer i,j,k,iti
+ double precision KEt_p,KEt_sc,KEr_p,KEr_sc,incr(3),&
+ mag1,mag2,v(3)
+
+ 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
+!c write (iout,*) "Kinetic trp:",i,(incr(j),j=1,3
+!c Skip dummy peptide groups
+ if (itype(i).ne.ntyp1 .and. itype(i+1).ne.ntyp1) then
+ do j=1,3
+ v(j)=incr(j)+0.5d0*d_t(j,i)
+ enddo
+!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+(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
+ iti=iabs(itype(i))
+ if (itype(i).eq.10 .and. itype(i).ne.ntyp1) 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
+!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)*(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
+ 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
+ if (itype(i).ne.ntyp1.and.itype(i+1).ne.ntyp1) then
+! 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+(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
+ iti=iabs(itype(i))
+ if (itype(i).ne.10.and.itype(i).ne.ntyp1) then
+ do j=1,3
+ incr(j)=d_t(j,nres+i)
+ enddo
+!c write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3)
+ KEr_sc=KEr_sc+Isc(iti)*(incr(1)*incr(1)+incr(2)*incr(2)+&
+ incr(3)*incr(3))
+ endif
+ enddo
+!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*(mp*KEt_p+KEt_sc+0.25d0*Ip*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
+ 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
+!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+(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
+ iti=iabs(itype(i))
+ if (iti.eq.10) 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)*(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
+ do j=1,3
+ incr(j)=d_t(j,i+1)-d_t(j,i)
+ enddo
+!c write (iout,*) i,(incr(j),j=1,3)
+!c write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3)
+ KEr_p=KEr_p+(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
+ iti=iabs(itype(i))
+ if (iti.ne.10) 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)*(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*(mp*KEt_p+KEt_sc+0.25d0*Ip*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
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,1).ne.ntyp1.and.mnum.lt.4) then
do j=1,3
ind=ind+1
v_work(ind)=d_t(j,i+nres)
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')
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
! 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)
+#else
call ginv_mult(fric_work, d_af_work)
call ginv_mult(stochforcvec, d_as_work)
+#endif
+
return
end subroutine sddir_precalc
!-----------------------------------------------------------------------------
! 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
! Compute the acceleration due to friction forces (d_af_work) and stochastic
! forces (d_as_work)
!
+#ifdef FIVEDIAG
+ call fivediaginv_mult(maxres6,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)
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"
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
! include 'COMMON.NAMES'
! include 'COMMON.TIME1'
real(kind=8) :: xv,sigv,lowb,highb ,Ek1
+#ifdef FIVEDIAG
+ integer ichain,n,innt,inct,ibeg,ierr
+ double precision work(48*nres)
+ integer iwork(maxres6)
+ 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
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)
+ if (mnum.ge.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
+ .or.(mnum.ge.4)) then
j=ii+3
else
j=ii+6
enddo
endif
if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
- .and.(mnum.ne.5)) ii=ii+3
+ .and.(mnum.lt.4)) ii=ii+3
write (iout,*) "i",i," itype",itype(i,mnum)," mass",msc(itype(i,mnum),mnum)
write (iout,*) "ii",ii
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))&
+ .and.(mnum.lt.4))&
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
enddo
enddo
mnum=molnum(j)
if (itype(j,1).ne.10 .and. itype(j,mnum).ne.ntyp1_molec(mnum)&
- .and.(mnum.ne.5)) then
+ .and.(mnum.lt.4)) then
do k=1,3
d_t(k,j+nres)=d_t_work(ind)
ind=ind+1
! 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
+ .or.(mnum.ge.4)) then
do j=1,3
d_t(j,i)=d_t(j,i+1)-d_t(j,i)
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
do j=1,3
ind=ind+1
d_t(j,i+nres)=d_t_work(ind)
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)
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))
! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
inres=i+nres
do j=1,3
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
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
+#ifdef FIVEDIAG
+ integer iposc,ichain,n,innt,inct
+ double precision rs(nres+2)
+#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
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
+ call fivediagmult(n,DMfric(iposc),DU1fric(iposc),
+ & DU2fric(iposc),vvec(iposc),rs)
+#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
+c write (iout,*) "ichain",ichain," i",i," j",j,
+c & "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
!-----------------------------------------------------------------------------
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
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(iint)).eq.10) 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,molnum(i)).eq.10) 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) 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
!-----------------------------------------------------------------------------
projects / unres4.git / blobdiff