#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)
#endif
endif
if (mod(itime,ntwx).eq.0) then
+ write(iout,*) 'time=',itime
+C call check_ecartint
+ call returnbox
write (tytul,'("time",f8.2)') totT
if(mdpdb) then
call hairpin(.true.,nharp,iharp)
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
potE=potEcomp(0)-potEcomp(20)
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
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
+C do i=0,nres
if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
inres=i+nres
do j=1,3
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)
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))
- write (iout,*) i,dacc
+c write (iout,*) i,dacc
if (dacc.gt.amax) amax=dacc
endif
enddo
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))
- write (iout,*) "side-chain",i,dacc
+c write (iout,*) "side-chain",i,dacc
if (dacc.gt.amax) amax=dacc
endif
enddo
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 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
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
#endif
call zerograd
call etotal(potEcomp)
+ call enerprint(potEcomp)
if (large) call enerprint(potEcomp)
#ifdef TIMING_ENE
#ifdef MPI
& "Time step reduced to",d_time,
& " because of too large initial acceleration."
endif
- if(me.eq.king.or..not.out1file)then
- write(iout,*) "Potential energy and its components"
- call enerprint(potEcomp)
+C if(me.eq.king.or..not.out1file)then
+C write(iout,*) "Potential energy and its components"
+C call enerprint(potEcomp)
c write(iout,*) (potEcomp(i),i=0,n_ene)
- endif
+C endif
potE=potEcomp(0)-potEcomp(20)
totE=EK+potE
itime=0
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