include 'COMMON.FFIELD'
include 'COMMON.FREE'
include 'COMMON.INTERACT'
+ include "COMMON.SPLITELE"
character*320 controlcard,ucase
#ifdef MPL
include 'COMMON.INFO'
call readi(controlcard,'RESCALE',rescale_mode,2)
call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
write (iout,*) "DISTCHAINMAX",distchainmax
+C Reading the dimensions of box in x,y,z coordinates
+ call reada(controlcard,'BOXX',boxxsize,100.0d0)
+ call reada(controlcard,'BOXY',boxysize,100.0d0)
+ call reada(controlcard,'BOXZ',boxzsize,100.0d0)
+c Cutoff range for interactions
+ call reada(controlcard,"R_CUT",r_cut,15.0d0)
+ call reada(controlcard,"LAMBDA",rlamb,0.3d0)
call readi(controlcard,'PDBOUT',outpdb,0)
call readi(controlcard,'MOL2OUT',outmol2,0)
refstr=(index(controlcard,'REFSTR').gt.0)
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)
C
C Read weights of the subsequent energy terms.
call card_concat(weightcard)
+ write(iout,*) weightcard
call reada(weightcard,'WSC',wsc,1.0d0)
call reada(weightcard,'WLONG',wsc,wsc)
+ write(iout,*) "WLONG=",wsc
call reada(weightcard,'WSCP',wscp,1.0d0)
call reada(weightcard,'WELEC',welec,1.0D0)
call reada(weightcard,'WVDWPP',wvdwpp,welec)
call reada(weightcard,"V2SS",v2ss,7.61d0)
call reada(weightcard,"V3SS",v3ss,13.7d0)
call reada(weightcard,"EBR",ebr,-5.50D0)
+ write (iout,*) "atriss",atriss
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)
+ write(iout,*) "after"
write (iout,*) "ATRISS=", atriss
write (iout,*) "BTRISS=", btriss
write (iout,*) "CTRISS=", ctriss
& wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
& wturn4,wturn6,wsccor
10 format (/'Energy-term weights (unscaled):'//
- & 'WSCC= ',f10.6,' (SC-SC)'/
+ & 'WSC= ',f10.6,' (SC-SC)'/
& 'WSCP= ',f10.6,' (SC-p)'/
& 'WELEC= ',f10.6,' (p-p electr)'/
& 'WVDWPP= ',f10.6,' (p-p VDW)'/
read (inp,*) ndih_constr
if (ndih_constr.gt.0) then
- read (inp,*) ftors
- write (iout,*) 'FTORS',ftors
+C read (inp,*) ftors
+C write (iout,*) 'FTORS',ftors
C ftors is the force constant for torsional quartic constrains
- read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
+ read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
+ & i=1,ndih_constr)
write (iout,*)
& 'There are',ndih_constr,' constraints on phi angles.'
do i=1,ndih_constr
- write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
+ write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
+ & ftors(i)
enddo
do i=1,ndih_constr
phi0(i)=deg2rad*phi0(i)
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