X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2FMD_A-MTS.F;h=2635a2f9fc586846b39e11946fa2a0af81f69fbb;hb=d96a47c690dabd9c300b08916cd2811a37c3e5eb;hp=0fe17b077a108be8493c5dadc17014878e044154;hpb=2d221fe0481d0fa3c71cf9f72d6b89f553d55a5c;p=unres.git diff --git a/source/unres/src_MD-M/MD_A-MTS.F b/source/unres/src_MD-M/MD_A-MTS.F index 0fe17b0..2635a2f 100644 --- a/source/unres/src_MD-M/MD_A-MTS.F +++ b/source/unres/src_MD-M/MD_A-MTS.F @@ -196,7 +196,11 @@ c Variable time step algorithm. #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) @@ -230,6 +234,9 @@ c Variable time step algorithm. #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) @@ -513,6 +520,8 @@ c Second step of the velocity Verlet algorithm 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 @@ -910,6 +919,7 @@ c Compute the complete potential energy 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 @@ -996,6 +1006,8 @@ c Applying velocity Verlet algorithm - step 1 to coordinates 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 @@ -1005,6 +1017,7 @@ c Applying velocity Verlet algorithm - step 1 to coordinates 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 @@ -1087,12 +1100,25 @@ c Applying velocity Verlet algorithm - step 1 to coordinates 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) @@ -1510,11 +1536,13 @@ c statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat' 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 @@ -1544,6 +1572,7 @@ c inquire(file=mremd_rst_name,exist=file_exist) endif call random_vel totT=0.0d0 + totTafm=totT endif else c Generate initial velocities @@ -1551,6 +1580,8 @@ 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 @@ -1608,7 +1639,7 @@ c Removing the velocity of the center of mass 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 @@ -1640,11 +1671,11 @@ c Removing the velocity of the center of mass & "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 @@ -1773,7 +1804,7 @@ c----------------------------------------------------------- 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 @@ -1787,10 +1818,27 @@ c call flush(iout) 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