X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fcluster%2Fwham%2Fsrc-M%2Freadrtns.F;h=21a2b7658ec0f659c85e02198f22913424684012;hb=e12c8db731e266c6e38b2f5883f3ae7c15bbfb98;hp=8624355d04078e176477e1d0bbb8f2c1c8bfb689;hpb=c65d296f3cdc480d3737eae335a510332b4f0dc8;p=unres.git

diff --git a/source/cluster/wham/src-M/readrtns.F b/source/cluster/wham/src-M/readrtns.F
index 8624355..21a2b76 100644
--- a/source/cluster/wham/src-M/readrtns.F
+++ b/source/cluster/wham/src-M/readrtns.F
@@ -15,12 +15,14 @@ C
       include 'COMMON.FFIELD'
       include 'COMMON.FREE'
       include 'COMMON.INTERACT'
+      include "COMMON.SPLITELE"
+      include 'COMMON.SHIELD'
       character*320 controlcard,ucase
 #ifdef MPL
       include 'COMMON.INFO'
 #endif
       integer i,i1,i2,it1,it2
-
+      double precision pi
       read (INP,'(a80)') titel
       call card_concat(controlcard)
 
@@ -28,6 +30,63 @@ C
       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 reada(controlcard,"LIPTHICK",lipthick,0.0d0)
+      call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
+      if (lipthick.gt.0.0d0) then
+       bordliptop=(boxzsize+lipthick)/2.0
+       bordlipbot=bordliptop-lipthick
+C      endif
+      if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
+     & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
+      buflipbot=bordlipbot+lipbufthick
+      bufliptop=bordliptop-lipbufthick
+      if ((lipbufthick*2.0d0).gt.lipthick)
+     &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
+      endif
+      write(iout,*) "bordliptop=",bordliptop
+      write(iout,*) "bordlipbot=",bordlipbot
+      write(iout,*) "bufliptop=",bufliptop
+      write(iout,*) "buflipbot=",buflipbot
+C Shielding mode
+      call readi(controlcard,'SHIELD',shield_mode,0)
+      write (iout,*) "SHIELD MODE",shield_mode
+      call readi(controlcard,'TUBEMOD',tubelog,0)
+      write (iout,*) TUBElog,"TUBEMODE"
+      if (shield_mode.gt.0) then
+      pi=3.141592d0
+C VSolvSphere the volume of solving sphere
+C      print *,pi,"pi"
+C rpp(1,1) is the energy r0 for peptide group contact and will be used for it 
+C there will be no distinction between proline peptide group and normal peptide
+C group in case of shielding parameters
+      VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
+      VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
+      write (iout,*) VSolvSphere,VSolvSphere_div
+C long axis of side chain 
+      do i=1,ntyp
+      long_r_sidechain(i)=vbldsc0(1,i)
+      short_r_sidechain(i)=sigma0(i)
+      enddo
+      buff_shield=1.0d0
+      endif
+      if (TUBElog.gt.0) then
+       call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
+       call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
+       call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
+       call reada(controlcard,"RTUBE",tubeR0,0.0d0)
+       call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
+       call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
+       call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
+       buftubebot=bordtubebot+tubebufthick
+       buftubetop=bordtubetop-tubebufthick
+      endif
       call readi(controlcard,'PDBOUT',outpdb,0)
       call readi(controlcard,'MOL2OUT',outmol2,0)
       refstr=(index(controlcard,'REFSTR').gt.0)
@@ -39,6 +98,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)
@@ -93,6 +154,7 @@ C
       include 'COMMON.CONTACTS'
       include 'COMMON.TIME1'
       include 'COMMON.TORCNSTR'
+      include 'COMMON.SHIELD'
 #ifdef MPL
       include 'COMMON.INFO'
 #endif
@@ -108,7 +170,9 @@ C Body
 C
 C Read weights of the subsequent energy terms.
       call card_concat(weightcard)
-      call reada(weightcard,'WSC',wsc,1.0d0)
+      write(iout,*) weightcard
+C      call reada(weightcard,'WSC',wsc,1.0d0)
+      write(iout,*) wsc
       call reada(weightcard,'WLONG',wsc,wsc)
       call reada(weightcard,'WSCP',wscp,1.0d0)
       call reada(weightcard,'WELEC',welec,1.0D0)
@@ -140,6 +204,12 @@ C Read weights of the subsequent energy terms.
       call reada(weightcard,"V2SS",v2ss,7.61d0)
       call reada(weightcard,"V3SS",v3ss,13.7d0)
       call reada(weightcard,"EBR",ebr,-5.50D0)
+      call reada(weightcard,'WSHIELD',wshield,1.0d0)
+      write(iout,*) 'WSHIELD',wshield
+      call reada(weightcard,'WTUBE',wtube,0.0d0)
+      write(iout,*) 'WTUBE',wtube
+      call reada(weightcard,"LIPSCALE",lipscale,1.3D0)
+       call reada(weightcard,'WLT',wliptran,0.0D0)
       call reada(weightcard,"ATRISS",atriss,0.301D0)
       call reada(weightcard,"BTRISS",btriss,0.021D0)
       call reada(weightcard,"CTRISS",ctriss,1.001D0)
@@ -148,6 +218,10 @@ C Read weights of the subsequent energy terms.
       write (iout,*) "BTRISS=", btriss
       write (iout,*) "CTRISS=", ctriss
       write (iout,*) "DTRISS=", dtriss
+       if (shield_mode.gt.0) then
+        lipscale=0.0d0
+        write (iout,*) "WARNING:liscale not used in shield mode"
+       endif
       dyn_ss=(index(weightcard,'DYN_SS').gt.0)
       do i=1,maxres
         dyn_ss_mask(i)=.false.
@@ -306,6 +380,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
@@ -706,6 +818,11 @@ C Get parameter filenames and open the parameter files.
       open (isidep1,file=sidepname,status="old")
       call getenv('SCCORPAR',sccorname)
       open (isccor,file=sccorname,status="old")
+      call getenv('LIPTRANPAR',liptranname)
+      open (iliptranpar,file=liptranname,status='old')
+      call getenv('TUBEPAR',tubename)
+      open (itube,file=tubename,status='old')
+
 #ifndef OLDSCP
 C
 C 8/9/01 In the newest version SCp interaction constants are read from a file