X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2FMD_A-MTS.F;h=3a452319fa8326d49b9629c158d9e35ba169ed6f;hb=14b55af23c7da34bfdc10f315551c8e01a3a906a;hp=d3c3cb0f025943c6d805ccc08c32bd13175f6900;hpb=395aa58d0681d84b591ddaab57bbee4b9e9d5dad;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 d3c3cb0..3a45231 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) @@ -209,7 +213,7 @@ c Variable time step algorithm. 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) @@ -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 @@ -955,7 +965,7 @@ c forces). 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 @@ -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,7 +1017,8 @@ c Applying velocity Verlet algorithm - step 1 to coordinates 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 @@ -1049,7 +1062,7 @@ c Step 2 of the velocity Verlet algorithm: update velocities 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 @@ -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) @@ -1148,7 +1174,7 @@ c 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 @@ -1213,7 +1239,7 @@ c 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) @@ -1260,6 +1286,7 @@ c accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i)) 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 @@ -1278,7 +1305,7 @@ c accel(j)=aux(j) 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) @@ -1289,6 +1316,7 @@ c accel(j)=accel(j)+d_a(j,i+nres)-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 @@ -1331,7 +1359,7 @@ c write (iout,*) "back",i,j,epdriftij 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)) @@ -1378,7 +1406,7 @@ c write(iout,*) "fact", fact 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) @@ -1433,7 +1461,8 @@ c if the friction coefficients do not depend on surface area 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 @@ -1507,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 @@ -1541,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 @@ -1548,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 @@ -1619,6 +1653,7 @@ c Removing the velocity of the center of mass #endif call zerograd call etotal(potEcomp) + call enerprint(potEcomp) if (large) call enerprint(potEcomp) #ifdef TIMING_ENE #ifdef MPI @@ -1637,11 +1672,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 @@ -1770,7 +1805,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 @@ -1784,10 +1819,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 @@ -1824,7 +1876,7 @@ c Transfer to the d_t vector 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) @@ -2053,7 +2105,7 @@ c enddo 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) @@ -2162,7 +2214,7 @@ c 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) @@ -2319,7 +2371,7 @@ c enddo 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) @@ -2368,7 +2420,7 @@ c enddo 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) @@ -2429,7 +2481,7 @@ c ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(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)