#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)
+C call enerprint(potEcomp)
+C print *,itime,'AFM',Eafmforc,etot
+ 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.21) then
+ if (itype(i).ne.10 .and. itype(i).ne.ntyp1) 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
+c write(iout,*) 'time=',itime
+C call check_ecartint
+ call returnbox
write (tytul,'("time",f8.2)') totT
if(mdpdb) then
call hairpin(.true.,nharp,iharp)
c 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
+C print *,totTafm,"TU?"
if (d_time.ne.d_time0) then
d_time=d_time0
#ifndef LANG0
include 'COMMON.IOUNITS'
include 'COMMON.NAMES'
include 'COMMON.TIME1'
+ logical lprint_short
+ common /shortcheck/ lprint_short
double precision energia_short(0:n_ene),
& energia_long(0:n_ene)
double precision cm(3),L(3),vcm(3),incr(3)
write (iout,*) "***************** RESPA itime",itime
write (iout,*) "Cartesian and internal coordinates: step 0"
c call cartprint
- call pdbout(0.0d0,"cipiszcze",iout)
+ call pdbout(0.0d0,
+ & "cipiszcze ",iout)
call intout
write (iout,*) "Accelerations from long-range forces"
do i=0,nres
if (large.and. mod(itime,ntwe).eq.0) then
write (iout,*) "***** ITSPLIT",itsplit
write (iout,*) "Cartesian and internal coordinates: step 1"
- call pdbout(0.0d0,"cipiszcze",iout)
+ call pdbout(0.0d0,
+ & "cipiszcze ",iout)
c call cartprint
call intout
write (iout,*) "Velocities, step 1"
tt0 = tcpu()
#endif
c Calculate energy and forces
+c if (large.and. mod(itime,ntwe).eq.0) lprint_short=.true.
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
#endif
call zerograd
call etotal_long(energia_long)
+! AL 4/17/2017 Handling NaNs
+ 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 enerprint(energia_long)
+ call flush(iout)
+ call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
+#else
+ write (iout,*) "Infinitied/NaNs in energia_long, terminating."
+ call enerprint(energia_long)
+ stop
+#endif
+ endif
+! end change
+c lprint_short=.false.
if (large.and. mod(itime,ntwe).eq.0)
& call enerprint(energia_long)
#ifdef TIMING_ENE
write (iout,*) "energia_long",energia_long(0)
write (iout,*) "Cartesian and internal coordinates: step 2"
c call cartprint
- call pdbout(0.0d0,"cipiszcze",iout)
+ call pdbout(0.0d0,
+ & cipiszcze ,iout)
call intout
write (iout,*) "Accelerations from long-range forces"
do i=0,nres
potEcomp(i)=energia_short(i)+energia_long(i)
enddo
potE=potEcomp(0)-potEcomp(20)
+ if (ntwe.ne.0) then
+ if (large.and. mod(itime,ntwe).eq.0) then
+ call enerprint(potEcomp)
+ write (iout,*) "potE",potD
+ endif
+ endif
c potE=energia_short(0)+energia_long(0)
totT=totT+d_time
+ totTafm=totT
c Calculate the kinetic and the total energy and the kinetic temperature
call kinetic(EK)
totE=EK+potE
enddo
enddo
do i=nnt,nct
- if (itype(i).ne.10 .and. itype(i).ne.21) then
+ if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
inres=i+nres
do j=1,3
d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
d_t(j,0)=d_t_old(j,0)+adt
enddo
do i=nnt,nct-1
+C SPYTAC ADAMA
+C do i=0,nres
do j=1,3
adt=d_a_old(j,i)*d_time
adt2=0.5d0*adt
enddo
enddo
do i=nnt,nct
- if (itype(i).ne.10 .and. itype(i).ne.21) then
+C do i=0,nres
+ if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
inres=i+nres
do j=1,3
adt=d_a_old(j,inres)*d_time
enddo
enddo
do i=nnt,nct
- if (itype(i).ne.10 .and. itype(i).ne.21) then
+ if (itype(i).ne.10 .and. itype(i).ne.ntyp1) 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
c
c Compute friction and stochastic forces
c
+#ifdef MPI
time00=MPI_Wtime()
+#else
+ time00=tcpu()
+#endif
call friction_force
+#ifdef MPI
time_fric=time_fric+MPI_Wtime()-time00
time00=MPI_Wtime()
+#else
+ time_fric=time_fric+tcpu()-time00
+ time00=tcpu()
+#endif
call stochastic_force(stochforcvec)
+#ifdef MPI
time_stoch=time_stoch+MPI_Wtime()-time00
+#else
+ time_stoch=time_stoch+tcpu()-time00
+#endif
c
c Compute the acceleration due to friction forces (d_af_work) and stochastic
c forces (d_as_work)
ind=ind+3
enddo
do i=nnt,nct
- if (itype(i).ne.10 .and. itype(i).ne.21) then
+ if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
inres=i+nres
do j=1,3
adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
ind=ind+3
enddo
do i=nnt,nct
- if (itype(i).ne.10 .and. itype(i).ne.21) then
+ if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
inres=i+nres
do j=1,3
d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres)
c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
if (dabs(accel(j)).gt.dabs(accel_old(j))) then
dacc=dabs(accel(j)-accel_old(j))
+c write (iout,*) i,dacc
if (dacc.gt.amax) amax=dacc
endif
enddo
enddo
endif
do i=nnt,nct
- if (itype(i).ne.10 .and. itype(i).ne.21) then
+ if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
do j=1,3
c 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)
c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
if (dabs(accel(j)).gt.dabs(accel_old(j))) then
dacc=dabs(accel(j)-accel_old(j))
+c write (iout,*) "side-chain",i,dacc
if (dacc.gt.amax) amax=dacc
endif
enddo
enddo
endif
c Side chains
- if (itype(i).ne.10 .and. itype(i).ne.21) then
+ if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
do j=1,3
epdriftij=
& dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
enddo
enddo
do i=nnt,nct
- if (itype(i).ne.10 .and. itype(i).ne.21) then
+ if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
inres=i+nres
do j=1,3
d_t(j,inres)=fact*d_t(j,inres)
stdforcp(i)=stdfp*dsqrt(gamp)
enddo
do i=nnt,nct
- stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamsc(itype(i)))
+ stdforcsc(i)=stdfsc(iabs(itype(i)))
+ & *dsqrt(gamsc(iabs(itype(i))))
enddo
endif
c Open the pdb file for snapshotshots
if (restart1file) then
if (me.eq.king)
& inquire(file=mremd_rst_name,exist=file_exist)
+#ifdef MPI
write (*,*) me," Before broadcast: file_exist",file_exist
call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
& IERR)
write (*,*) me," After broadcast: file_exist",file_exist
c inquire(file=mremd_rst_name,exist=file_exist)
+#endif
if(me.eq.king.or..not.out1file)
& write(iout,*) "Initial state read by master and distributed"
else
endif
call rescale_weights(t_bath)
else
+ rest=.false.
if(me.eq.king.or..not.out1file)then
if (restart1file) then
write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),
endif
call random_vel
totT=0.0d0
+ totTafm=totT
endif
else
c Generate initial velocities
& write(iout,*) "Initial velocities randomly generated"
call random_vel
totT=0.0d0
+CtotTafm is the variable for AFM time which eclipsed during
+ totTafm=totT
endif
c rest2name = prefix(:ilen(prefix))//'.rst'
if(me.eq.king.or..not.out1file)then
c Zeroing the total angular momentum of the system
write(iout,*) "Calling the zero-angular
& momentum subroutine"
+ call flush(iout)
endif
call inertia_tensor
+c write (iout,*) "After inertia"
+c call flush(iout)
c Getting the potential energy and forces and velocities and accelerations
+ if(me.eq.king.or..not.out1file)then
+ write(iout,*) "Calling the vcm_vel"
+ call flush(iout)
+ endif
call vcm_vel(vcm)
-c write (iout,*) "velocity of the center of the mass:"
-c write (iout,*) (vcm(j),j=1,3)
+ write (iout,*) "velocity of the center of the mass:"
+ write (iout,*) (vcm(j),j=1,3)
+ call flush(iout)
do j=1,3
d_t(j,0)=d_t(j,0)-vcm(j)
enddo
c Removing the velocity of the center of mass
+ if(me.eq.king.or..not.out1file)then
+ write(iout,*) "Calling the vcm_vel"
+ call flush(iout)
+ endif
call vcm_vel(vcm)
if(me.eq.king.or..not.out1file)then
write (iout,*) "vcm right after adjustment:"
write (iout,*) (vcm(j),j=1,3)
+ call flush(iout)
endif
if (.not.rest) then
- call chainbuild
if(iranconf.ne.0) then
+c 8/22/17 AL Loop to produce a low-energy random conformation
+ do iranmin=1,10
+ write (iout,*) "iranmin",iranmin
+ call chainbuild
if (overlapsc) then
print *, 'Calling OVERLAP_SC'
call overlap_sc(fail)
endif
if(me.eq.king.or..not.out1file)
& write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
+ if (isnan(etot) .or. etot.gt.1.0d4) then
+ write (iout,*) "Energy too large",etot,
+ & " trying another random conformation"
+ do itrial=1,100
+ itmp=1
+ nrestmp=nres
+ call gen_rand_conf(itmp,nrestmp,*30)
+ goto 40
+ 30 write (iout,*) 'Failed to generate random conformation',
+ & ', itrial=',itrial
+ write (*,*) 'Processor:',me,
+ & ' Failed to generate random conformation',
+ & ' itrial=',itrial
+ call intout
+
+#ifdef AIX
+ call flush_(iout)
+#else
+ call flush(iout)
+#endif
+ enddo
+ write (iout,'(a,i3,a)') 'Processor:',me,
+ & ' error in generating random conformation.'
+ write (*,'(a,i3,a)') 'Processor:',me,
+ & ' error in generating random conformation.'
+ call flush(iout)
+#ifdef MPI
+ call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
+#else
+ stop
+#endif
+ 40 continue
+ else
+ goto 44
+ endif
+ enddo
+ write (iout,'(a,i3,a)') 'Processor:',me,
+ & ' failed to generate a low-energy random conformation.'
+ write (*,'(a,i3,a)') 'Processor:',me,
+ & ' failed to generate a low-energy random conformation.'
+ call flush(iout)
+#ifdef MPI
+ call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
+#else
+ stop
+#endif
+ 44 continue
+ else if (indpdb.gt.0) then
+C 8/22/17 AL Minimize initial PDB structure
+ if(dccart)then
+ print *, 'Calling MINIM_DC'
+ call minim_dc(etot,iretcode,nfun)
+ else
+ call geom_to_var(nvar,varia)
+ print *,'Calling MINIMIZE.'
+ call minimize(etot,varia,iretcode,nfun)
+ call var_to_geom(nvar,varia)
+ endif
+ if(me.eq.king.or..not.out1file)
+ & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
endif
endif
call chainbuild_cart
call kinetic(EK)
if (tbf) then
- call verlet_bath(EK)
+ call verlet_bath
endif
kinetic_T=2.0d0/(dimen3*Rb)*EK
if(me.eq.king.or..not.out1file)then
if(me.eq.king.or..not.out1file)then
write(iout,*) "Potential energy and its components"
call enerprint(potEcomp)
-c write(iout,*) (potEcomp(i),i=0,n_ene)
+ write(iout,*) (potEcomp(i),i=0,n_ene)
endif
potE=potEcomp(0)-potEcomp(20)
totE=EK+potE
include 'COMMON.IOUNITS'
include 'COMMON.NAMES'
include 'COMMON.TIME1'
- double precision xv,sigv,lowb,highb
+ double precision xv,sigv,lowb,highb,vec_afm(3)
c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
c First generate velocities in the eigenspace of the G matrix
c 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)
+
c write (iout,*) "i",i," ii",ii," geigen",geigen(i),
c & " d_t_work_new",d_t_work_new(ii)
enddo
enddo
+C if (SELFGUIDE.gt.0) then
+C distance=0.0
+C do j=1,3
+C vec_afm(j)=c(j,afmend)-c(j,afmbeg)
+C distance=distance+vec_afm(j)**2
+C enddo
+C distance=dsqrt(distance)
+C do j=1,3
+C d_t_work_new(j+(afmbeg-1)*3)=-velAFMconst*vec_afm(j)/distance
+C d_t_work_new(j+(afmend-1)*3)=velAFMconst*vec_afm(j)/distance
+C write(iout,*) "myvel",d_t_work_new(j+(afmbeg-1)*3),
+C & d_t_work_new(j+(afmend-1)*3)
+C enddo
+
+C endif
+
c diagnostics
c Ek1=0.0d0
c ii=0
enddo
enddo
do i=nnt,nct
- if (itype(i).ne.10 .and. itype(i).ne.21) then
+ if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
do j=1,3
ind=ind+1
d_t(j,i+nres)=d_t_work(ind)
ind=ind+3
enddo
do i=nnt,nct
- if (itype(i).ne.10 .and. itype(i).ne.21) then
+ if (itype(i).ne.10 .and. itype(i).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.21) then
+ if (itype(i).ne.10 .and. itype(i).ne.ntyp1) 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.21) then
+ if (itype(i).ne.10 .and. itype(i).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.21) then
+ if (itype(i).ne.10 .and. itype(i).ne.ntyp1) 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.21) then
+ if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
inres=i+nres
do j=1,3
d_t(j,inres)=d_t_work(ind+j)
projects / unres.git / blobdiff