real(kind=8),dimension(:,:),allocatable :: ginvfric !(2*nres,2*nres) !maxres2=2*maxres
!-----------------------------------------------------------------------------
! common /przechowalnia/ subroutine: setup_fricmat
+#ifndef LBFGS
real(kind=8),dimension(:,:),allocatable :: fcopy !(2*nres,2*nres)
+#endif
!-----------------------------------------------------------------------------
!
!
do j=1,3
incr(j)=d_t(j,0)
enddo
- do i=nnt,nct-1
+ 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
do j=1,3
v(j)=incr(j)+0.5d0*d_t(j,i)
enddo
- if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
+ 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))
do i=nnt,nct
mnum=molnum(i)
iti=iabs(itype(i,mnum))
+ if (mnum.eq.5) iti=itype(i,mnum)
+! if (itype(i,mnum).eq.ntyp1_molec(mnum)) then
+! do j=1,3
+! incr(j)=d_t(j,i)
+! enddo
+! endif
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)+d_t(j,nres+i)
enddo
endif
- write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3)
- write (iout,*) "i",i," msc",msc(iti,1)," v",(v(j),j=1,3)
+! 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
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)
+ if (mnum.eq.5) Ip(mnum)=0.0
! write (iout,*) "i",i
! write (iout,*) "i",i," mag1",mag1," mag2",mag2
do j=1,3
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
+! 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
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)=0.0d0
+! 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))
!c to the velocities of the first Calpha.
do i=innt,inct
mnum=molnum(i)
+ if (mnum.eq.5) then
+ iti=itype(i,mnum)
+ else
iti=iabs(itype(i,mnum))
+ endif
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
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)
+ if (mnum.eq.5) Ip(mnum)=0.0d0
!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)&
!c The rotational part of the side chain virtual bond
do i=innt,inct
mnum=molnum(i)
- iti=iabs(itype(i,mnum))
+! iti=iabs(itype(i,mnum))
+ if (mnum.eq.5) then
+ iti=itype(i,mnum)
+ else
+ iti=iabs(itype(i,mnum))
+ endif
+
! 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
! include 'COMMON.NAMES'
! include 'COMMON.TIME1'
! include 'COMMON.HAIRPIN'
- real(kind=8),dimension(3) :: L,vcm
+ real(kind=8),dimension(3) :: L,vcm,boxx
#ifdef VOUT
real(kind=8),dimension(6*nres) :: v_work,v_transf !(maxres6) maxres6=6*maxres
#endif
!el common /gucio/ cm
integer :: i,j,nharp
integer,dimension(4,nres) :: iharp !(4,nres/3)(4,maxres/3)
+
! logical :: ovrtim
real(kind=8) :: tt0,scalfac
integer :: nres2,itime
nres2=2*nres
print *, "ENTER MD"
+ boxx(1)=boxxsize
+ boxx(2)=boxysize
+ boxx(3)=boxzsize
+
!
#ifdef MPI
print *,"MY tmpdir",tmpdir,ilen(tmpdir)
endif
if (reset_vel .and. tbf .and. lang.eq.0 &
.and. mod(itime,count_reset_vel).eq.0) then
+ !WARP WATER
+ do i=1,nres
+ if (molnum(i).eq.5) then
+ call to_box(c(1,i),c(2,i),c(3,i))
+ do j=1,3
+ if (c(j,i).le.0) c(j,i)=c(j,i)+boxx(j)
+ enddo
+ endif
+! write(iout,*) "COORD",c(1,i),c(2,i),c(3,i)
+ enddo
call random_vel
+
write(iout,'(a,f20.2)') &
"Velocities reset to random values, time",totT
do i=0,2*nres
endif
if (mod(itime,ntwx).eq.0) then
call returnbox
+ call enerprint(potEcomp)
+
write (tytul,'("time",f8.2)') totT
if(mdpdb) then
+ write(iout,*) "before hairpin"
call hairpin(.true.,nharp,iharp)
+ write(iout,*) "before secondary"
call secondary2(.true.)
+ write(iout,*) "before pdbout"
call pdbout(potE,tytul,ipdb)
- call enerprint(potEcomp)
+! call enerprint(potEcomp)
else
call cartout(totT)
endif
! call ginv_mult(fric_work, d_af_work)
! call ginv_mult(stochforcvec, d_as_work)
#ifdef FIVEDIAG
+ write(iout,*) "forces before fivediaginv"
+ do i=1,dimen*3
+ write(iout,*) "fricwork",i,fric_work(i)
+ enddo
call fivediaginv_mult(dimen,fric_work, d_af_work)
call fivediaginv_mult(dimen,stochforcvec, d_as_work)
if (large) then
! include 'COMMON.INTERACT'
! include 'COMMON.IOUNITS'
! include 'COMMON.NAMES'
- real(kind=8),dimension(6*nres) :: stochforcvec,d_as_work1 !(MAXRES6) maxres6=6*maxres
+ real(kind=8),dimension(:),allocatable :: stochforcvec,d_as_work1 !(MAXRES6) maxres6=6*maxres
real(kind=8) :: cos60 = 0.5d0, sin60 = 0.86602540378443864676d0
integer :: i,j,ind,inres,mnum
+ if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres))
+ if (.not.allocated(d_as_work1)) allocate(d_as_work1(6*nres))
! Revised 3/31/05 AL: correlation between random contributions to
! position and velocity increments included.
! The correlation coefficients are calculated at low-friction limit.
character(len=50) :: tytul
logical :: file_exist
!el common /gucio/ cm
- integer :: i,j,ipos,iq,iw,nft_sc,iretcode,nfun,ierr,mnum,itime
+ integer :: i,j,ipos,iq,iw,nft_sc,iretcode,ierr,mnum,itime
+#ifndef LBFGS
+ integer :: nfun
+#endif
real(kind=8) :: etot,tt0
logical :: fail
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
call inertia_tensor
! Getting the potential energy and forces and velocities and accelerations
call vcm_vel(vcm)
-! write (iout,*) "velocity of the center of the mass:"
-! write (iout,*) (vcm(j),j=1,3)
+ write (iout,*) "velocity of the center of the mass:"
+ write (iout,*) (vcm(j),j=1,3)
do j=1,3
d_t(j,0)=d_t(j,0)-vcm(j)
enddo
if (me.eq.king.or..not.out1file) write(iout,*)&
'Minimizing initial PDB structure: Calling MINIM_DC'
call minim_dc(etot,iretcode,nfun)
+#ifdef LBFGS
+ if (me.eq.king.or..not.out1file)&
+ write(iout,*) 'LBFGS return code is ',statusbf,' eval ',nfun
+#endif
else
call geom_to_var(nvar,varia)
if(me.eq.king.or..not.out1file) write (iout,*)&
call var_to_geom(nvar,varia)
#ifdef LBFGS
if (me.eq.king.or..not.out1file)&
- write(iout,*) 'LBFGS return code is ',status,' eval ',nfun
+ write(iout,*) 'LBFGS return code is ',statusbf,' eval ',nfun
if(me.eq.king.or..not.out1file)&
- write(iout,*) 'LBFGS return code is ',status,' eval ',nfun
+ write(iout,*) 'LBFGS return code is ',statusbf,' eval ',nfun
#else
if (me.eq.king.or..not.out1file)&
write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
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.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(6*nres)
+ integer ichain,n,innt,inct,ibeg,ierr,innt_org
+ 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
+!#ifdef DEBUG
+! double precision inertia(maxres2_chain,maxres2_chain)
+!#endif
#endif
-#define DEBUG
+!#define DEBUG
#ifdef FIVEDIAG
- real(kind=8) ,allocatable, dimension(:) :: DDU1,DDU2,DL2,DL1,xsolv,DML,rs
- real(kind=8) :: sumx,Ek2,Ek3,aux
+ real(kind=8) ,allocatable, dimension(:) :: xsolv,DML,rs
+ real(kind=8) :: sumx,Ek2,Ek3,aux,masinv
#ifdef DEBUG
real(kind=8) ,allocatable, dimension(:) :: rsold
real (kind=8),allocatable,dimension(:,:) :: matold,inertia
integer :: iti
#endif
- integer :: i,j,ii,k,ind,mark,imark,mnum
+#endif
+ 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)
-!#define DEBUG
+#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,*) "Random_vel, fivediag"
+ flush(iout)
allocate(inertia(2*nres,2*nres))
#endif
d_t=0.0d0
EK=0.0d0
Ek3=0.0d0
#ifdef DEBUG
- write(iout,*), nchain
+ write(iout,*), "nchain",nchain
#endif
do ichain=1,nchain
ind=0
- ghalf=0.0d0
+! 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)
+! innt_org=
+ innt_org=chain_border(1,ichain)
+ if ((molnum(innt_org).eq.5).or.(molnum(innt_org).eq.4)) go to 137
+ if(.not.allocated(ghalf)) print *,"COCO"
+ if(.not.allocated(Ghalf)) allocate(Ghalf(1300*(1300+1)/2))
+ ghalf=0.0d0
inct=innt+n-1
#ifdef DEBUG
write (iout,*) "Chain",ichain," n",n," start",innt
enddo
enddo
#endif
+137 continue
+ write(iout,*) "HERE,",n,innt
+ innt_org=chain_border(1,ichain)
xv=0.0d0
ii=0
do i=1,n
do k=1,3
ii=ii+1
+ mnum=molnum(innt_org)
+ if (molnum(innt_org).ge.4) geigen(i)=3.0/msc(itype(innt_org+i-1,mnum),mnum)
+! if (molnum(innt).eq.5) write(iout,*) "typ",i,innt-1+i,itype(innt+i-1,5)
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)
+ write (iout,*) "i",i," ii",ii," geigen",geigen(i), &
+ " d_t_work_new",d_t_work_new(ii),innt_org+i-1
+ enddo
+ enddo
+ if (molnum(innt_org).ge.4) then
+ mnum=molnum(innt_org)
+ do k=1,3
+ do i=1,n
+ ind=(i-1)*3+k
+ d_t_work(ind)=0.0d0
+ masinv=1.0d0/msc(itype(innt_org+i-1,mnum),mnum)
+ d_t_work(ind)=d_t_work(ind)&
+ +masinv*d_t_work_new((i-1)*3+k)
enddo
enddo
+
+ else
do k=1,3
do i=1,n
ind=(i-1)*3+k
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
+ endif
#ifdef DEBUG
aux=0.0d0
do k=1,3
d_t(:,0)=d_t(:,1)
d_t(:,1)=0.0d0
endif
+ if (large) then
+ write (iout,*)
+ write (iout,*) "Random vel after 1st transf the Calpha,SC space"
+ write (iout,'(3hORG,1h(,i5,1h),3f10.5)') 0,(d_t(j,0),j=1,3)
+ 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
+
!c d_a(:,0)=d_a(:,1)
!c d_a(:,1)=0.0d0
!c write (iout,*) "Shifting accelerations"
do ichain=2,nchain
+ write(iout,*) "nchain",ichain,chain_border1(1,ichain),molnum(chain_border1(1,ichain))
+ if (molnum(chain_border1(1,ichain)+1).eq.5) cycle
!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))
enddo
!c write (iout,*) "Adding accelerations"
do ichain=2,nchain
+ if (molnum(chain_border1(1,ichain)+1).eq.5) cycle
!c write (iout,*) "chain",ichain,chain_border1(1,ichain)-1,
!c & chain_border(2,ichain-1)
d_t(:,chain_border1(1,ichain)-1)=&
do ichain=2,nchain
write (iout,*) "chain",ichain,chain_border1(1,ichain)-1,&
chain_border(2,ichain-1)
+ if (molnum(chain_border1(1,ichain)+1).eq.5) cycle
+
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
+ if (large) then
+ write (iout,*)
+ write (iout,*) "Random vel after 2nd transf the Calpha,SC space"
+ write (iout,'(3hORG,1h(,i5,1h),3f10.5)') 0,(d_t(j,0),j=1,3)
+ 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
+
#else
ibeg=0
!c do j=1,3
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)=msc(itype(i,mnum),mnum)
+ if (mnum.ge.5) mp(mnum)=0.0d0
M_PEP=M_PEP+mp(mnum)
do j=1,3
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)=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
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)=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
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
+! if (mnum.gt.4) then
+! mscab=0.0d0
+! else
mscab=msc(iti,mnum)
- endif
+! endif
do j=1,3
L(j)=L(j)+mscab*vp(j)
enddo
summas=0.0d0
do i=nnt,nct
mnum=molnum(i)
- if ((mnum.ge.5).or.(mnum.eq.3))&
+ if ((mnum.gt.5).or.(mnum.eq.3))&
mp(mnum)=msc(itype(i,mnum),mnum)
if (i.lt.nct) then
summas=summas+mp(mnum)
vcm(j)=vcm(j)+mp(mnum)*(vv(j)+0.5d0*d_t(j,i))
enddo
endif
- if (mnum.ne.4) then
+! if (mnum.ne.4) then
amas=msc(iabs(itype(i,mnum)),mnum)
- else
- amas=0.0d0
- endif
+! else
+! amas=0.0d0
+! endif
! amas=msc(iabs(itype(i)))
summas=summas+amas
! if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
vv(j)=vv(j)+d_t(j,i)
enddo
enddo
-!c write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas
+ write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas
do j=1,3
vcm(j)=vcm(j)/summas
enddo
M_PEP=0.0d0
do i=nnt,nct-1
mnum=molnum(i)
- if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
+! 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)
M_SC=0.0d0
do i=nnt,nct
mnum=molnum(i)
- if (mnum.ge.5) cycle
+! if (mnum.ge.5) cycle
iti=iabs(itype(i,mnum))
M_SC=M_SC+msc(iabs(iti),mnum)
inres=i+nres
+ if (mnum.ge.4) inres=i
do j=1,3
cm(j)=cm(j)+msc(iabs(iti),mnum)*c(j,inres)
enddo
! write(iout,*) "Center of mass:",cm
do i=nnt,nct-1
mnum=molnum(i)
- if (mnum.ge.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
do i=nnt,nct
mnum=molnum(i)
iti=iabs(itype(i,mnum))
- if (mnum.ge.5) cycle
+! if (mnum.ge.5) cycle
inres=i+nres
+ if (mnum.ge.4) inres=i
do j=1,3
pr(j)=c(j,inres)-cm(j)
enddo
summas=0.0d0
do i=nnt,nct
mnum=molnum(i)
- if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
+! if (mnum.ge.4) mp(mnum)=msc(itype(i,mnum),mnum)
if (i.lt.nct) then
summas=summas+mp(mnum)
do j=1,3
! print *,i,j,vv(j),d_t(j,i)
enddo
endif
- if (mnum.ne.4) then
+! if (mnum.ne.4) then
amas=msc(iabs(itype(i,mnum)),mnum)
- else
- amas=0.0d0
- endif
+! else
+! amas=0.0d0
+! endif
summas=summas+amas
if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
.and.(mnum.lt.4)) then
vv(j)=vv(j)+d_t(j,i)
enddo
enddo
-! write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas
+ write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas
do j=1,3
vcm(j)=vcm(j)/summas
enddo
!c innt=chain_border(1,1)
!c inct=chain_border(2,1)
do i=innt,inct
+ mnum=molnum(i)
vvec(ind+1)=v_work(j,i)
ind=ind+1
! if (iabs(itype(i)).ne.10) then
! include 'COMMON.IOUNITS'
integer :: IERROR
integer :: i,j,ind,ind1,m,ichain,innt,inct
- logical :: lprn = .false.
+ logical :: lprn = .true.
real(kind=8) :: dtdi !el ,gamvec(2*nres)
!el real(kind=8),dimension(2*nres,2*nres) :: ginvfric,fcopy
! real(kind=8),allocatable,dimension(:,:) :: fcopy
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
+ mnum=molnum(i)
if (iabs(itype(i,1)).ne.10.and.mnum.le.2) then
ind=ind+2
else
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
+ mnum=molnum(i)
if (iabs(itype(i,1)).ne.10.and.mnum.le.2) then
DU2fric(ind)=gamvec(i-nnt+1)/4
DU2fric(ind+1)=0.0d0
! 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))
- allocate(Gvec(nres2,nres2))
+!#ifdef DEBUG
+ allocate(Gvec(1300,1300))!maximum number of protein in one chain
+!#endif
#else
write (iout,*) "Before A Allocation",nfgtasks-1
call flush(iout)
projects / unres4.git / blobdiff