adding some param and control
[unres.git] / source / unres / src_MD-M / readrtns_CSA.F
index 91d7c07..fa4c75a 100644 (file)
@@ -81,6 +81,7 @@ C
       include 'COMMON.INTERACT'
       include 'COMMON.SETUP'
       include 'COMMON.SPLITELE'
+      include 'COMMON.SHIELD'
       COMMON /MACHSW/ KDIAG,ICORFL,IXDR
       character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
       character*80 ucase
@@ -103,7 +104,7 @@ C constrains on theta angles WITH_THETA_CONSTR is the keyword
       with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
       write (iout,*) "with_theta_constr ",with_theta_constr
       call readi(controlcard,'SYM',symetr,1)
-      call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
+      call reada(controlcard,'TIMLIM',timlim,2800.0D0) ! default 16 hours
       unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
       call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
       call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
@@ -144,7 +145,23 @@ C constrains on theta angles WITH_THETA_CONSTR is the keyword
       AFMlog=(index(controlcard,'AFM'))
       selfguide=(index(controlcard,'SELFGUIDE'))
       print *,'AFMlog',AFMlog,selfguide,"KUPA"
+      call readi(controlcard,'TUBEMOD',tubelog,0)
+      write (iout,*) TUBElog,"TUBEMODE"
+      call readi(controlcard,'GENCONSTR',genconstr,0)
+C      write (iout,*) TUBElog,"TUBEMODE"
       call readi(controlcard,'IPRINT',iprint,0)
+C SHIELD keyword sets if the shielding effect of side-chains is used
+C 0 denots no shielding is used all peptide are equally despite the 
+C solvent accesible area
+C 1 the newly introduced function
+C 2 reseved for further possible developement
+      call readi(controlcard,'SHIELD',shield_mode,0)
+C      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
+        write(iout,*) "shield_mode",shield_mode
+C      endif
+      call readi(controlcard,'TORMODE',tor_mode,0)
+C      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
+        write(iout,*) "torsional and valence angle mode",tor_mode
       call readi(controlcard,'MAXGEN',maxgen,10000)
       call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
       call readi(controlcard,"KDIAG",kdiag,0)
@@ -166,6 +183,7 @@ C constrains on theta angles WITH_THETA_CONSTR is the keyword
         else
           icheckgrad=3
         endif
+        call reada(controlcard,'DELTA',aincr,1.0d-7)
       elseif (index(controlcard,'THREAD').gt.0) then
         modecalc=2
         call readi(controlcard,'THREAD',nthread,0)
@@ -237,7 +255,7 @@ c Cutoff range for interactions
        bordliptop=(boxzsize+lipthick)/2.0
        bordlipbot=bordliptop-lipthick
 C      endif
-      if ((bordliptop.gt.boxzsize).or.(borlipbot.lt.0.0)) 
+      if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) 
      & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
       buflipbot=bordlipbot+lipbufthick
       bufliptop=bordliptop-lipbufthick
@@ -248,8 +266,34 @@ C      endif
       write(iout,*) "bordlipbot=",bordlipbot
       write(iout,*) "bufliptop=",bufliptop
       write(iout,*) "buflipbot=",buflipbot
-
-
+      write (iout,*) "SHIELD MODE",shield_mode
+      if (TUBElog.gt.0) then
+       call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
+       call reada(controlcard,"YTUBE",tubecenter(2),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
+      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 (me.eq.king .or. .not.out1file ) 
      &  write (iout,*) "DISTCHAINMAX",distchainmax
       
@@ -571,6 +615,7 @@ C
       include 'COMMON.BOUNDS'
       include 'COMMON.MD'
       include 'COMMON.SETUP'
+      include 'COMMON.SHIELD'
       character*4 sequence(maxres)
       integer rescode
       double precision x(maxvar)
@@ -612,7 +657,9 @@ C Read weights of the subsequent energy terms.
        call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
        call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
        call reada(weightcard,'TEMP0',temp0,300.0d0)
+       call reada(weightcard,'WSHIELD',wshield,1.0d0)
        call reada(weightcard,'WLT',wliptran,0.0D0)
+       call reada(weightcard,'WTUBE',wtube,1.0D0)
        if (index(weightcard,'SOFT').gt.0) ipot=6
 C 12/1/95 Added weight for the multi-body term WCORR
        call reada(weightcard,'WCORRH',wcorr,1.0D0)
@@ -721,10 +768,16 @@ C 12/1/95 Added weight for the multi-body term WCORR
       call reada(weightcard,"BTRISS",btriss,0.021D0)
       call reada(weightcard,"CTRISS",ctriss,1.001D0)
       call reada(weightcard,"DTRISS",dtriss,1.001D0)
+      call reada(weightcard,"LIPSCALE",lipscale,1.3D0)
       write (iout,*) "ATRISS=", atriss
       write (iout,*) "BTRISS=", btriss
       write (iout,*) "CTRISS=", ctriss
       write (iout,*) "DTRISS=", dtriss
+       if (shield_mode.gt.0) then
+        lipscale=0.0d0
+        write (iout,*) "liscale not used in shield mode"
+       endif
+      write (iout,*) "lipid scaling factor", lipscale
       dyn_ss=(index(weightcard,'DYN_SS').gt.0)
       do i=1,maxres
         dyn_ss_mask(i)=.false.
@@ -738,12 +791,12 @@ C 12/1/95 Added weight for the multi-body term WCORR
       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
+        akcm=akcm/wsc
+        akth=akth/wsc
+        akct=akct/wsc
+        v1ss=v1ss/wsc
+        v2ss=v2ss/wsc
+        v3ss=v3ss/wsc
       else
         if (wstrain.ne.0.0) then
          ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
@@ -751,7 +804,7 @@ C 12/1/95 Added weight for the multi-body term WCORR
           ss_depth=0.0
         endif
       endif
-
+      write (iout,*) "wshield,", wshield
       if(me.eq.king.or..not.out1file) then
        write (iout,*) "Parameters of the SS-bond potential:"
        write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
@@ -1009,7 +1062,7 @@ c----------------------
           call MPI_Finalize(MPI_COMM_WORLD,IERROR)
           stop 'Error reading reference structure'
 #endif
-   39     call chainbuild
+   39     call chainbuild_extconf
           call setup_var
 czscore          call geom_to_var(nvar,coord_exp_zs(1,1))
           nstart_sup=nnt
@@ -1082,6 +1135,7 @@ C initial geometry.
             return
           else
             call read_angles(inp,*36)
+            call chainbuild_extconf
           endif
           goto 37
    36     write (iout,'(a)') 'Error reading angle file.'
@@ -1106,6 +1160,7 @@ C initial geometry.
           omeg(i)=-120d0*deg2rad
           if (itype(i).le.0) omeg(i)=-omeg(i)
          enddo
+        call chainbuild_extconf
         else
           if(me.eq.king.or..not.out1file)
      &     write (iout,'(a)') 'Random-generated initial geometry.'
@@ -1474,6 +1529,31 @@ cd      print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
 cd      write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
 cd     & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
 cd     & ' nsup',nsup
+      if (constr_dist.eq.11) then
+             do i=nstart_sup,nstart_sup+nsup-1
+              do j=i+2,nstart_sup+nsup-1
+                 distance=dist(i,j)
+                 if (distance.le.15.0) then
+                 nhpb=nhpb+1
+                 ihpb(nhpb)=i+nstart_seq-nstart_sup
+                 jhpb(nhpb)=j+nstart_seq-nstart_sup
+                 forcon(nhpb)=sqrt(0.04*distance)
+                 fordepth(nhpb)=sqrt(40.0/distance)
+                 dhpb(nhpb)=distance-0.1d0
+                 dhpb1(nhpb)=distance+0.1d0
+                
+#ifdef MPI
+          if (.not.out1file .or. me.eq.king)
+     &    write (iout,'(a,3i5,f8.2,f10.2)') "+dist.constr ",
+     &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+#else
+          write (iout,'(a,3i5,f8.2,f10.2)') "+dist.constr ",
+     &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+#endif
+            endif
+             enddo
+           enddo
+      else
       do i=nstart_sup,nstart_sup+nsup-1
 cd      write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
 cd     &    ' seq_pdb', restyp(itype_pdb(i))
@@ -1484,7 +1564,8 @@ cd     &    ' seq_pdb', restyp(itype_pdb(i))
           forcon(nhpb)=weidis
           dhpb(nhpb)=dist(i,j)
         enddo
-      enddo 
+      enddo
+      endif 
 cd      write (iout,'(a)') 'Distance constraints:' 
 cd      do i=nss+1,nhpb
 cd        ii=ihpb(i)
@@ -1973,12 +2054,18 @@ C Get parameter filenames and open the parameter files.
       open (itorp,file=torname,status='old',readonly,shared)
       call getenv_loc('TORDPAR',tordname)
       open (itordp,file=tordname,status='old',readonly,shared)
+      call getenv_loc('TORKCC',torkccname)
+      open (itorkcc,file=torkccname,status='old',readonly,shared)
+      call getenv_loc('THETKCC',thetkccname)
+      open (ithetkcc,file=thetkccname,status='old',readonly,shared)
       call getenv_loc('FOURIER',fouriername)
       open (ifourier,file=fouriername,status='old',readonly,shared)
       call getenv_loc('ELEPAR',elename)
       open (ielep,file=elename,status='old',readonly,shared)
       call getenv_loc('SIDEPAR',sidename)
       open (isidep,file=sidename,status='old',readonly,shared)
+      call getenv_loc('LIPTRANPAR',liptranname)
+      open (iliptranpar,file=liptranname,status='old',readonly,shared)
 #elif (defined CRAY) || (defined AIX)
       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
      &  action='read')
@@ -2001,6 +2088,10 @@ c      print *,"Processor",myrank," opened file IROTAM"
 c      print *,"Processor",myrank," opened file ITORP" 
       call getenv_loc('TORDPAR',tordname)
       open (itordp,file=tordname,status='old',action='read')
+      call getenv_loc('TORKCC',torkccname)
+      open (itorkcc,file=torkccname,status='old',action='read')
+      call getenv_loc('THETKCC',thetkccname)
+      open (ithetkcc,file=thetkccname,status='old',action='read')
 c      print *,"Processor",myrank," opened file ITORDP" 
       call getenv_loc('SCCORPAR',sccorname)
       open (isccor,file=sccorname,status='old',action='read')
@@ -2013,6 +2104,8 @@ c      print *,"Processor",myrank," opened file IFOURIER"
 c      print *,"Processor",myrank," opened file IELEP" 
       call getenv_loc('SIDEPAR',sidename)
       open (isidep,file=sidename,status='old',action='read')
+      call getenv_loc('LIPTRANPAR',liptranname)
+      open (iliptranpar,file=liptranname,status='old',action='read')
 c      print *,"Processor",myrank," opened file ISIDEP" 
 c      print *,"Processor",myrank," opened parameter files" 
 #elif (defined G77)
@@ -2030,6 +2123,10 @@ C Get parameter filenames and open the parameter files.
       open (itorp,file=torname,status='old')
       call getenv_loc('TORDPAR',tordname)
       open (itordp,file=tordname,status='old')
+      call getenv_loc('TORKCC',torkccname)
+      open (itorkcc,file=torkccname,status='old')
+      call getenv_loc('THETKCC',thetkccname)
+      open (ithetkcc,file=thetkccname,status='old')
       call getenv_loc('SCCORPAR',sccorname)
       open (isccor,file=sccorname,status='old')
       call getenv_loc('FOURIER',fouriername)
@@ -2038,6 +2135,8 @@ C Get parameter filenames and open the parameter files.
       open (ielep,file=elename,status='old')
       call getenv_loc('SIDEPAR',sidename)
       open (isidep,file=sidename,status='old')
+      call getenv_loc('LIPTRANPAR',liptranname)
+      open (iliptranpar,file=liptranname,status='old')
 #else
       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
      &  readonly)
@@ -2054,6 +2153,10 @@ C Get parameter filenames and open the parameter files.
       open (itorp,file=torname,status='old',readonly)
       call getenv_loc('TORDPAR',tordname)
       open (itordp,file=tordname,status='old',readonly)
+      call getenv_loc('TORKCC',torkccname)
+      open (itorkcc,file=torkccname,status='old',readonly)
+      call getenv_loc('THETKCC',thetkccname)
+      open (ithetkcc,file=thetkccname,status='old',readonly)
       call getenv_loc('SCCORPAR',sccorname)
       open (isccor,file=sccorname,status='old',readonly)
 #ifndef CRYST_THETA
@@ -2075,6 +2178,8 @@ C Get parameter filenames and open the parameter files.
       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
 #endif
 #endif
+      call getenv_loc('TUBEPAR',tubename)
+      open (itube,file=tubename,status='old',readonly)
 #ifndef OLDSCP
 C
 C 8/9/01 In the newest version SCp interaction constants are read from a file
@@ -2379,6 +2484,10 @@ c-------------------------------------------------------------------------------
       write (iout,*) "Calling read_dist_constr"
 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
 c      call flush(iout)
+      if ((genconstr.gt.0).and.(constr_dist.eq.11)) then
+      call gen_dist_constr
+      go to 1712
+      endif
       call card_concat(controlcard)
       call readi(controlcard,"NFRAG",nfrag_,0)
       call readi(controlcard,"NPAIR",npair_,0)
@@ -2504,6 +2613,7 @@ C          dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
 #endif
 
       enddo
+ 1712 continue
       call flush(iout)
       return
       end