real(kind=8),dimension(:,:),allocatable :: ginvfric !(2*nres,2*nres) !maxres2=2*maxres
!-----------------------------------------------------------------------------
! common /przechowalnia/ subroutine: setup_fricmat
-!el real(kind=8),dimension(:,:),allocatable :: fcopy !(2*nres,2*nres)
+ real(kind=8),dimension(:,:),allocatable :: fcopy !(2*nres,2*nres)
!-----------------------------------------------------------------------------
!
!
real(kind=8),dimension(6*nres,2*nres) :: Bmat,GBmat,Tmat !(MAXRES6,MAXRES2) (maxres2=2*maxres,maxres6=6*maxres)
real(kind=8),dimension(2*nres,2*nres) :: Cmat_,Cinv !(maxres2,maxres2) maxres2=2*maxres
real(kind=8),dimension(6*nres,6*nres) :: Pmat !(maxres6,maxres6) maxres6=6*maxres
+! real(kind=8),dimension(:,:),allocatable :: Bmat,GBmat,Tmat !(MAXRES6,MAXRES2) (maxres2=2*maxres,maxres6=6*maxres)
+! real(kind=8),dimension(:,:),allocatable :: Cmat_,Cinv !(maxres2,maxres2) maxres2=2*maxres
+! real(kind=8),dimension(:,:),allocatable :: Pmat !(maxres6,maxres6) maxres6=6*maxres
real(kind=8),dimension(6*nres) :: Td !(maxres6) maxres6=6*maxres
real(kind=8),dimension(2*nres) :: ppvec !(maxres2) maxres2=2*maxres
!el common /stochcalc/ stochforcvec
logical :: osob
nres2=2*nres
nres6=6*nres
+! if (.not.allocated(Bmat)) allocate(Bmat(nres6,nres2))
+! if (.not.allocated(GBmat)) allocate (GBmat(nres6,nres2))
+! if (.not.allocated(Tmat)) allocate (Tmat(nres6,nres2))
+! if (.not.allocated(Cmat_)) allocate(Cmat_(nres2,nres2))
+! if (.not.allocated(Cinv)) allocate (Cinv(nres2,nres2))
+! if (.not.allocated(Pmat)) allocate(Pmat(6*nres,6*nres))
if (.not.allocated(stochforcvec)) allocate(stochforcvec(nres6)) !(MAXRES6) maxres6=6*maxres
! include 'COMMON.INTERACT'
! include 'COMMON.MD'
! include 'COMMON.IOUNITS'
- real(kind=8) :: KE_total
+ real(kind=8) :: KE_total,mscab
integer :: i,j,k,iti,mnum
real(kind=8) :: KEt_p,KEt_sc,KEr_p,KEr_sc,incr(3),&
enddo
do i=nnt,nct-1
mnum=molnum(i)
+ if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
! write (iout,*) "Kinetic trp:",i,(incr(j),j=1,3)
do j=1,3
v(j)=incr(j)+0.5d0*d_t(j,i)
do i=nnt,nct
mnum=molnum(i)
iti=iabs(itype(i,mnum))
+ if (mnum.eq.5) then
+ mscab=0.0d0
+ else
+ mscab=msc(iti,mnum)
+ endif
! 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
+ if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
+ .or.(mnum.eq.5)) then
do j=1,3
v(j)=incr(j)
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)
- KEt_sc=KEt_sc+msc(iti,mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(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
incr(j)=incr(j)+d_t(j,i)
real(kind=8) :: tt0,scalfac
integer :: nres2
nres2=2*nres
+ print *, "ENTER MD"
!
#ifdef MPI
+ print *,"MY tmpdir",tmpdir,ilen(tmpdir)
if (ilen(tmpdir).gt.0) &
call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_" &
//liczba(:ilen(liczba))//'.rst')
#else
tt0 = tcpu()
#endif
+ print *,"just befor setup matix",nres
! Determine the inverse of the inertia matrix.
call setup_MD_matrices
! Initialize MD
+ print *,"AFTER SETUP MATRICES"
call init_MD
+ print *,"AFTER INIT MD"
+
#ifdef MPI
t_MDsetup = MPI_Wtime()-tt0
#else
stop
#endif
else if (lang.eq.1 .or. lang.eq.4) then
+ print *,"before setup_fricmat"
call setup_fricmat
+ print *,"after setup_fricmat"
endif
#ifdef MPI
t_langsetup=MPI_Wtime()-tt0
#endif
endif
if (ntwe.ne.0) then
- if (mod(itime,ntwe).eq.0) call statout(itime)
+ if (mod(itime,ntwe).eq.0) then
+ call statout(itime)
+ call returnbox
+! call check_ecartint
+ endif
#ifdef VOUT
do j=1,3
v_work(j)=d_t(j,0)
enddo
enddo
do i=nnt,nct
- if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
+ mnum=molnum(i)
+ if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1.and.mnum.ne.5) then
do j=1,3
ind=ind+1
v_work(ind)=d_t(j,i+nres)
#endif
endif
if (mod(itime,ntwx).eq.0) then
+ call returnbox
write (tytul,'("time",f8.2)') totT
if(mdpdb) then
call hairpin(.true.,nharp,iharp)
enddo
endif
! Side chains
- if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
+ if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1.and.&
+ molnum(i).ne.5) then
do j=1,3
epdriftij= &
dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
write (iout,*) "vcm right after adjustment:"
write (iout,*) (vcm(j),j=1,3)
endif
- if ((.not.rest).and.(indpdb.eq.0)) then
- call chainbuild
- if(iranconf.ne.0) then
+ if (.not.rest) then
+! call chainbuild
+ if(iranconf.ne.0 .or.indpdb.gt.0.and..not.unres_pdb .or.preminim) then
if (overlapsc) then
print *, 'Calling OVERLAP_SC'
call overlap_sc(fail)
! include 'COMMON.NAMES'
! include 'COMMON.TIME1'
real(kind=8) :: xv,sigv,lowb,highb ,Ek1
+!#define DEBUG
#ifdef FIVEDIAG
real(kind=8) ,allocatable, dimension(:) :: DDU1,DDU2,DL2,DL1,xsolv,DML,rs
real(kind=8) :: sumx
#ifdef DEBUG
real(kind=8) ,allocatable, dimension(:) :: rsold
real (kind=8),allocatable,dimension(:,:) :: matold
+ integer :: iti
#endif
#endif
integer :: i,j,ii,k,ind,mark,imark,mnum
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).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
- .and.(mnum.ne.5)) 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
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*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
write (iout,*) "i",i," ii",ii
d_t(k,j)=d_t_work(ind)
ind=ind+1
enddo
- mnum=molnum(i)
- if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
- .and.(mnum.ne.5))
+ 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)
ind=ind+1
do i=nnt,nct
! if (itype(i,1).eq.10) then
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).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
+ .or.(mnum.eq.5)) then
do j=1,3
d_t(j,i)=d_t(j,i+1)-d_t(j,i)
enddo
! write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",&
! 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
! call flush(iout)
+! write(iout,*) "end init MD"
return
-#undef DEBUG
end subroutine random_vel
!-----------------------------------------------------------------------------
#ifndef LANG0
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
M_PEP=M_PEP+mp(mnum)
do j=1,3
M_SC=0.0d0
do i=nnt,nct
mnum=molnum(i)
+ if (mnum.eq.5) cycle
iti=iabs(itype(i,mnum))
M_SC=M_SC+msc(iabs(iti),mnum)
inres=i+nres
do i=nnt,nct-1
mnum=molnum(i)
+ if (mnum.eq.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.eq.5) cycle
inres=i+nres
do j=1,3
pr(j)=c(j,inres)-cm(j)
call angmom(cm,L)
! write(iout,*) "The angular momentum before adjustment:"
-! write(iout,*) (L(j),j=1,3)
+! write(iout,*) (L(j),j=1,3)
Im(2,1)=Im(1,2)
Im(3,1)=Im(1,3)
do i=1,3
do j=1,3
Imcp(i,j)=Im(i,j)
- Id(i,j)=0.0d0
+ Id(i,j)=0.0d0
enddo
enddo
enddo
call angmom(cm,L)
! write(iout,*) "The angular momentum after adjustment:"
-! write(iout,*) (L(j),j=1,3)
+! write(iout,*) (L(j),j=1,3)
return
end subroutine inertia_tensor
! include 'COMMON.INTERACT'
! include 'COMMON.IOUNITS'
! include 'COMMON.NAMES'
-
+ real(kind=8) :: mscab
real(kind=8),dimension(3) :: L,cm,pr,vp,vrot,incr,v,pp
integer :: iti,inres,i,j,mnum
! Calculate the angular momentum
enddo
do i=nnt,nct-1
mnum=molnum(i)
+ if (mnum.eq.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
mnum=molnum(i)
iti=iabs(itype(i,mnum))
inres=i+nres
+ if (mnum.eq.5) then
+ mscab=0.0d0
+ else
+ mscab=msc(iti,mnum)
+ endif
do j=1,3
pr(j)=c(j,inres)-cm(j)
enddo
enddo
endif
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)
+! write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),&
+! " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
do j=1,3
- L(j)=L(j)+msc(iabs(iti),mnum)*vp(j)
+ L(j)=L(j)+mscab*vp(j)
enddo
! write (iout,*) "L",(l(j),j=1,3)
if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
summas=0.0d0
do i=nnt,nct
mnum=molnum(i)
+ if (mnum.eq.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))
enddo
endif
+ if (mnum.ne.5) 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
logical :: lprn = .false.
real(kind=8) :: dtdi !el ,gamvec(2*nres)
!el real(kind=8),dimension(2*nres,2*nres) :: ginvfric,fcopy
- real(kind=8),dimension(2*nres,2*nres) :: fcopy
+! real(kind=8),allocatable,dimension(:,:) :: fcopy
!el real(kind=8),dimension(2*nres*(2*nres+1)/2) :: Ghalf !(mmaxres2) (mmaxres2=(maxres2*(maxres2+1)/2))
!el common /syfek/ gamvec
real(kind=8) :: work(8*2*nres)
integer :: iwork(2*nres)
!el common /przechowalnia/ ginvfric,Ghalf,fcopy
integer :: ii,iti,k,l,nzero,nres2,nres6,ierr,mnum
+ nres2=2*nres
+ nres6=6*nres
#ifdef MPI
+ 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
- nres2=2*nres
- nres6=6*nres
+! nres2=2*nres
+! nres6=6*nres
if(.not.allocated(gamvec)) allocate(gamvec(nres2)) !(MAXRES2)
if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres
-!el if(.not.allocated(fcopy)) allocate(fcopy(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
-!el if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
+ if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
! Zeroing out fricmat
do i=1,dimen
do j=1,dimen
! Load the friction coefficients corresponding to side chains
m=nct-nnt
ind=0
- do j=1,5
+ do j=1,2
gamsc(ntyp1_molec(j),j)=1.0d0
enddo
do i=nnt,nct
allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
#endif
-!el if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2))
+ if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2))
!----------------------
! commom.hairpin in CSA module
!----------------------
projects / unres4.git / blobdiff