!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
integer :: i,j,k,iti
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)
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
+! 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
! 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
! 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
+#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
+#endif
+#endif
+ integer :: i,j,ii,k,ind,mark,imark
! 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).eq.10) 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*((d_t_work(ii+j+k)+d_t_work_new(ii+k))/2)**2+&
+ 0.5d0*0.25d0*IP*(d_t_work(ii+j+k)-d_t_work_new(ii+k))**2
+ enddo
+ endif
+ if (itype(i).ne.10) ii=ii+3
+ write (iout,*) "i",i," itype",itype(i)," mass",msc(itype(i))
+ write (iout,*) "ii",ii
+ do k=1,3
+ ii=ii+1
+ write (iout,*) "k",k," ii",ii,"EK1",EK1
+ if (iabs(itype(i)).ne.10) Ek1=Ek1+0.5d0*Isc(iabs(itype(i)))*(d_t_work(ii)-d_t_work(ii-3))**2
+ Ek1=Ek1+0.5d0*msc(iabs(itype(i)))*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
+ if (itype(j).ne.10 .and. itype(j).ne.ntyp1) 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).eq.10) 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
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)
return
+#undef DEBUG
end subroutine random_vel
!-----------------------------------------------------------------------------
#ifndef LANG0
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
!
! 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
projects / unres4.git / blobdiff