X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fio.f90;h=379978e31142cd3ddc053513c50d58113d00eb46;hb=8d2ab9ba185dbbc31bfe3d1c66d7e1c9d632463b;hp=243c8b62acd3bb43e9411004a82352433f27afad;hpb=05fce52032feb59f7e96cd897fb82ca2aa90a888;p=unres4.git diff --git a/source/unres/io.f90 b/source/unres/io.f90 index 243c8b6..379978e 100644 --- a/source/unres/io.f90 +++ b/source/unres/io.f90 @@ -487,7 +487,7 @@ subroutine statout(itime) use energy_data - use control_data, only:refstr + use control_data use MD_data use MPI_data use compare, only:rms_nac_nnc @@ -513,8 +513,8 @@ character(len=4) :: format1,format2 character(len=30) :: format !el local variables - integer :: i,ii1,ii2 - real(kind=8) :: rms,frac,frac_nn,co + integer :: i,ii1,ii2,j + real(kind=8) :: rms,frac,frac_nn,co,distance #ifdef AIX if(itime.eq.0) then @@ -527,6 +527,46 @@ open(istat,file=statname,access="append") #endif #endif + if (AFMlog.gt.0) then + if (refstr) then + call rms_nac_nnc(rms,frac,frac_nn,co,.false.) + write (line1,'(i10,f15.2,3f12.3,f7.2,2f6.3,4f12.3,i5,$)')& + itime,totT,EK,potE,totE,& + rms,frac,frac_nn,kinetic_T,t_bath,gyrate(),& + potEcomp(23),me + format1="a133" + else +!C print *,'A CHUJ',potEcomp(23) + write (line1,'(i10,f15.2,7f12.3,i5,$)') & + itime,totT,EK,potE,totE,& + kinetic_T,t_bath,gyrate(),& + potEcomp(23),me + format1="a114" + endif + else if (selfguide.gt.0) then + distance=0.0 + do j=1,3 + distance=distance+(c(j,afmend)-c(j,afmbeg))**2 + enddo + distance=dsqrt(distance) + if (refstr) then + call rms_nac_nnc(rms,frac,frac_nn,co,.false.) + write (line1,'(i10,f15.2,3f12.3,f7.2,2f6.3,f12.3,f10.1,2f8.2, & + f9.3,i5,$)') & + itime,totT,EK,potE,totE,& + rms,frac,frac_nn,kinetic_T,t_bath,gyrate(),& + distance,potEcomp(23),me + format1="a133" +!C print *,"CHUJOWO" + else +!C print *,'A CHUJ',potEcomp(23) + write (line1,'(i10,f15.2,8f12.3,i5,$)')& + itime,totT,EK,potE,totE, & + kinetic_T,t_bath,gyrate(),& + distance,potEcomp(23),me + format1="a114" + endif + else if (refstr) then call rms_nac_nnc(rms,frac,frac_nn,co,.false.) write (line1,'(i10,f15.2,3f12.3,f7.2,4f6.3,3f12.3,i5,$)') & @@ -539,6 +579,7 @@ amax,kinetic_T,t_bath,gyrate(),me format1="a114" endif + ENDIF if(usampl.and.totT.gt.eq_time) then write(line2,'(i5,2f9.4,300f7.4)') iset,uconst,uconst_back,& (qfrag(ii1),ii1=1,nfrag),(qpair(ii2),ii2=1,npair),& @@ -695,7 +736,7 @@ integer,dimension(maxres) :: itype_alloc integer :: iti,nsi,maxsi,itrial,itmp - real(kind=8) :: wlong,scalscp,co + real(kind=8) :: wlong,scalscp,co,ssscale allocate(weights(n_ene)) !----------------------------- allocate(c(3,2*maxres+2)) !(3,maxres2+2) maxres2=2*maxres @@ -730,6 +771,10 @@ call reada(weightcard,'WBOND',wbond,1.0D0) call reada(weightcard,'WTOR',wtor,1.0D0) call reada(weightcard,'WTORD',wtor_d,1.0D0) + call reada(weightcard,'WSHIELD',wshield,0.05D0) + call reada(weightcard,'LIPSCALE',lipscale,1.0D0) + call reada(weightcard,'WLT',wliptran,1.0D0) + call reada(weightcard,'WTUBE',wtube,1.0d0) call reada(weightcard,'WANG',wang,1.0D0) call reada(weightcard,'WSCLOC',wscloc,1.0D0) call reada(weightcard,'SCAL14',scal14,0.4D0) @@ -841,18 +886,23 @@ call reada(weightcard,"V2SS",v2ss,7.61d0) call reada(weightcard,"V3SS",v3ss,13.7d0) call reada(weightcard,"EBR",ebr,-5.50D0) + call reada(weightcard,"ATRISS",atriss,0.301D0) + call reada(weightcard,"BTRISS",btriss,0.021D0) + call reada(weightcard,"CTRISS",ctriss,1.001D0) + call reada(weightcard,"DTRISS",dtriss,1.001D0) + call reada(weightcard,"SSSCALE",ssscale,1.0D0) dyn_ss=(index(weightcard,'DYN_SS').gt.0) call reada(weightcard,"HT",Ht,0.0D0) if (dyn_ss) then - ss_depth=ebr/wsc-0.25*eps(1,1) - Ht=Ht/wsc-0.25*eps(1,1) - akcm=akcm*wstrain/wsc - akth=akth*wstrain/wsc - akct=akct*wstrain/wsc - v1ss=v1ss*wstrain/wsc - v2ss=v2ss*wstrain/wsc - v3ss=v3ss*wstrain/wsc + ss_depth=(ebr/wsc-0.25*eps(1,1))*ssscale + Ht=(Ht/wsc-0.25*eps(1,1))*ssscale + akcm=akcm/wsc*ssscale + akth=akth/wsc*ssscale + akct=akct/wsc*ssscale + v1ss=v1ss/wsc*ssscale + v2ss=v2ss/wsc*ssscale + v3ss=v3ss/wsc*ssscale else ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain endif @@ -1019,6 +1069,7 @@ if (ndih_constr.gt.0) then allocate(idih_constr(ndih_constr),idih_nconstr(ndih_constr)) !(maxdih_constr) allocate(phi0(ndih_constr),drange(ndih_constr)) !(maxdih_constr) + read (inp,*) ftors read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr) if(me.eq.king.or..not.out1file)then @@ -1040,6 +1091,47 @@ phibound(2,ii) = phi0(i)+drange(i) enddo endif + if (with_theta_constr) then +!C with_theta_constr is keyword allowing for occurance of theta constrains + read (inp,*) ntheta_constr +!C ntheta_constr is the number of theta constrains + if (ntheta_constr.gt.0) then +!C read (inp,*) ftors + allocate(itheta_constr(ntheta_constr)) + allocate(theta_constr0(ntheta_constr)) + allocate(theta_drange(ntheta_constr),for_thet_constr(ntheta_constr)) + read (inp,*) (itheta_constr(i),theta_constr0(i), & + theta_drange(i),for_thet_constr(i), & + i=1,ntheta_constr) +!C the above code reads from 1 to ntheta_constr +!C itheta_constr(i) residue i for which is theta_constr +!C theta_constr0 the global minimum value +!C theta_drange is range for which there is no energy penalty +!C for_thet_constr is the force constant for quartic energy penalty +!C E=k*x**4 + if(me.eq.king.or..not.out1file)then + write (iout,*) & + 'There are',ntheta_constr,' constraints on phi angles.' + do i=1,ntheta_constr + write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i), & + theta_drange(i), & + for_thet_constr(i) + enddo + endif + do i=1,ntheta_constr + theta_constr0(i)=deg2rad*theta_constr0(i) + theta_drange(i)=deg2rad*theta_drange(i) + enddo +!C if(me.eq.king.or..not.out1file) +!C & write (iout,*) 'FTORS',ftors +!C do i=1,ntheta_constr +!C ii = itheta_constr(i) +!C thetabound(1,ii) = phi0(i)-drange(i) +!C thetabound(2,ii) = phi0(i)+drange(i) +!C enddo + endif ! ntheta_constr.gt.0 + endif! with_theta_constr + nnt=1 #ifdef MPI if (me.eq.king) then @@ -1130,6 +1222,7 @@ call flush(iout) if (constr_dist.gt.0) call read_dist_constr write (iout,*) "After read_dist_constr nhpb",nhpb + if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp call hpb_partition if(me.eq.king.or..not.out1file) & write (iout,*) 'Contact order:',co