X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2FMD_A-MTS.F;h=84c5a187e2de8f95f9b63a4eeb0d8210afb12fa3;hb=3a2758de24121fb0f72bd604e2163b7a4b73a27e;hp=56b3ea89b487f9e074a1044bab92074b2bbf405c;hpb=af72f8e89a5d33f0d86ba898d6c5bbbda4b25b84;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 56b3ea8..84c5a18 100644 --- a/source/unres/src_MD-M/MD_A-MTS.F +++ b/source/unres/src_MD-M/MD_A-MTS.F @@ -36,6 +36,7 @@ c------------------------------------------------ common /gucio/ cm integer itime logical ovrtim + integer nharp,iharp(4,maxres/3) c #ifdef MPI if (ilen(tmpdir).gt.0) @@ -51,6 +52,33 @@ c t_enegrad=0.0d0 t_sdsetup=0.0d0 write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started" + write (iout,'(/a)') + & "Cartesian coordinates of the initial structure" + write (iout,'(a,3(3x,a5),5x,3(3x,a5))') + & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)" + do ires=1,nres + write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)') + & restyp(itype(ires)),ires,(c(j,ires),j=1,3), + & (c(j,ires+nres),j=1,3) + enddo + write (iout,'(/a)') + & "Initial dC vectors of the chain" + write (iout,'(a,3(3x,a5),5x,3(3x,a5))') + & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)" + do ires=1,nres + write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)') + & restyp(itype(ires)),ires,(dc(j,ires),j=1,3), + & (dc(j,ires+nres),j=1,3) + enddo + write (iout,'(/a)') + & "Initial dC_norm vectors of the chain" + write (iout,'(a,3(3x,a5),5x,3(3x,a5))') + & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)" + do ires=1,nres + write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)') + & restyp(itype(ires)),ires,(dc_norm(j,ires),j=1,3), + & (dc_norm(j,ires+nres),j=1,3) + enddo #ifdef MPI tt0=MPI_Wtime() #else @@ -196,7 +224,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 +241,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 +262,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 +548,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 +947,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 +993,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 +1034,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 +1045,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 +1090,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 @@ -1161,7 +1202,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 @@ -1226,7 +1267,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) @@ -1273,6 +1314,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 @@ -1291,7 +1333,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) @@ -1302,6 +1344,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 @@ -1344,7 +1387,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)) @@ -1391,7 +1434,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) @@ -1428,7 +1471,7 @@ c Set up the initial conditions of a MD simulation include 'COMMON.NAMES' include 'COMMON.REMD' real*8 energia_long(0:n_ene), - & energia_short(0:n_ene),vcm(3),incr(3) + & energia_short(0:n_ene),energia(0:n_ene),vcm(3),incr(3) double precision cm(3),L(3),xv,sigv,lowb,highb double precision varia(maxvar) character*256 qstr @@ -1446,7 +1489,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 @@ -1556,6 +1600,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 @@ -1563,6 +1608,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 @@ -1590,10 +1637,38 @@ c Removing the velocity of the center of mass write (iout,*) (vcm(j),j=1,3) endif if (.not.rest) then - call chainbuild - if(iranconf.ne.0) then + if (start_from_model) then + i_model=iran_num(1,constr_homology) + write (iout,*) 'starting from model ',i_model + do i=1,2*nres + do j=1,3 + c(j,i)=chomo(j,i,i_model) + enddo + enddo + + call int_from_cart(.true.,.false.) + call sc_loc_geom(.false.) + do i=1,nres-1 + do j=1,3 + dc(j,i)=c(j,i+1)-c(j,i) + dc_norm(j,i)=dc(j,i)*vbld_inv(i+1) + enddo + enddo + do i=2,nres-1 + do j=1,3 + dc(j,i+nres)=c(j,i+nres)-c(j,i) + dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres) + enddo + enddo + call etotal(energia(0)) + endif +c call chainbuild_cart + write (iout,*) "Initial energies" + call enerprint(energia(0)) + write (iout,*) "PREMINIM ",preminim + if(iranconf.ne.0 .or. preminim) then if (overlapsc) then - print *, 'Calling OVERLAP_SC' + write (iout,*) 'Calling OVERLAP_SC' call overlap_sc(fail) endif @@ -1613,8 +1688,87 @@ c Removing the velocity of the center of mass call minimize(etot,varia,iretcode,nfun) call var_to_geom(nvar,varia) endif + if(me.eq.king.or..not.out1file) then + write(iout,*) "Minimized energy is",etot + write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun + call etotal(potEcomp) + call enerprint(potEcomp) + endif + endif + if(iranconf.ne.0) then +c 8/22/17 AL Loop to produce a low-energy random conformation + DO iranmin=1,10 + if (overlapsc) then + if(me.eq.king.or..not.out1file) + & write (iout,*) 'Calling OVERLAP_SC' + call overlap_sc(fail) + endif + + if (searchsc) then + call sc_move(2,nres-1,10,1d10,nft_sc,etot) + print *,'SC_move',nft_sc,etot + if(me.eq.king.or..not.out1file) + & write(iout,*) 'SC_move',nft_sc,etot + endif + + 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 + + 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 + call gen_rand_conf(itmp,*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 endif endif call chainbuild_cart @@ -1785,7 +1939,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 @@ -1799,10 +1953,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 @@ -1839,7 +2010,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) @@ -2068,7 +2239,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) @@ -2177,7 +2348,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) @@ -2334,7 +2505,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) @@ -2383,7 +2554,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) @@ -2444,7 +2615,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)