From: Adam Sieradzan Date: Mon, 3 Aug 2015 19:39:58 +0000 (+0200) Subject: 1) debug of Lorentzian like constrains in cluster 2) introduction of constrains on... X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?p=unres.git;a=commitdiff_plain;h=29c2a95f7c3afbc83690ae6148c805a3572dc3fd 1) debug of Lorentzian like constrains in cluster 2) introduction of constrains on theta angle --- diff --git a/source/cluster/wham/src-M/COMMON.CONTROL b/source/cluster/wham/src-M/COMMON.CONTROL index e65a5b2..f339ae9 100644 --- a/source/cluster/wham/src-M/COMMON.CONTROL +++ b/source/cluster/wham/src-M/COMMON.CONTROL @@ -3,9 +3,9 @@ & constr_dist logical refstr,pdbref,punch_dist,print_dist,caonly,lside, & lprint_cart,lprint_int,from_cart,efree,from_bx,from_cx, - & with_dihed_constr + & with_dihed_constr,with_theta_constr common /cntrl/ betaT,iscode,indpdb,refstr,pdbref,outpdb,outmol2, & punch_dist,print_dist,caonly,lside,lprint_cart,lprint_int, - & from_cart,from_bx,from_cx, with_dihed_constr, + & from_cart,from_bx,from_cx, with_dihed_constr,with_theta_constr, & efree,iopt,nstart,nend,symetr, & constr_dist diff --git a/source/cluster/wham/src-M/DIMENSIONS b/source/cluster/wham/src-M/DIMENSIONS index 68ac673..29e8f24 100644 --- a/source/cluster/wham/src-M/DIMENSIONS +++ b/source/cluster/wham/src-M/DIMENSIONS @@ -60,7 +60,7 @@ C Max number of symetric chains parameter (maxperm=120) C Max. number of energy components integer max_ene - parameter (max_ene=21) + parameter (max_ene=24) C Max. number of temperatures integer maxt parameter (maxT=5) diff --git a/source/cluster/wham/src-M/energy_p_new.F b/source/cluster/wham/src-M/energy_p_new.F index 39d35dd..f78e2e9 100644 --- a/source/cluster/wham/src-M/energy_p_new.F +++ b/source/cluster/wham/src-M/energy_p_new.F @@ -67,7 +67,7 @@ cd print *,'EHPB exitted succesfully.' C C Calculate the virtual-bond-angle energy. C - call ebend(ebe) + call ebend(ebe,ethetacnstr) cd print *,'Bend energy finished.' C C Calculate the SC local energy. @@ -111,7 +111,7 @@ c write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t & +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4 & +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6 & +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d - & +wbond*estr+wsccor*fact(1)*esccor + & +wbond*estr+wsccor*fact(1)*esccor+ethetacnstr #else etot=wsc*(evdw+fact(6)*evdw_t)+wscp*evdw2 & +welec*fact(1)*(ees+evdw1) @@ -120,7 +120,7 @@ c write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t & +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4 & +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6 & +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d - & +wbond*estr+wsccor*fact(1)*esccor + & +wbond*estr+wsccor*fact(1)*esccor+ethetacnstr #endif energia(0)=etot energia(1)=evdw @@ -154,6 +154,7 @@ c write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t energia(19)=esccor energia(20)=edihcnstr energia(21)=evdw_t + energia(24)=ethetacnstr c detecting NaNQ #ifdef ISNAN #ifdef AIX @@ -270,6 +271,7 @@ C------------------------------------------------------------------------ esccor=energia(19) edihcnstr=energia(20) estr=energia(18) + ethetacnstr=energia(24) #ifdef SPLITELE write (iout,10) evdw,wsc,evdw2,wscp,ees,welec*fact(1),evdw1, & wvdwpp, @@ -278,7 +280,7 @@ C------------------------------------------------------------------------ & ecorr,wcorr*fact(3),ecorr5,wcorr5*fact(4),ecorr6,wcorr6*fact(5), & eel_loc,wel_loc*fact(2),eello_turn3,wturn3*fact(2), & eello_turn4,wturn4*fact(3),eello_turn6,wturn6*fact(5), - & esccor,wsccor*fact(1),edihcnstr,ebr*nss,etot + & esccor,wsccor*fact(1),edihcnstr,ethetacnstr,ebr*nss,etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ @@ -300,6 +302,7 @@ C------------------------------------------------------------------------ & 'ETURN6=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 6th order)'/ & 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/ & 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/ + & 'ETHETC= ',1pE16.6,' (valence angle constraints)'/ & 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ & 'ETOT= ',1pE16.6,' (total)') #else @@ -309,7 +312,7 @@ C------------------------------------------------------------------------ & ecorr6,wcorr6*fact(5),eel_loc,wel_loc*fact(2), & eello_turn3,wturn3*fact(2),eello_turn4,wturn4*fact(3), & eello_turn6,wturn6*fact(5),esccor*fact(1),wsccor, - & edihcnstr,ebr*nss,etot + & edihcnstr,ethetacnstr,ebr*nss,etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ @@ -330,6 +333,7 @@ C------------------------------------------------------------------------ & 'ETURN6=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 6th order)'/ & 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/ & 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/ + & 'ETHETC= ',1pE16.6,' (valence angle constraints)'/ & 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ & 'ETOT= ',1pE16.6,' (total)') #endif @@ -2909,6 +2913,7 @@ C include 'COMMON.DERIV' include 'COMMON.VAR' include 'COMMON.INTERACT' + include 'COMMON.CONTROL' dimension ggg(3) ehpb=0.0D0 cd print *,'edis: nhpb=',nhpb,' fbr=',fbr @@ -2952,6 +2957,7 @@ C & *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i)) & *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd C write (iout,'(a6,2i5,3f8.3)') "edisl",ii,jj, C & ehpb,fordepth(i),dd +C print *,"TUTU" C write(iout,*) ehpb,"atu?" C ehpb,"tu?" C fac=fordepth(i)**4.0d0 @@ -3228,7 +3234,7 @@ c & AKSC(j,iti),abond0(j,iti),u(j),j=1,nbi) end #ifdef CRYST_THETA C-------------------------------------------------------------------------- - subroutine ebend(etheta) + subroutine ebend(etheta,ethetacnstr) C C Evaluate the virtual-bond-angle energy given the virtual-bond dihedral C angles gamma and its derivatives in consecutive thetas and gammas. @@ -3245,6 +3251,7 @@ C include 'COMMON.IOUNITS' include 'COMMON.NAMES' include 'COMMON.FFIELD' + include 'COMMON.TORCNSTR' common /calcthet/ term1,term2,termm,diffak,ratak, & ak,aktc,termpre,termexp,sigc,sig0i,time11,time12,sigcsq, & delthe0,sig0inv,sigtc,sigsqtc,delthec,it @@ -3363,6 +3370,34 @@ c & rad2deg*phii,rad2deg*phii1,ethetai c 1215 continue enddo C Ufff.... We've done all this!!! +C now constrains + ethetacnstr=0.0d0 +C print *,ithetaconstr_start,ithetaconstr_end,"TU" + do i=1,ntheta_constr + itheta=itheta_constr(i) + thetiii=theta(itheta) + difi=pinorm(thetiii-theta_constr0(i)) + if (difi.gt.theta_drange(i)) then + difi=difi-theta_drange(i) + ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4 + gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) + & +for_thet_constr(i)*difi**3 + else if (difi.lt.-drange(i)) then + difi=difi+drange(i) + ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4 + gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) + & +for_thet_constr(i)*difi**3 + else + difi=0.0 + endif +C if (energy_dec) then +C write (iout,'(a6,2i5,4f8.3,2e14.5)') "ethetc", +C & i,itheta,rad2deg*thetiii, +C & rad2deg*theta_constr0(i), rad2deg*theta_drange(i), +C & rad2deg*difi,0.25d0*for_thet_constr(i)*difi**4, +C & gloc(itheta+nphi-2,icg) +C endif + enddo return end C--------------------------------------------------------------------------- @@ -3475,7 +3510,7 @@ C "Thank you" to MAPLE (probably spared one day of hand-differentiation). end #else C-------------------------------------------------------------------------- - subroutine ebend(etheta) + subroutine ebend(etheta,ethetacnstr) C C Evaluate the virtual-bond-angle energy given the virtual-bond dihedral C angles gamma and its derivatives in consecutive thetas and gammas. @@ -3495,6 +3530,7 @@ C include 'COMMON.NAMES' include 'COMMON.FFIELD' include 'COMMON.CONTROL' + include 'COMMON.TORCNSTR' double precision coskt(mmaxtheterm),sinkt(mmaxtheterm), & cosph1(maxsingle),sinph1(maxsingle),cosph2(maxsingle), & sinph2(maxsingle),cosph1ph2(maxdouble,maxdouble), @@ -3675,6 +3711,34 @@ c call flush(iout) c gloc(nphi+i-2,icg)=wang*dethetai gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wang*dethetai enddo +C now constrains + ethetacnstr=0.0d0 +C print *,ithetaconstr_start,ithetaconstr_end,"TU" + do i=1,ntheta_constr + itheta=itheta_constr(i) + thetiii=theta(itheta) + difi=pinorm(thetiii-theta_constr0(i)) + if (difi.gt.theta_drange(i)) then + difi=difi-theta_drange(i) + ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4 + gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) + & +for_thet_constr(i)*difi**3 + else if (difi.lt.-drange(i)) then + difi=difi+drange(i) + ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4 + gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) + & +for_thet_constr(i)*difi**3 + else + difi=0.0 + endif +C if (energy_dec) then +C write (iout,'(a6,2i5,4f8.3,2e14.5)') "ethetc", +C & i,itheta,rad2deg*thetiii, +C & rad2deg*theta_constr0(i), rad2deg*theta_drange(i), +C & rad2deg*difi,0.25d0*for_thet_constr(i)*difi**4, +C & gloc(itheta+nphi-2,icg) +C endif + enddo return end #endif diff --git a/source/cluster/wham/src-M/gnmr1.f b/source/cluster/wham/src-M/gnmr1.f index 905e746..2357e6d 100644 --- a/source/cluster/wham/src-M/gnmr1.f +++ b/source/cluster/wham/src-M/gnmr1.f @@ -41,3 +41,34 @@ c------------------------------------------------------------------------------- return end c--------------------------------------------------------------------------------- +c--------------------------------------------------------------------------------- + double precision function rlornmr1(y,ymin,ymax,sigma) + implicit none + double precision y,ymin,ymax,sigma + double precision wykl /4.0d0/ + if (y.lt.ymin) then + rlornmr1=(ymin-y)**wykl/((ymin-y)**wykl+sigma**wykl) + else if (y.gt.ymax) then + rlornmr1=(y-ymax)**wykl/((y-ymax)**wykl+sigma**wykl) + else + rlornmr1=0.0d0 + endif + return + end +c------------------------------------------------------------------------------ + double precision function rlornmr1prim(y,ymin,ymax,sigma) + implicit none + double precision y,ymin,ymax,sigma + double precision wykl /4.0d0/ + if (y.lt.ymin) then + rlornmr1prim=-(ymin-y)**(wykl-1)*sigma**wykl*wykl/ + & ((ymin-y)**wykl+sigma**wykl)**2 + else if (y.gt.ymax) then + rlornmr1prim=(y-ymax)**(wykl-1)*sigma**wykl*wykl/ + & ((y-ymax)**wykl+sigma**wykl)**2 + else + rlornmr1prim=0.0d0 + endif + return + end + diff --git a/source/cluster/wham/src-M/include_unres/COMMON.TORCNSTR b/source/cluster/wham/src-M/include_unres/COMMON.TORCNSTR index 5c26487..845c44d 100644 --- a/source/cluster/wham/src-M/include_unres/COMMON.TORCNSTR +++ b/source/cluster/wham/src-M/include_unres/COMMON.TORCNSTR @@ -1,6 +1,14 @@ - integer ndih_constr,idih_constr(maxdih_constr) + integer ndih_constr,idih_constr(maxdih_constr),ntheta_constr, + & itheta_constr(maxdih_constr) integer ndih_nconstr,idih_nconstr(maxdih_constr) + integer idihconstr_start,idihconstr_end,ithetaconstr_start, + & ithetaconstr_end double precision phi0(maxdih_constr),drange(maxdih_constr), - & ftors(maxdih_constr) - common /torcnstr/ phi0,drange,ftors,ndih_constr,idih_constr, - & ndih_nconstr,idih_nconstr + & ftors(maxdih_constr),theta_constr0(maxdih_constr), + & theta_drange(maxdih_constr),for_thet_constr(maxdih_constr) + common /torcnstr/ phi0,drange,ftors,theta_constr0,theta_drange, + & for_thet_constr, + & ndih_constr,idih_constr, + & ndih_nconstr,idih_nconstr,idihconstr_start,idihconstr_end, + & ntheta_constr,itheta_constr,ithetaconstr_start, + & ithetaconstr_end diff --git a/source/cluster/wham/src-M/initialize_p.F b/source/cluster/wham/src-M/initialize_p.F index 4ebced6..f385f1a 100644 --- a/source/cluster/wham/src-M/initialize_p.F +++ b/source/cluster/wham/src-M/initialize_p.F @@ -256,14 +256,16 @@ c------------------------------------------------------------------------- & "EVDW SC-SC","EVDW2 SC-p","EES p-p","ECORR4 ","ECORR5 ", & "ECORR6 ","EELLO ","ETURN3 ","ETURN4 ","ETURN6 ", & "EBE bend","ESC SCloc","ETORS ","ETORSD ","EHPB","EVDWPP", - & "ESTR","ESCCOR","EVDW2_14","","EVDW_T"/ + & "EVDW2_14","ESTR","ESCCOR","EDIHC","EVDW_T","ELIPTRAN", + & "EAFM","ETHETC"/ data wname / & "WSC","WSCP","WELEC","WCORR","WCORR5","WCORR6","WEL_LOC", & "WTURN3","WTURN4","WTURN6","WANG","WSCLOC","WTOR","WTORD", - & "WHPB","WVDWPP","WBOND","WSCCOR","WSCP14","","WSC"/ - data nprint_ene /21/ + & "WHPB","WVDWPP","WSCP14","WBOND","WSCCOR","WDIHC","WSC", + & "WLIPTRAN","WAFM","WTHETC"/ + data nprint_ene /22/ data print_order /1,2,3,18,11,12,13,14,4,5,6,7,8,9,10,19, - & 16,15,17,20,21/ + & 16,15,17,20,21,24,22,23/ end c--------------------------------------------------------------------------- subroutine init_int_table diff --git a/source/cluster/wham/src-M/readrtns.F b/source/cluster/wham/src-M/readrtns.F index 8624355..aa85ca1 100644 --- a/source/cluster/wham/src-M/readrtns.F +++ b/source/cluster/wham/src-M/readrtns.F @@ -39,6 +39,8 @@ C call readi(controlcard,'CONSTR_DIST',constr_dist,0) write (iout,*) "with_dihed_constr ",with_dihed_constr, & " CONSTR_DIST",constr_dist + with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0 + write (iout,*) "with_theta_constr ",with_theta_constr call flush(iout) min_var=(index(controlcard,'MINVAR').gt.0) plot_tree=(index(controlcard,'PLOT_TREE').gt.0) @@ -306,6 +308,44 @@ C ftors is the force constant for torsional quartic constrains enddo endif ! endif ndif_constr.gt.0 endif ! with_dihed_constr + 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 + 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 +C 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 +C 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 nct=nres print *,'NNT=',NNT,' NCT=',NCT diff --git a/source/unres/src_MD-M/COMMON.CONTROL b/source/unres/src_MD-M/COMMON.CONTROL index 40346c7..4d05560 100644 --- a/source/unres/src_MD-M/COMMON.CONTROL +++ b/source/unres/src_MD-M/COMMON.CONTROL @@ -3,11 +3,12 @@ logical minim,refstr,pdbref,outpdb,outmol2,overlapsc,energy_dec, & sideadd,lsecondary,read_cart,unres_pdb, & vdisulf,searchsc,lmuca,dccart,extconf,out1file, - & gnorm_check,gradout,split_ene + & gnorm_check,gradout,split_ene,with_theta_constr common /cntrl/ modecalc,iscode,indpdb,indback,indphi,iranconf, & icheckgrad,minim,i2ndstr,refstr,pdbref,outpdb,outmol2,iprint, & overlapsc,energy_dec,sideadd,lsecondary,read_cart,unres_pdb & ,vdisulf,searchsc,lmuca,dccart,mucadyn,extconf,out1file, - & constr_dist,gnorm_check,gradout,split_ene,symetr + & constr_dist,gnorm_check,gradout,split_ene,with_theta_constr, + & symetr C... minim = .true. means DO minimization. C... energy_dec = .true. means print energy decomposition matrix diff --git a/source/unres/src_MD-M/COMMON.TORCNSTR b/source/unres/src_MD-M/COMMON.TORCNSTR index 4d77fc6..845c44d 100644 --- a/source/unres/src_MD-M/COMMON.TORCNSTR +++ b/source/unres/src_MD-M/COMMON.TORCNSTR @@ -6,5 +6,9 @@ double precision phi0(maxdih_constr),drange(maxdih_constr), & ftors(maxdih_constr),theta_constr0(maxdih_constr), & theta_drange(maxdih_constr),for_thet_constr(maxdih_constr) - common /torcnstr/ phi0,drange,ftors,ndih_constr,idih_constr, - & ndih_nconstr,idih_nconstr,idihconstr_start,idihconstr_end + common /torcnstr/ phi0,drange,ftors,theta_constr0,theta_drange, + & for_thet_constr, + & ndih_constr,idih_constr, + & ndih_nconstr,idih_nconstr,idihconstr_start,idihconstr_end, + & ntheta_constr,itheta_constr,ithetaconstr_start, + & ithetaconstr_end diff --git a/source/unres/src_MD-M/DIMENSIONS b/source/unres/src_MD-M/DIMENSIONS index c8a95df..bda06b4 100644 --- a/source/unres/src_MD-M/DIMENSIONS +++ b/source/unres/src_MD-M/DIMENSIONS @@ -95,7 +95,7 @@ C Max. number of conformations in the pool parameter (max_pool=10) C Number of energy components integer n_ene,n_ene2 - parameter (n_ene=21,n_ene2=2*n_ene) + parameter (n_ene=25,n_ene2=2*n_ene) C Number of threads in deformation integer max_thread,max_thread2 parameter (max_thread=4,max_thread2=2*max_thread) diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F index ec0aef9..9414f1c 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -193,6 +193,7 @@ C call ebend(ebe,ethetacnstr) else ebe=0 + ethetacnstr=0 endif c print *,"Processor",myrank," computed UB" C @@ -301,7 +302,9 @@ C energia(17)=estr energia(20)=Uconst+Uconst_back energia(21)=esccor -C energia(22)= +C energia(22)=eliptrans (the energy for lipid transfere implemented in lipid branch) +C energia(23)= ... (energy for AFM, steered molecular dynamics) + energia(24)=ethetacnstr c Here are the energies showed per procesor if the are more processors c per molecule then we sum it up in sum_energy subroutine c print *," Processor",myrank," calls SUM_ENERGY" @@ -393,20 +396,22 @@ cMS$ATTRIBUTES C :: proc_proc estr=energia(17) Uconst=energia(20) esccor=energia(21) + ethetacnstr=energia(24) + #ifdef SPLITELE etot=wsc*evdw+wscp*evdw2+welec*ees+wvdwpp*evdw1 & +wang*ebe+wtor*etors+wscloc*escloc & +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5 & +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3 & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d - & +wbond*estr+Uconst+wsccor*esccor + & +wbond*estr+Uconst+wsccor*esccor+ethetacnstr #else etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1) & +wang*ebe+wtor*etors+wscloc*escloc & +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5 & +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3 & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d - & +wbond*estr+Uconst+wsccor*esccor + & +wbond*estr+Uconst+wsccor*esccor+ethetacnstr #endif energia(0)=etot c detecting NaNQ @@ -971,6 +976,7 @@ C------------------------------------------------------------------------ estr=energia(17) Uconst=energia(20) esccor=energia(21) + ethetacnstr=energia(24) #ifdef SPLITELE write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,evdw1,wvdwpp, & estr,wbond,ebe,wang, @@ -4405,7 +4411,7 @@ c end #ifdef CRYST_THETA C-------------------------------------------------------------------------- - subroutine ebend(etheta) + subroutine ebend(etheta,ethetacnstr) C C Evaluate the virtual-bond-angle energy given the virtual-bond dihedral C angles gamma and its derivatives in consecutive thetas and gammas. @@ -4422,6 +4428,7 @@ C include 'COMMON.NAMES' include 'COMMON.FFIELD' include 'COMMON.CONTROL' + include 'COMMON.TORCNSTR' common /calcthet/ term1,term2,termm,diffak,ratak, & ak,aktc,termpre,termexp,sigc,sig0i,time11,time12,sigcsq, & delthe0,sig0inv,sigtc,sigsqtc,delthec,it @@ -4525,18 +4532,26 @@ C Derivatives of the "mean" values in gamma1 and gamma2. & E_theta,E_tc) endif etheta=etheta+ethetai + if (energy_dec) write (iout,'(a6,i5,0pf7.3)') + & 'ebend',i,ethetai + if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*E_tc*dthetg1 + if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*E_tc*dthetg2 + gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett) + enddo + ethetacnstr=0.0d0 +C print *,ithetaconstr_start,ithetaconstr_end,"TU" do i=ithetaconstr_start,ithetaconstr_end itheta=itheta_constr(i) - thetiii=theta(itori) + thetiii=theta(itheta) difi=pinorm(thetiii-theta_constr0(i)) if (difi.gt.theta_drange(i)) then difi=difi-theta_drange(i) - ethetcnstr=ethetcnstr+0.25d0*for_thet_constr(i)*difi**4 + ethetacnstr=ethetcnstr+0.25d0*for_thet_constr(i)*difi**4 gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) & +for_thet_constr(i)*difi**3 else if (difi.lt.-drange(i)) then difi=difi+drange(i) - ethetcnstr=ethetcnstr+0.25d0*for_thet_constr(i)*difi**4 + ethetacnstr=ethetcnstr+0.25d0*for_thet_constr(i)*difi**4 gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) & +for_thet_constr(i)*difi**3 else @@ -4550,12 +4565,7 @@ C Derivatives of the "mean" values in gamma1 and gamma2. & gloc(itheta+nphi-2,icg) endif enddo - if (energy_dec) write (iout,'(a6,i5,0pf7.3)') - & 'ebend',i,ethetai - if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*E_tc*dthetg1 - if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*E_tc*dthetg2 - gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett) - enddo + C Ufff.... We've done all this!!! return end @@ -4669,7 +4679,7 @@ C "Thank you" to MAPLE (probably spared one day of hand-differentiation). end #else C-------------------------------------------------------------------------- - subroutine ebend(etheta) + subroutine ebend(etheta,ethetacnstr) C C Evaluate the virtual-bond-angle energy given the virtual-bond dihedral C angles gamma and its derivatives in consecutive thetas and gammas. @@ -4688,6 +4698,7 @@ C include 'COMMON.NAMES' include 'COMMON.FFIELD' include 'COMMON.CONTROL' + include 'COMMON.TORCNSTR' double precision coskt(mmaxtheterm),sinkt(mmaxtheterm), & cosph1(maxsingle),sinph1(maxsingle),cosph2(maxsingle), & sinph2(maxsingle),cosph1ph2(maxdouble,maxdouble), @@ -4862,19 +4873,33 @@ C print *,"tu",cosph1(k),sinph1(k),cosph2(k),sinph2(k) enddo enddo 10 continue -C now we have the theta_constrains +c lprn1=.true. +C print *,ethetai + if (lprn1) + & write (iout,'(i2,3f8.1,9h ethetai ,f10.5)') + & i,theta(i)*rad2deg,phii*rad2deg, + & phii1*rad2deg,ethetai +c lprn1=.false. + etheta=etheta+ethetai + if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*dephii + if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*dephii1 + gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wang*dethetai + enddo +C now constrains + ethetacnstr=0.0d0 +C print *,ithetaconstr_start,ithetaconstr_end,"TU" do i=ithetaconstr_start,ithetaconstr_end itheta=itheta_constr(i) - thetiii=theta(itori) + thetiii=theta(itheta) difi=pinorm(thetiii-theta_constr0(i)) if (difi.gt.theta_drange(i)) then difi=difi-theta_drange(i) - ethetcnstr=ethetcnstr+0.25d0*for_thet_constr(i)*difi**4 + ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4 gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) & +for_thet_constr(i)*difi**3 else if (difi.lt.-drange(i)) then difi=difi+drange(i) - ethetcnstr=ethetcnstr+0.25d0*for_thet_constr(i)*difi**4 + ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4 gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) & +for_thet_constr(i)*difi**3 else @@ -4889,18 +4914,6 @@ C now we have the theta_constrains endif enddo -c lprn1=.true. -C print *,ethetai - if (lprn1) - & write (iout,'(i2,3f8.1,9h ethetai ,f10.5)') - & i,theta(i)*rad2deg,phii*rad2deg, - & phii1*rad2deg,ethetai -c lprn1=.false. - etheta=etheta+ethetai - if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*dephii - if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*dephii1 - gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wang*dethetai - enddo return end #endif diff --git a/source/unres/src_MD-M/initialize_p.F b/source/unres/src_MD-M/initialize_p.F index 6c5f4eb..bdc3269 100644 --- a/source/unres/src_MD-M/initialize_p.F +++ b/source/unres/src_MD-M/initialize_p.F @@ -670,7 +670,7 @@ c nlen=nres-nnt+1 & ' ivec_start',ivec_start,' ivec_end',ivec_end, & ' iset_start',iset_start,' iset_end',iset_end, & ' idihconstr_start',idihconstr_start,' idihconstr_end', - & idihconstr_end + & idihconstr_end, & ' ithetaconstr_start',ithetaconstr_start,' ithetaconstr_end', & ithetaconstr_end diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F index 286cc57..e10b3ba 100644 --- a/source/unres/src_MD-M/readrtns_CSA.F +++ b/source/unres/src_MD-M/readrtns_CSA.F @@ -856,11 +856,12 @@ C & write (iout,*) 'FTORS',ftors enddo endif C first setting the theta boundaries to 0 to pi -C this mean that there is no energy penalty for any angle occuring - do i=1,nres - thetabound(1,i)=0 - thetabound(2,i)=pi - enddo +C this mean that there is no energy penalty for any angle occuring this can be applied +C for generate random conformation but is not implemented in this way +C do i=1,nres +C thetabound(1,i)=0 +C thetabound(2,i)=pi +C enddo C begin reading theta constrains this is quartic constrains allowing to C have smooth second derivative if (with_theta_constr) then @@ -887,17 +888,17 @@ C E=k*x**4 & for_thet_constr(i) enddo endif - do i=1,nthet_constr + 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 - do i=1,ntheta_constr - ii = itheta_constr(i) - thetabound(1,ii) = phi0(i)-drange(i) - thetabound(2,ii) = phi0(i)+drange(i) - enddo +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 C diff --git a/source/wham/src-M/COMMON.CONTROL b/source/wham/src-M/COMMON.CONTROL index 6c87389..122ff8d 100644 --- a/source/wham/src-M/COMMON.CONTROL +++ b/source/wham/src-M/COMMON.CONTROL @@ -2,9 +2,10 @@ & ensembles,constr_dist,symetr logical refstr,pdbref,punch_dist,print_rms,caonly,verbose, & merge_helices,bxfile,cxfile,histfile,entfile,zscfile, - & rmsrgymap,with_dihed_constr,check_conf,histout + & rmsrgymap,with_dihed_constr,check_conf,histout,with_theta_constr common /cntrl/ iscode,indpdb,refstr,pdbref,outpdb,outmol2, & punch_dist,print_rms,caonly,verbose,icomparfunc,pdbint, & merge_helices,bxfile,cxfile,histfile,entfile,zscfile,rmsrgymap, & ensembles,with_dihed_constr,constr_dist,check_conf,histout, + & with_theta_constr, &symetr diff --git a/source/wham/src-M/DIMENSIONS.ZSCOPT b/source/wham/src-M/DIMENSIONS.ZSCOPT index cfc74a4..dbeffcd 100644 --- a/source/wham/src-M/DIMENSIONS.ZSCOPT +++ b/source/wham/src-M/DIMENSIONS.ZSCOPT @@ -3,7 +3,7 @@ c Maximum number of structures in the database, energy components, proteins, c and structural classes c#ifdef JUBL - parameter (maxstr=200000,max_ene=21,maxprot=7,maxclass=5000) + parameter (maxstr=200000,max_ene=24,maxprot=7,maxclass=5000) parameter (maxclass1=10) c Maximum number of structures to be dealt with by one processor parameter (maxstr_proc=10000) diff --git a/source/wham/src-M/energy_p_new.F b/source/wham/src-M/energy_p_new.F index 33e670a..44e5cba 100644 --- a/source/wham/src-M/energy_p_new.F +++ b/source/wham/src-M/energy_p_new.F @@ -67,7 +67,7 @@ cd print *,'EHPB exitted succesfully.' C C Calculate the virtual-bond-angle energy. C - call ebend(ebe) + call ebend(ebe,ethetacnstr) cd print *,'Bend energy finished.' C C Calculate the SC local energy. @@ -111,7 +111,7 @@ c write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t & +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4 & +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6 & +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d - & +wbond*estr+wsccor*fact(1)*esccor + & +wbond*estr+wsccor*fact(1)*esccor+ethetacnstr #else etot=wsc*(evdw+fact(6)*evdw_t)+wscp*evdw2 & +welec*fact(1)*(ees+evdw1) @@ -120,7 +120,7 @@ c write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t & +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4 & +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6 & +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d - & +wbond*estr+wsccor*fact(1)*esccor + & +wbond*estr+wsccor*fact(1)*esccor+ethetacnstr #endif energia(0)=etot energia(1)=evdw @@ -154,6 +154,7 @@ c write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t energia(19)=esccor energia(20)=edihcnstr energia(21)=evdw_t + energia(24)=ethetacnstr c detecting NaNQ #ifdef ISNAN #ifdef AIX @@ -270,6 +271,7 @@ C------------------------------------------------------------------------ esccor=energia(19) edihcnstr=energia(20) estr=energia(18) + ethetacnstr=energia(24) #ifdef SPLITELE write (iout,10) evdw,wsc,evdw2,wscp,ees,welec*fact(1),evdw1, & wvdwpp, @@ -278,7 +280,7 @@ C------------------------------------------------------------------------ & ecorr,wcorr*fact(3),ecorr5,wcorr5*fact(4),ecorr6,wcorr6*fact(5), & eel_loc,wel_loc*fact(2),eello_turn3,wturn3*fact(2), & eello_turn4,wturn4*fact(3),eello_turn6,wturn6*fact(5), - & esccor,wsccor*fact(1),edihcnstr,ebr*nss,etot + & esccor,wsccor*fact(1),edihcnstr,ethetacnstr,ebr*nss,etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ @@ -300,6 +302,7 @@ C------------------------------------------------------------------------ & 'ETURN6=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 6th order)'/ & 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/ & 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/ + & 'ETHETC= ',1pE16.6,' (valence angle constraints)'/ & 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ & 'ETOT= ',1pE16.6,' (total)') #else @@ -309,7 +312,7 @@ C------------------------------------------------------------------------ & ecorr6,wcorr6*fact(5),eel_loc,wel_loc*fact(2), & eello_turn3,wturn3*fact(2),eello_turn4,wturn4*fact(3), & eello_turn6,wturn6*fact(5),esccor*fact(1),wsccor, - & edihcnstr,ebr*nss,etot + & edihcnstr,ethetacnstr,ebr*nss,etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ @@ -330,6 +333,7 @@ C------------------------------------------------------------------------ & 'ETURN6=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 6th order)'/ & 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/ & 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/ + & 'ETHETC= ',1pE16.6,' (valence angle constraints)'/ & 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ & 'ETOT= ',1pE16.6,' (total)') #endif @@ -3265,7 +3269,7 @@ c & AKSC(j,iti),abond0(j,iti),u(j),j=1,nbi) end #ifdef CRYST_THETA C-------------------------------------------------------------------------- - subroutine ebend(etheta) + subroutine ebend(etheta,ethetacnstr) C C Evaluate the virtual-bond-angle energy given the virtual-bond dihedral C angles gamma and its derivatives in consecutive thetas and gammas. @@ -3282,6 +3286,7 @@ C include 'COMMON.IOUNITS' include 'COMMON.NAMES' include 'COMMON.FFIELD' + include 'COMMON.TORCNSTR' common /calcthet/ term1,term2,termm,diffak,ratak, & ak,aktc,termpre,termexp,sigc,sig0i,time11,time12,sigcsq, & delthe0,sig0inv,sigtc,sigsqtc,delthec,it @@ -3400,6 +3405,33 @@ c & rad2deg*phii,rad2deg*phii1,ethetai gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett) c 1215 continue enddo + ethetacnstr=0.0d0 +C print *,ithetaconstr_start,ithetaconstr_end,"TU" + do i=1,ntheta_constr + itheta=itheta_constr(i) + thetiii=theta(itheta) + difi=pinorm(thetiii-theta_constr0(i)) + if (difi.gt.theta_drange(i)) then + difi=difi-theta_drange(i) + ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4 + gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) + & +for_thet_constr(i)*difi**3 + else if (difi.lt.-drange(i)) then + difi=difi+drange(i) + ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4 + gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) + & +for_thet_constr(i)*difi**3 + else + difi=0.0 + endif +C if (energy_dec) then +C write (iout,'(a6,2i5,4f8.3,2e14.5)') "ethetc", +C & i,itheta,rad2deg*thetiii, +C & rad2deg*theta_constr0(i), rad2deg*theta_drange(i), +C & rad2deg*difi,0.25d0*for_thet_constr(i)*difi**4, +C & gloc(itheta+nphi-2,icg) + endif + enddo C Ufff.... We've done all this!!! return end @@ -3513,7 +3545,7 @@ C "Thank you" to MAPLE (probably spared one day of hand-differentiation). end #else C-------------------------------------------------------------------------- - subroutine ebend(etheta) + subroutine ebend(etheta,ethetacnstr) C C Evaluate the virtual-bond-angle energy given the virtual-bond dihedral C angles gamma and its derivatives in consecutive thetas and gammas. @@ -3533,6 +3565,7 @@ C include 'COMMON.NAMES' include 'COMMON.FFIELD' include 'COMMON.CONTROL' + include 'COMMON.TORCNSTR' double precision coskt(mmaxtheterm),sinkt(mmaxtheterm), & cosph1(maxsingle),sinph1(maxsingle),cosph2(maxsingle), & sinph2(maxsingle),cosph1ph2(maxdouble,maxdouble), @@ -3712,6 +3745,34 @@ c call flush(iout) c gloc(nphi+i-2,icg)=wang*dethetai gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wang*dethetai enddo +C now constrains + ethetacnstr=0.0d0 +C print *,ithetaconstr_start,ithetaconstr_end,"TU" + do i=1,ntheta_constr + itheta=itheta_constr(i) + thetiii=theta(itheta) + difi=pinorm(thetiii-theta_constr0(i)) + if (difi.gt.theta_drange(i)) then + difi=difi-theta_drange(i) + ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4 + gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) + & +for_thet_constr(i)*difi**3 + else if (difi.lt.-drange(i)) then + difi=difi+drange(i) + ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4 + gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) + & +for_thet_constr(i)*difi**3 + else + difi=0.0 + endif +C if (energy_dec) then +C write (iout,'(a6,2i5,4f8.3,2e14.5)') "ethetc", +C & i,itheta,rad2deg*thetiii, +C & rad2deg*theta_constr0(i), rad2deg*theta_drange(i), +C & rad2deg*difi,0.25d0*for_thet_constr(i)*difi**4, +C & gloc(itheta+nphi-2,icg) +C endif + enddo return end #endif diff --git a/source/wham/src-M/include_unres/COMMON.TORCNSTR b/source/wham/src-M/include_unres/COMMON.TORCNSTR index bef2264..845c44d 100644 --- a/source/wham/src-M/include_unres/COMMON.TORCNSTR +++ b/source/wham/src-M/include_unres/COMMON.TORCNSTR @@ -1,6 +1,14 @@ - integer ndih_constr,idih_constr(maxdih_constr) + integer ndih_constr,idih_constr(maxdih_constr),ntheta_constr, + & itheta_constr(maxdih_constr) integer ndih_nconstr,idih_nconstr(maxdih_constr) + integer idihconstr_start,idihconstr_end,ithetaconstr_start, + & ithetaconstr_end double precision phi0(maxdih_constr),drange(maxdih_constr), - & ftors(maxdih_constr) - common /torcnstr/ phi0,drange,ftors,ndih_constr,idih_constr, - & ndih_nconstr,idih_nconstr + & ftors(maxdih_constr),theta_constr0(maxdih_constr), + & theta_drange(maxdih_constr),for_thet_constr(maxdih_constr) + common /torcnstr/ phi0,drange,ftors,theta_constr0,theta_drange, + & for_thet_constr, + & ndih_constr,idih_constr, + & ndih_nconstr,idih_nconstr,idihconstr_start,idihconstr_end, + & ntheta_constr,itheta_constr,ithetaconstr_start, + & ithetaconstr_end diff --git a/source/wham/src-M/initialize_p.F b/source/wham/src-M/initialize_p.F index c890e09..aea09dd 100644 --- a/source/wham/src-M/initialize_p.F +++ b/source/wham/src-M/initialize_p.F @@ -262,17 +262,19 @@ c------------------------------------------------------------------------- & "EVDW SC-SC","EVDW2 SC-p","EES p-p","ECORR4 ","ECORR5 ", & "ECORR6 ","EELLO ","ETURN3 ","ETURN4 ","ETURN6 ", & "EBE bend","ESC SCloc","ETORS ","ETORSD ","EHPB","EVDWPP", - & "EVDW2_14","ESTR","ESCCOR","EDIHC","EVDW_T"/ + & "EVDW2_14","ESTR","ESCCOR","EDIHC","EVDW_T","ELIPTRAN", + & "EAFM","ETHETC"/ data wname / & "WSC","WSCP","WELEC","WCORR","WCORR5","WCORR6","WEL_LOC", & "WTURN3","WTURN4","WTURN6","WANG","WSCLOC","WTOR","WTORD", - & "WHPB","WVDWPP","WSCP14","WBOND","WSCCOR","WDIHC","WSC"/ + & "WHPB","WVDWPP","WSCP14","WBOND","WSCCOR","WDIHC","WSC", + & "WLIPTRAN","WAFM","WTHETC"/ data ww0 /1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0, & 1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,0.4d0,1.0d0,1.0d0, - & 0.0d0,0.0/ - data nprint_ene /21/ + & 0.0d0,0.0,0.0d0,0.0d0,0.0d0/ + data nprint_ene /22/ data print_order /1,2,3,18,11,12,13,14,4,5,6,7,8,9,10,19, - & 16,15,17,20,21/ + & 16,15,17,20,21,24,22,23/ end c--------------------------------------------------------------------------- subroutine init_int_table diff --git a/source/wham/src-M/molread_zs.F b/source/wham/src-M/molread_zs.F index 061ae05..2b9fe8c 100644 --- a/source/wham/src-M/molread_zs.F +++ b/source/wham/src-M/molread_zs.F @@ -94,7 +94,43 @@ C write (iout,*) 'FTORS',ftors endif 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 + 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 +C 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 +C 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 nct=nres if (itype(1).eq.ntyp1) nnt=2