d_t_work_new,d_af_work,d_as_work,kinetic_force !(MAXRES6)
real(kind=8),dimension(:,:),allocatable :: d_t_new,&
d_a_old,d_a_short!,d_a !(3,0:MAXRES2)
+! real(kind=8),dimension(:),allocatable :: d_a_work !(6*MAXRES)
+! real(kind=8),dimension(:,:),allocatable :: Gmat,Ginv,A,&
+! Gsqrp,Gsqrm,Gvec !(maxres2,maxres2)
+! real(kind=8),dimension(:),allocatable :: Geigen !(maxres2)
+! integer :: dimen,dimen1,dimen3
+! integer :: lang,count_reset_moment,count_reset_vel
+! logical :: reset_moment,reset_vel,rattle,RESPA
! common /inertia/
! common /langevin/
+! real(kind=8) :: rwat,etawat,stdfp,cPoise
+! real(kind=8),dimension(:),allocatable :: gamsc !(ntyp1)
+! real(kind=8),dimension(:),allocatable :: stdfsc !(ntyp)
real(kind=8),dimension(:),allocatable :: stdforcp,stdforcsc !(MAXRES)
!-----------------------------------------------------------------------------
! 'sizes.i'
!el integer maxamino,maxnuc,maxbnd
!el integer maxang,maxtors,maxpi
!el integer maxpib,maxpit
- integer :: maxatm !=2*nres !maxres2 maxres2=2*maxres
- integer,parameter :: maxval=8
- integer,parameter :: maxgrp=1000
- integer,parameter :: maxtyp=3000
- integer,parameter :: maxclass=500
- integer,parameter :: maxkey=10000
- integer,parameter :: maxrot=1000
- integer,parameter :: maxopt=1000
- integer,parameter :: maxhess=1000000
- integer :: maxlight !=8*maxatm
- integer,parameter :: maxvib=1000
- integer,parameter :: maxgeo=1000
- integer,parameter :: maxcell=10000
- integer,parameter :: maxring=10000
- integer,parameter :: maxfix=10000
- integer,parameter :: maxbio=10000
- integer,parameter :: maxamino=31
- integer,parameter :: maxnuc=12
- integer :: maxbnd !=2*maxatm
- integer :: maxang !=3*maxatm
- integer :: maxtors !=4*maxatm
- integer,parameter :: maxpi=100
- integer,parameter :: maxpib=2*maxpi
- integer,parameter :: maxpit=4*maxpi
+! integer :: maxatm !=2*nres !maxres2 maxres2=2*maxres
+! integer,parameter :: maxval=8
+! integer,parameter :: maxgrp=1000
+! integer,parameter :: maxtyp=3000
+! integer,parameter :: maxclass=500
+! integer,parameter :: maxkey=10000
+! integer,parameter :: maxrot=1000
+! integer,parameter :: maxopt=1000
+! integer,parameter :: maxhess=1000000
+! integer :: maxlight !=8*maxatm
+! integer,parameter :: maxvib=1000
+! integer,parameter :: maxgeo=1000
+! integer,parameter :: maxcell=10000
+! integer,parameter :: maxring=10000
+! integer,parameter :: maxfix=10000
+! integer,parameter :: maxbio=10000
+! integer,parameter :: maxamino=31
+! integer,parameter :: maxnuc=12
+! integer :: maxbnd !=2*maxatm
+! integer :: maxang !=3*maxatm
+! integer :: maxtors !=4*maxatm
+! integer,parameter :: maxpi=100
+! integer,parameter :: maxpib=2*maxpi
+! integer,parameter :: maxpit=4*maxpi
!-----------------------------------------------------------------------------
! Maximum number of seed
- integer,parameter :: max_seed=1
+! integer,parameter :: max_seed=1
!-----------------------------------------------------------------------------
real(kind=8),dimension(:),allocatable :: stochforcvec !(MAXRES6) maxres6=6*maxres
! common /stochcalc/ stochforcvec
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)
!-----------------------------------------------------------------------------
!
!
use control, only: tcpu
use control_data
use energy_data
+! use io_conf, only:cartprint
! include 'DIMENSIONS'
#ifdef MPI
include 'mpif.h'
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
nbond=nct-nnt
do i=nnt,nct
- if (itype(i).ne.10) nbond=nbond+1
+ if (itype(i,1).ne.10) nbond=nbond+1
enddo
!
if (lprn1) then
ind=ind+3
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
ind1=ind1+1
do j=1,3
Bmat(ind+j,ind1)=dC_norm(j,i+nres)
Td(i)=Td(i)+vbl*Tmat(i,ind)
enddo
do k=nnt,nct
- if (itype(k).ne.10) then
+ if (itype(k,1).ne.10) then
ind=ind+1
- Td(i)=Td(i)+vbldsc0(1,itype(k))*Tmat(i,ind)
+ Td(i)=Td(i)+vbldsc0(1,itype(k,1))*Tmat(i,ind)
endif
enddo
enddo
enddo
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
do j=1,3
ind=ind+1
zapas(ind)=-gxcart(j,i)+stochforcvec(ind)
i,(dC(j,i),j=1,3),xx
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
ind=ind+1
- xx=vbld(i+nres)-vbldsc0(1,itype(i))
+ xx=vbld(i+nres)-vbldsc0(1,itype(i,1))
write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
i,(dC(j,i+nres),j=1,3),xx
endif
endif
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
ind=ind+1
blen2 = scalar(dc(1,i+nres),dc(1,i+nres))
- ppvec(ind)=2*vbldsc0(1,itype(i))**2-blen2
- diffbond=dabs(vbldsc0(1,itype(i))-dsqrt(blen2))
+ ppvec(ind)=2*vbldsc0(1,itype(i,1))**2-blen2
+ diffbond=dabs(vbldsc0(1,itype(i,1))-dsqrt(blen2))
if (diffbond.gt.diffmax) diffmax=diffbond
if (ppvec(ind).gt.0.0d0) then
ppvec(ind)=dsqrt(ppvec(ind))
ind=ind+3
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
do j=1,3
dc(j,i+nres)=zapas(ind+j)
dc_work(ind+j)=zapas(ind+j)
i,(dC(j,i),j=1,3),xx
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
ind=ind+1
- xx=vbld(i+nres)-vbldsc0(1,itype(i))
+ xx=vbld(i+nres)-vbldsc0(1,itype(i,1))
write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
i,(dC(j,i+nres),j=1,3),xx
endif
potE=potEcomp(0)-potEcomp(20)
call cartgrad
totT=totT+d_time
+ totTafm=totT
! Calculate the kinetic and total energy and the kinetic temperature
call kinetic(EK)
#ifdef MPI
! 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
+ integer :: i,j,k,iti,mnum
real(kind=8) :: KEt_p,KEt_sc,KEr_p,KEr_sc,incr(3),&
mag1,mag2,v(3)
-
+#ifdef DEBUG
+ write (iout,*) "Velocities, kietic"
+ do i=0,nres
+ write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
+ (d_t(j,i+nres),j=1,3)
+ enddo
+#endif
KEt_p=0.0d0
KEt_sc=0.0d0
-! write (iout,*) "ISC",(isc(itype(i)),i=1,nres)
+! write (iout,*) "ISC",(isc(itype(i,1)),i=1,nres)
! The translational part for peptide virtual bonds
do j=1,3
incr(j)=d_t(j,0)
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)
enddo
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))
+ KEt_p=KEt_p+mp(mnum)*(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
incr(j)=d_t(j,0)
enddo
do i=nnt,nct
- iti=iabs(itype(i))
- if (itype(i).eq.10) then
+ 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).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)*(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)
! The part due to stretching and rotation of the peptide groups
KEr_p=0.0D0
do i=nnt,nct-1
+ mnum=molnum(i)
! write (iout,*) "i",i
! write (iout,*) "i",i," mag1",mag1," mag2",mag2
do j=1,3
incr(j)=d_t(j,i)
enddo
! write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3)
- KEr_p=KEr_p+(incr(1)*incr(1)+incr(2)*incr(2) &
+ KEr_p=KEr_p+Ip(mnum)*(incr(1)*incr(1)+incr(2)*incr(2) &
+incr(3)*incr(3))
enddo
! goto 111
! The rotational part of the side chain virtual bond
KEr_sc=0.0D0
do i=nnt,nct
- iti=iabs(itype(i))
- if (itype(i).ne.10) then
+ mnum=molnum(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
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)*(incr(1)*incr(1)+incr(2)*incr(2)+ &
+ KEr_sc=KEr_sc+Isc(iti,mnum)*(incr(1)*incr(1)+incr(2)*incr(2)+ &
incr(3)*incr(3))
endif
enddo
! The total kinetic energy
111 continue
! write(iout,*) 'KEr_sc', KEr_sc
- KE_total=0.5d0*(mp*KEt_p+KEt_sc+0.25d0*Ip*KEr_p+KEr_sc)
+ KE_total=0.5d0*(KEt_p+KEt_sc+0.25d0*KEr_p+KEr_sc)
! write (iout,*) "KE_total",KE_total
return
end subroutine kinetic
! The driver for molecular dynamics subroutines
!------------------------------------------------
use comm_gucio
+! use MPI
use control, only:tcpu,ovrtim
+! use io_comm, only:ilen
use control_data
use compare, only:secondary2,hairpin
use io, only:cartout,statout
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).ne.10 .and. itype(i).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)
if (rstcount.eq.1000.or.itime.eq.n_timestep) then
open(irest2,file=rest2name,status='unknown')
write(irest2,*) totT,EK,potE,totE,t_bath
- do i=1,2*nres
+ totTafm=totT
+! AL 4/17/17: Now writing d_t(0,:) too
+ do i=0,2*nres
write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
enddo
- do i=1,2*nres
+! AL 4/17/17: Now writing d_c(0,:) too
+ do i=0,2*nres
write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
enddo
close(irest2)
! Calculate energy and forces
call zerograd
call etotal(potEcomp)
+! AL 4/17/17: Reduce the steps if NaNs occurred.
+ if (potEcomp(0).gt.0.99e20 .or. isnan(potEcomp(0)).gt.0) then
+ d_time=d_time/2
+ cycle
+ endif
+! end change
if (large.and. mod(itime,ntwe).eq.0) &
call enerprint(potEcomp)
#ifdef TIMING_ENE
endif
if (rattle) call rattle2
totT=totT+d_time
+ totTafm=totT
if (d_time.ne.d_time0) then
d_time=d_time0
#ifndef LANG0
! include 'DIMENSIONS'
use comm_gucio
use comm_cipiszcze
+! use MPI
use control, only:tcpu
use control_data
+! use io_conf, only:cartprint
#ifdef MPI
include 'mpif.h'
integer :: IERROR,ERRCODE
! Calculate energy and forces
call zerograd
call etotal_short(energia_short)
+! AL 4/17/17: Exit itime_split loop when energy goes infinite
+ if (energia_short(0).gt.0.99e20 .or. isnan(energia_short(0)) ) then
+ if (PRINT_AMTS_MSG) &
+ write (iout,*) "Infinities/NaNs in energia_short",energia_short(0),"; increasing ntime_split to",ntime_split
+ ntime_split=ntime_split*2
+ if (ntime_split.gt.maxtime_split) then
+#ifdef MPI
+ write (iout,*) &
+ "Cannot rescue the run; aborting job. Retry with a smaller time step"
+ call flush(iout)
+ call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
+#else
+ write (iout,*) &
+ "Cannot rescue the run; terminating. Retry with a smaller time step"
+#endif
+ endif
+ exit
+ endif
+! End change
if (large.and. mod(itime,ntwe).eq.0) &
call enerprint(energia_short)
#ifdef TIMING_ENE
if (ntime_split.lt.maxtime_split) then
scale=.true.
ntime_split=ntime_split*2
+! AL 4/17/17: We should exit the itime_split loop when acceleration change is too big
+ exit
do i=0,2*nres
do j=1,3
dc_old(j,i)=dc_old0(j,i)
#endif
call zerograd
call etotal_long(energia_long)
+ if (energia_long(0).gt.0.99e20 .or. isnan(energia_long(0))) then
+#ifdef MPI
+ write (iout,*) &
+ "Infinitied/NaNs in energia_long, Aborting MPI job."
+ call flush(iout)
+ call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
+#else
+ write (iout,*) "Infinitied/NaNs in energia_long, terminating."
+ stop
+#endif
+ endif
if (large.and. mod(itime,ntwe).eq.0) &
call enerprint(energia_long)
#ifdef TIMING_ENE
potE=potEcomp(0)-potEcomp(20)
! potE=energia_short(0)+energia_long(0)
totT=totT+d_time
+ totTafm=totT
! Calculate the kinetic and the total energy and the kinetic temperature
call kinetic(EK)
totE=EK+potE
! include 'COMMON.INTERACT'
! include 'COMMON.IOUNITS'
! include 'COMMON.NAMES'
- integer :: i,j,inres
+ integer :: i,j,inres,mnum
do j=1,3
d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
enddo
enddo
do i=nnt,nct
- if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ 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
inres=i+nres
do j=1,3
d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
! include 'COMMON.IOUNITS'
! include 'COMMON.NAMES'
real(kind=8) :: adt,adt2
- integer :: i,j,inres
+ integer :: i,j,inres,mnum
#ifdef DEBUG
write (iout,*) "VELVERLET1 START: DC"
enddo
enddo
do i=nnt,nct
- if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ 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
inres=i+nres
do j=1,3
adt=d_a_old(j,inres)*d_time
! include 'COMMON.INTERACT'
! include 'COMMON.IOUNITS'
! include 'COMMON.NAMES'
- integer :: i,j,inres
+ integer :: i,j,inres,mnum
do j=1,3
d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
enddo
enddo
do i=nnt,nct
- if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ mnum=molnum(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
inres=i+nres
do j=1,3
d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
! position and velocity increments included.
real(kind=8) :: sqrt13 = 0.57735026918962576451d0 ! 1/sqrt(3)
real(kind=8) :: adt,adt2
- integer :: i,j,ind,inres
+ integer :: i,j,ind,inres,mnum
!
! Add the contribution from BOTH friction and stochastic force to the
! coordinates, but ONLY the contribution from the friction forces to velocities
ind=ind+3
enddo
do i=nnt,nct
- if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ mnum=molnum(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
inres=i+nres
do j=1,3
adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
! include 'COMMON.NAMES'
real(kind=8),dimension(6*nres) :: stochforcvec,d_as_work1 !(MAXRES6) maxres6=6*maxres
real(kind=8) :: cos60 = 0.5d0, sin60 = 0.86602540378443864676d0
- integer :: i,j,ind,inres
+ integer :: i,j,ind,inres,mnum
! Revised 3/31/05 AL: correlation between random contributions to
! position and velocity increments included.
! The correlation coefficients are calculated at low-friction limit.
ind=ind+3
enddo
do i=nnt,nct
- if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ mnum=molnum(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
inres=i+nres
do j=1,3
d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres) &
! include 'COMMON.IOUNITS'
real(kind=8),dimension(3) :: aux,accel,accel_old
real(kind=8) :: dacc
- integer :: i,j
+ integer :: i,j,mnum
do j=1,3
! aux(j)=d_a(j,0)-d_a_old(j,0)
enddo
endif
do i=nnt,nct
- if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ mnum=molnum(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
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)
enddo
endif
! Side chains
- if (itype(i).ne.10 .and. itype(i).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))
! include 'COMMON.IOUNITS'
! include 'COMMON.NAMES'
real(kind=8) :: T_half,fact
- integer :: i,j,inres
+ integer :: i,j,inres,mnum
!
T_half=2.0d0/(dimen3*Rb)*EK
fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
enddo
enddo
do i=nnt,nct
- if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ mnum=molnum(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
inres=i+nres
do j=1,3
d_t(j,inres)=fact*d_t(j,inres)
character(len=50) :: tytul
logical :: file_exist
!el common /gucio/ cm
- integer :: i,j,ipos,iq,iw,nft_sc,iretcode,nfun,itime,ierr
+ integer :: i,j,ipos,iq,iw,nft_sc,iretcode,nfun,itime,ierr,mnum
real(kind=8) :: etot,tt0
logical :: fail
! if the friction coefficients do not depend on surface area
if (lang.gt.0 .and. .not.surfarea) then
do i=nnt,nct-1
- stdforcp(i)=stdfp*dsqrt(gamp)
+ mnum=molnum(i)
+ stdforcp(i)=stdfp(mnum)*dsqrt(gamp(mnum))
enddo
do i=nnt,nct
- stdforcsc(i)=stdfsc(iabs(itype(i))) &
- *dsqrt(gamsc(iabs(itype(i))))
+ mnum=molnum(i)
+ stdforcsc(i)=stdfsc(iabs(itype(i,mnum)),mnum) &
+ *dsqrt(gamsc(iabs(itype(i,mnum)),mnum))
enddo
endif
! Open the pdb file for snapshotshots
endif
call random_vel
totT=0.0d0
+ totTafm=totT
endif
else
! Generate initial velocities
write(iout,*) "Initial velocities randomly generated"
call random_vel
totT=0.0d0
+ totTafm=totT
endif
! rest2name = prefix(:ilen(prefix))//'.rst'
if(me.eq.king.or..not.out1file)then
write (iout,*) (vcm(j),j=1,3)
endif
if (.not.rest) then
- call chainbuild
- if(iranconf.ne.0) 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)
endif
endif
call chainbuild_cart
-write(iout,*) "przed kinetic, EK=",EK
call kinetic(EK)
if (tbf) then
call verlet_bath
endif
-write(iout,*) "dimen3",dimen3,"Rb",Rb,"EK",EK
kinetic_T=2.0d0/(dimen3*Rb)*EK
-write(iout,*) "kinetic_T",kinetic_T
if(me.eq.king.or..not.out1file)then
call cartprint
call intout
#endif
potE=potEcomp(0)
call cartgrad
-write(iout,*) "kinetic_T if large",kinetic_T
call lagrangian
-write(iout,*) "kinetic_T if large",kinetic_T
call max_accel
if (amax*d_time .gt. dvmax) then
d_time=d_time*dvmax/amax
#endif
#endif
call cartgrad
-write(iout,*) "przed lagrangian"
call lagrangian
-write(iout,*) "po lagrangian"
if(.not.out1file .and. large) then
write (iout,*) "energia_long",energia_long(0),&
" energia_short",energia_short(0),&
#endif
#endif
call cartgrad
-write(iout,*) "przed lagrangian2"
call lagrangian
-write(iout,*) "po lagrangian2"
if(.not.out1file .and. large) then
write (iout,*) "energia_long",energia_long(0)
write (iout,*) "Initial slow-force accelerations"
t_enegrad=t_enegrad+tcpu()-tt0
#endif
endif
-write(iout,*) "end init MD"
return
end subroutine init_MD
!-----------------------------------------------------------------------------
! implicit real*8 (a-h,o-z)
use energy_data
use random, only:anorm_distr
+ use MD_data
! include 'DIMENSIONS'
! include 'COMMON.CONTROL'
! include 'COMMON.VAR'
! include 'COMMON.NAMES'
! include 'COMMON.TIME1'
real(kind=8) :: xv,sigv,lowb,highb ,Ek1
- integer :: i,j,ii,k,ind
+!#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
! Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
! First generate velocities in the eigenspace of the G matrix
! write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
lowb=-5*sigv
highb=5*sigv
d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
-! write (iout,*) "i",i," ii",ii," geigen",geigen(i),&
-! " d_t_work_new",d_t_work_new(ii)
+#ifdef DEBUG
+ write (iout,*) "i",i," ii",ii," geigen",geigen(i),&
+ " d_t_work_new",d_t_work_new(ii)
+#endif
enddo
enddo
+#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
+ 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
+#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))
+#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)
+#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
+#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
+#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
+#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
+ endif
+ if (i.lt.nct) then
+ 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
+ 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
+ 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
+ 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
+ deallocate(DDU1)
+ deallocate(DDU2)
+ deallocate(DL2)
+ deallocate(DL1)
+ deallocate(xsolv)
+ deallocate(DML)
+ deallocate(rs)
+#ifdef DEBUG
+ deallocate(matold)
+ deallocate(rsold)
+#endif
+ ind=1
+ do j=nnt,nct
+ do k=1,3
+ d_t(k,j)=d_t_work(ind)
+ ind=ind+1
+ 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)
+ ind=ind+1
+ enddo
+ end if
+ 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)
+ enddo
+ do i=nnt,nct
+ write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3)
+ enddo
+#endif
+ do j=1,3
+ d_t(j,0)=d_t(j,nnt)
+ 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
+ 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
+#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)
+ enddo
+ do i=nnt,nct
+ write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3)
+ enddo
+ call kinetic(Ek1)
+ write (iout,*) "Ek from d_t_work",Ek1
+ write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
+#endif
+#else
do k=0,2
do i=1,dimen
ind=(i-1)*3+k+1
enddo
enddo
do i=nnt,nct
- if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ 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
do j=1,3
ind=ind+1
d_t(j,i+nres)=d_t_work(ind)
enddo
endif
enddo
+#endif
! call kinetic(EK)
! 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
end subroutine random_vel
!-----------------------------------------------------------------------------
ind=ind+3
enddo
do i=nnt,nct
- if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ 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
do j=1,3
dc_work(ind+j)=dc_old(j,i+nres)
d_t_work(ind+j)=d_t_old(j,i+nres)
ind=ind+3
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ 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
inres=i+nres
do j=1,3
dc(j,inres)=dc_work(ind+j)
ind=ind+3
enddo
do i=nnt,nct
- if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ 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
inres=i+nres
do j=1,3
d_t(j,inres)=d_t_work(ind+j)
ind=ind+3
enddo
do i=nnt,nct
- if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) 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)
ind=ind+3
enddo
do i=nnt,nct
- if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ 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
inres=i+nres
do j=1,3
dc(j,inres)=dc_work(ind+j)
ind=ind+3
enddo
do i=nnt,nct
- if (itype(i).ne.10 .and. itype(i).ne.ntyp1) 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).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)
real(kind=8),dimension(3,3) :: Im,Imcp,eigvec,Id
real(kind=8),dimension(3) :: pr,eigval,L,vp,vrot
- real(kind=8) :: M_SC,mag,mag2
+ real(kind=8) :: M_SC,mag,mag2,M_PEP
real(kind=8),dimension(3,0:nres) :: vpp !(3,0:MAXRES)
real(kind=8),dimension(3) :: vs_p,pp,incr,v
real(kind=8),dimension(3,3) :: pr1,pr2
!el common /gucio/ cm
- integer :: iti,inres,i,j,k
+ integer :: iti,inres,i,j,k,mnum
do i=1,3
do j=1,3
Im(i,j)=0.0d0
vrot(i)=0.0d0
enddo
! calculating the center of the mass of the protein
+ 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
- cm(j)=cm(j)+c(j,i)+0.5d0*dc(j,i)
+ 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
+! do j=1,3
+! cm(j)=mp(1)*cm(j)
+! enddo
M_SC=0.0d0
do i=nnt,nct
- iti=iabs(itype(i))
- M_SC=M_SC+msc(iabs(iti))
+ 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 j=1,3
- cm(j)=cm(j)+msc(iabs(iti))*c(j,inres)
+ cm(j)=cm(j)+msc(iabs(iti),mnum)*c(j,inres)
enddo
enddo
do j=1,3
- cm(j)=cm(j)/(M_SC+(nct-nnt)*mp)
+ cm(j)=cm(j)/(M_SC+M_PEP)
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
- Im(1,1)=Im(1,1)+mp*(pr(2)*pr(2)+pr(3)*pr(3))
- Im(1,2)=Im(1,2)-mp*pr(1)*pr(2)
- Im(1,3)=Im(1,3)-mp*pr(1)*pr(3)
- Im(2,3)=Im(2,3)-mp*pr(2)*pr(3)
- Im(2,2)=Im(2,2)+mp*(pr(3)*pr(3)+pr(1)*pr(1))
- Im(3,3)=Im(3,3)+mp*(pr(1)*pr(1)+pr(2)*pr(2))
+ 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
- iti=iabs(itype(i))
+ 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)
enddo
- Im(1,1)=Im(1,1)+msc(iabs(iti))*(pr(2)*pr(2)+pr(3)*pr(3))
- Im(1,2)=Im(1,2)-msc(iabs(iti))*pr(1)*pr(2)
- Im(1,3)=Im(1,3)-msc(iabs(iti))*pr(1)*pr(3)
- Im(2,3)=Im(2,3)-msc(iabs(iti))*pr(2)*pr(3)
- Im(2,2)=Im(2,2)+msc(iabs(iti))*(pr(3)*pr(3)+pr(1)*pr(1))
- Im(3,3)=Im(3,3)+msc(iabs(iti))*(pr(1)*pr(1)+pr(2)*pr(2))
+ 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
- Im(1,1)=Im(1,1)+Ip*(1-dc_norm(1,i)*dc_norm(1,i))* &
+ mnum=molnum(i)
+ 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*(-dc_norm(1,i)*dc_norm(2,i))* &
+ 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*(-dc_norm(1,i)*dc_norm(3,i))* &
+ 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*(-dc_norm(2,i)*dc_norm(3,i))* &
+ 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*(1-dc_norm(2,i)*dc_norm(2,i))* &
+ 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*(1-dc_norm(3,i)*dc_norm(3,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
do i=nnt,nct
- if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
- iti=iabs(itype(i))
+ 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
+ iti=iabs(itype(i,mnum))
inres=i+nres
- Im(1,1)=Im(1,1)+Isc(iti)*(1-dc_norm(1,inres)* &
+ 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)*(dc_norm(1,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)*(dc_norm(1,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)*(dc_norm(2,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)*(1-dc_norm(2,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)*(1-dc_norm(3,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)
! 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
enddo
do i=nnt,nct
- if(itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ mnum=molnum(i)
+ if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
+ .and.(mnum.ne.5)) 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
call vecpr(vrot(1),dc(1,inres),vp)
do j=1,3
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
+ integer :: iti,inres,i,j,mnum
! Calculate the angular momentum
do j=1,3
L(j)=0.0d0
incr(j)=d_t(j,0)
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
enddo
call vecpr(pr(1),v(1),vp)
do j=1,3
- L(j)=L(j)+mp*vp(j)
+ L(j)=L(j)+mp(mnum)*vp(j)
enddo
do j=1,3
pr(j)=0.5d0*dc(j,i)
enddo
call vecpr(pr(1),pp(1),vp)
do j=1,3
- L(j)=L(j)+Ip*vp(j)
+ L(j)=L(j)+Ip(mnum)*vp(j)
enddo
enddo
do j=1,3
incr(j)=d_t(j,0)
enddo
do i=nnt,nct
- iti=iabs(itype(i))
+ 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
- if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
+ .and.(mnum.ne.5)) then
do j=1,3
v(j)=incr(j)+d_t(j,inres)
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))*vp(j)
+ L(j)=L(j)+mscab*vp(j)
enddo
! write (iout,*) "L",(l(j),j=1,3)
- 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.ne.5)) 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)*vp(j)
+ L(j)=L(j)+Isc(iti,mnum)*vp(j)
enddo
endif
do j=1,3
! include 'COMMON.IOUNITS'
real(kind=8),dimension(3) :: vcm,vv
real(kind=8) :: summas,amas
- integer :: i,j
+ integer :: i,j,mnum
do j=1,3
vcm(j)=0.0d0
enddo
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
+ summas=summas+mp(mnum)
do j=1,3
- vcm(j)=vcm(j)+mp*(vv(j)+0.5d0*d_t(j,i))
+ vcm(j)=vcm(j)+mp(mnum)*(vv(j)+0.5d0*d_t(j,i))
enddo
endif
- amas=msc(iabs(itype(i)))
+ if (mnum.ne.5) then
+ amas=msc(iabs(itype(i,mnum)),mnum)
+ else
+ amas=0.0d0
+ endif
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.ne.5)) then
do j=1,3
vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres))
enddo
if (lprn) write (iout,*) "RATTLE1"
nbond=nct-nnt
do i=nnt,nct
- if (itype(i).ne.10) nbond=nbond+1
+ mnum=molnum(i)
+ if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
+ .and.(mnum.ne.5)) nbond=nbond+1
enddo
! Make a folded form of the Ginv-matrix
ind=0
enddo
enddo
do k=nnt,nct
- if (itype(k).ne.10) then
+ mnum=molnum(k)
+ if (itype(k,1).ne.10 .and. itype(k,mnum).ne.ntyp1_molec(mnum)&
+ .and.(mnum.ne.5)) then
jj=jj+1
do l=1,3
ind1=ind1+1
enddo
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ mnum=molnum(i)
+ if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
+ .and.(mnum.ne.5))
ii=ii+1
do j=1,3
ind=ind+1
enddo
enddo
do k=nnt,nct
- if (itype(k).ne.10) then
+ if (itype(k,1).ne.10) then
jj=jj+1
do l=1,3
ind1=ind1+1
enddo
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
ind1=ind1+1
do j=1,3
dC_uncor(j,ind1)=dC(j,i+nres)
enddo
enddo
do k=nnt,nct
- if (itype(k).ne.10) then
+ if (itype(k,1).ne.10) then
ind=ind+1
do j=1,3
gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
x(ind)=vbld(i+1)**2-vbl**2
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
ind=ind+1
- x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i))**2
+ x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
endif
enddo
if (lprn) then
i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ mnum=molnum(i)
+ if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
+ .and.(mnum.ne.5)) then
+
ind=ind+1
write (iout,'(2i5,3f10.5,5x,e15.5)') &
i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
enddo
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
ind=ind+1
do j=1,3
xx=0.0d0
i,(dC(j,i),j=1,3),x(ind),xx
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ mnum=molnum(i)
+ if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
+ .and.(mnum.ne.5))
ind=ind+1
- xx=vbld(i+nres)**2-vbldsc0(1,itype(i))**2
+ xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
i,(dC(j,i+nres),j=1,3),x(ind),xx
endif
i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
ind=ind+1
write (iout,'(2i5,3f10.5,5x,e15.5)') &
i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
enddo
enddo
do k=nnt,nct
- if (itype(k).ne.10) then
+ mnum=molnum(i)
+ if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
+ .and.(mnum.ne.5)) then
ind=ind+1
do j=1,3
gdc(j,i,ind)=GGinv(i,ind)*dC(j,k+nres)
enddo
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ mnum=molnum(i)
+ if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
+ .and.(mnum.ne.5)) then
ind=ind+1
do j=1,nbond
Cmat(ind,j)=0.0d0
x(ind)=scalar(d_t(1,i),dC(1,i))
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ mnum=molnum(i)
+ if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
+ .and.(mnum.ne.5)) then
ind=ind+1
x(ind)=scalar(d_t(1,i+nres),dC(1,i+nres))
endif
i,ind,(d_t(j,i),j=1,3),x(ind)
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ mnum=molnum(i)
+ if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
+ .and.(mnum.ne.5)) 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)
enddo
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ mnum=molnum(i)
+ if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
+ .and.(mnum.ne.5)) then
ind=ind+1
do j=1,3
xx=0.0d0
i,ind,(d_t(j,i),j=1,3),x(ind),scalar(d_t(1,i),dC(1,i))
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ mnum=molnum(i)
+ if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
+ .and.(mnum.ne.5))
ind=ind+1
write (iout,'(2i5,3f10.5,5x,2e15.5)') &
i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind),&
!el common /przechowalnia/ GGinv,gdc,Cmat,nbond
!el common /przechowalnia/ nbond
integer :: max_rattle = 5
- logical :: lprn = .true., lprn1 = .true., not_done
+ logical :: lprn = .false., lprn1 = .false., not_done
real(kind=8) :: tol_rattle = 1.0d-5
integer :: nres2
nres2=2*nres
if (lprn) write (iout,*) "RATTLE_BROWN"
nbond=nct-nnt
do i=nnt,nct
- if (itype(i).ne.10) nbond=nbond+1
+ if (itype(i,1).ne.10) nbond=nbond+1
enddo
! Make a folded form of the Ginv-matrix
ind=0
enddo
enddo
do k=nnt,nct
- if (itype(k).ne.10) then
+ if (itype(k,1).ne.10) then
jj=jj+1
do l=1,3
ind1=ind1+1
enddo
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
ii=ii+1
do j=1,3
ind=ind+1
enddo
enddo
do k=nnt,nct
- if (itype(k).ne.10) then
+ if (itype(k,1).ne.10) then
jj=jj+1
do l=1,3
ind1=ind1+1
enddo
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
ind1=ind1+1
do j=1,3
dC_uncor(j,ind1)=dC(j,i+nres)
enddo
enddo
do k=nnt,nct
- if (itype(k).ne.10) then
+ if (itype(k,1).ne.10) then
ind=ind+1
do j=1,3
gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
x(ind)=vbld(i+1)**2-vbl**2
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
ind=ind+1
- x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i))**2
+ x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
endif
enddo
if (lprn) then
i,ind,(d_t(j,i),j=1,3),scalar(d_t(1,i),dC_old(1,i))
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
ind=ind+1
write (iout,'(2i5,3f10.5,5x,e15.5)') &
i+nres,ind,(d_t(j,i+nres),j=1,3),&
enddo
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
ind=ind+1
do j=1,3
xx=0.0d0
i,(dC(j,i),j=1,3),x(ind),xx
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
ind=ind+1
- xx=vbld(i+nres)**2-vbldsc0(1,itype(i))**2
+ xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2
write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') &
i,(dC(j,i+nres),j=1,3),x(ind),xx
endif
i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
enddo
do i=nnt,nct
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
ind=ind+1
write (iout,'(2i5,3f10.5,5x,e15.5)') &
i+nres,ind,(d_t_new(j,i+nres),j=1,3),&
!el common /przechowalnia/ ginvfric
logical :: lprn = .false., checkmode = .false.
- integer :: i,j,ind,k,nres2,nres6
+ integer :: i,j,ind,k,nres2,nres6,mnum
nres2=2*nres
nres6=6*nres
ind=ind+3
enddo
do i=nnt,nct
- if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
+ mnum=molnum(i)
+ if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
+ .and.(mnum.ne.5)) then
do j=1,3
d_t_work(ind+j)=d_t(j,i+nres)
enddo
ind=ind+3
enddo
do i=nnt,nct
- if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
+ mnum=molnum(i)
+ if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
+ .and.(mnum.ne.5)) 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
!-----------------------------------------------------------------------------
subroutine setup_fricmat
+! use MPI
use energy_data
use control_data, only:time_Bcast
use control, only:tcpu
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
+ 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 peptide groups
ind1=0
do i=nnt,nct-1
+ mnum=molnum(i)
ind1=ind1+1
- gamvec(ind1)=gamp
+ gamvec(ind1)=gamp(mnum)
enddo
! Load the friction coefficients corresponding to side chains
m=nct-nnt
ind=0
- gamsc(ntyp1)=1.0d0
+ 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)
- gamvec(ii)=gamsc(iabs(iti))
+ iti=itype(i,mnum)
+ gamvec(ii)=gamsc(iabs(iti),mnum)
enddo
if (surfarea) call sdarea(gamvec)
! if (lprn) then
#else
time00=tcpu()
#endif
-write(iout,*)"przed MPI_Scatterv in fricmat"
! Scatter the friction matrix
call MPI_Scatterv(fricmat(1,1),nginv_counts(0),&
nginv_start(0),MPI_DOUBLE_PRECISION,fcopy(1,1),&
myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR)
-write(iout,*)"po MPI_Scatterv in fricmat"
#ifdef TIMING
#ifdef MPI
time_scatter=time_scatter+MPI_Wtime()-time00
time_scatter_fmat=time_scatter_fmat+tcpu()-time00
#endif
#endif
-write(iout,*)"po MPI_Scatterv in fricmat"
do i=1,dimen
do j=1,2*my_ng_count
fricmat(j,i)=fcopy(i,j)
enddo
enddo
-write(iout,*)"po MPI_Scatterv in fricmat"
! write (iout,*) "My chunk of fricmat"
! call MATOUT2(my_ng_count,dimen,maxres2,maxres2,fcopy)
endif
real(kind=8) :: highb2,sig2,forcvec(6*nres),stochforcvec(6*nres)
real(kind=8) :: time00
logical :: lprn = .false.
- integer :: i,j,ind
+ integer :: i,j,ind,mnum
do i=0,2*nres
do j=1,3
do j=1,3
ff(j)=ff(j)+force(j,i)
enddo
- if (itype(i+1).ne.ntyp1) then
+! if (itype(i+1,1).ne.ntyp1) then
+ mnum=molnum(i)
+ if (itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
do j=1,3
stochforc(j,i)=stochforc(j,i)+force(j,i+nres+1)
ff(j)=ff(j)+force(j,i+nres+1)
stochforc(j,0)=ff(j)+force(j,nnt+nres)
enddo
do i=nnt,nct
- if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
+ mnum=molnum(i)
+ if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
+ .and.(mnum.ne.5)) 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)
enddo
ind=ind+3
enddo
do i=nnt,nct
- if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
+ mnum=molnum(i)
+ if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))&
+ .and.(mnum.ne.5)) 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
real(kind=8),parameter :: twosix = 1.122462048309372981d0
logical :: lprn = .false.
real(kind=8) :: probe,area,ratio
- integer :: i,j,ind,iti
+ integer :: i,j,ind,iti,mnum
!
! determine new friction coefficients every few SD steps
!
enddo
! Load peptide group radii
do i=nnt,nct-1
- radius(i)=pstok
+ mnum=molnum(i)
+ radius(i)=pstok(mnum)
enddo
! Load side chain radii
do i=nnt,nct
- iti=itype(i)
- radius(i+nres)=restok(iti)
+ mnum=molnum(i)
+ iti=itype(i,mnum)
+ radius(i+nres)=restok(iti,mnum)
enddo
! do i=1,2*nres
! write (iout,*) "i",i," radius",radius(i)
ind=ind+1
gamvec(ind) = ratio * gamvec(ind)
enddo
- stdforcp(i)=stdfp*dsqrt(gamvec(ind))
+ mnum=molnum(i)
+ stdforcp(i)=stdfp(mnum)*dsqrt(gamvec(ind))
if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcp(i)
endif
enddo
ind=ind+1
gamvec(ind) = ratio * gamvec(ind)
enddo
- stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamvec(ind))
+ mnum=molnum(i)
+ stdforcsc(i)=stdfsc(itype(i,mnum),mnum)*dsqrt(gamvec(ind))
if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcsc(i)
endif
enddo
logical :: omit(maxarc)
!
! include 'sizes.i'
- maxatm = 2*nres !maxres2 maxres2=2*maxres
- maxlight = 8*maxatm
- maxbnd = 2*maxatm
- maxang = 3*maxatm
- maxtors = 4*maxatm
+! maxatm = 2*nres !maxres2 maxres2=2*maxres
+! maxlight = 8*maxatm
+! maxbnd = 2*maxatm
+! maxang = 3*maxatm
+! maxtors = 4*maxatm
!
! zero out the surface area for the sphere of interest
!
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
!----------------------
! common /lagrange/
allocate(d_t(3,0:nres2),d_a(3,0:nres2),d_t_old(3,0:nres2)) !(3,0:MAXRES2)
allocate(d_a_work(nres6)) !(6*MAXRES)
+#ifdef FIVEDIAG
+ allocate(DM(nres2),DU1(nres2),DU2(nres2))
+ allocate(DMorig(nres2),DU1orig(nres2),DU2orig(nres2))
+#else
allocate(Gmat(nres2,nres2),A(nres2,nres2))
if(.not.allocated(Ginv)) allocate(Ginv(nres2,nres2)) !in control: ergastulum
allocate(Gsqrp(nres2,nres2),Gsqrm(nres2,nres2),Gvec(nres2,nres2)) !(maxres2,maxres2)
+#endif
allocate(Geigen(nres2)) !(maxres2)
if(.not.allocated(vtot)) allocate(vtot(nres2)) !(maxres2)
! common /inertia/ in io_conf: parmread
if(.not.allocated(nginv_start)) allocate(nginv_start(0:nfgtasks)) !(0:MaxProcs)
!----------------------
! common.muca in read_muca
+! common /double_muca/
+! real(kind=8) :: elow,ehigh,factor,hbin,factor_min
+! real(kind=8),dimension(:),allocatable :: emuca,nemuca,&
+! nemuca2,hist !(4*maxres)
+! common /integer_muca/
+! integer :: nmuca,imtime,muca_smooth
+! common /mucarem/
+! real(kind=8),dimension(:),allocatable :: elowi,ehighi !(maxprocs)
!----------------------
! common.MD
! common /mdgrad/ in module.energy
projects / unres4.git / blobdiff