AFM+tube+kcc
[unres.git] / source / unres / src_MD-M / readrtns_CSA.F
index 286cc57..e40717f 100644 (file)
@@ -8,6 +8,7 @@
       include 'COMMON.CONTROL'
       include 'COMMON.SBRIDGE'
       include 'COMMON.IOUNITS'
+      include 'COMMON.SPLITELE'
       logical file_exist
 C Read force-field parameters except weights
       call parmread
@@ -79,6 +80,8 @@ C
       include 'COMMON.FFIELD'
       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
@@ -101,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)
@@ -139,7 +142,24 @@ C constrains on theta angles WITH_THETA_CONSTR is the keyword
       refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
       indpdb=index(controlcard,'PDBSTART')
       extconf=(index(controlcard,'EXTCONF').gt.0)
+      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,'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)
@@ -161,6 +181,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)
@@ -217,7 +238,60 @@ cfmc        modecalc=10
       i2ndstr=index(controlcard,'USE_SEC_PRED')
       gradout=index(controlcard,'GRADOUT').gt.0
       gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
+C  DISTCHAINMAX become obsolete for periodic boundry condition
       call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
+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
+      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
       
@@ -344,8 +418,8 @@ C
       ntime_split0=ntime_split
       call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
       ntime_split0=ntime_split
-      call reada(controlcard,"R_CUT",r_cut,2.0d0)
-      call reada(controlcard,"LAMBDA",rlamb,0.3d0)
+c      call reada(controlcard,"R_CUT",r_cut,2.0d0)
+c      call reada(controlcard,"LAMBDA",rlamb,0.3d0)
       rest = index(controlcard,"REST").gt.0
       tbf = index(controlcard,"TBF").gt.0
       usampl = index(controlcard,"USAMPL").gt.0
@@ -539,6 +613,7 @@ C
       include 'COMMON.BOUNDS'
       include 'COMMON.MD'
       include 'COMMON.SETUP'
+      include 'COMMON.SHIELD'
       character*4 sequence(maxres)
       integer rescode
       double precision x(maxvar)
@@ -580,6 +655,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)
@@ -718,7 +796,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,
@@ -856,11 +934,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 +966,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
@@ -975,7 +1054,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
@@ -995,6 +1074,7 @@ czscore          call geom_to_var(nvar,coord_exp_zs(1,1))
         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
@@ -1047,6 +1127,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.'
@@ -1071,6 +1152,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.'
@@ -1107,18 +1189,7 @@ C initial geometry.
    40       continue
           endif
 #else
-          do itrial=1,100
-            itmp=1
-            call gen_rand_conf(itmp,*30)
-            goto 40
-   30       write (iout,*) 'Failed to generate random conformation',
-     &        ', itrial=',itrial
-            write (*,*) 'Failed to generate random conformation',
-     &        ', itrial=',itrial
-          enddo
-          write (iout,'(a,i3,a)') 'Processor:',me,
-     &      ' error in generating random conformation.'
-          write (*,'(a,i3,a)') 'Processor:',me,
+          write (*,'(a)') 
      &      ' error in generating random conformation.'
           stop
    40     continue
@@ -1949,12 +2020,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')
@@ -1977,6 +2054,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')
@@ -1989,6 +2070,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)
@@ -2006,6 +2089,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)
@@ -2014,6 +2101,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)
@@ -2030,6 +2119,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
@@ -2044,11 +2137,15 @@ C Get parameter filenames and open the parameter files.
       open (ielep,file=elename,status='old',readonly)
       call getenv_loc('SIDEPAR',sidename)
       open (isidep,file=sidename,status='old',readonly)
+      call getenv_loc('LIPTRANPAR',liptranname)
+      open (iliptranpar,file=liptranname,status='old',action='read')
 #ifndef CRYST_SC
       call getenv_loc('ROTPARPDB',rotname_pdb)
       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
@@ -2118,7 +2215,7 @@ c      print *,"Processor",myrank," fg_rank",fg_rank
       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
       statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
       if (lentmp.gt.0)
-     &  call copy_to_tmp(pref_orig(:ile(pref_orig))//'_'//pot(:lenpot)//
+     &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
      &    //'.stat')
       rest2name=prefix(:ilen(prefix))//'.rst'
       if(usampl) then 
@@ -2248,6 +2345,7 @@ c-------------------------------------------------------------------------------
       include 'COMMON.MD'
       open(irest2,file=rest2name,status='unknown')
       read(irest2,*) totT,EK,potE,totE,t_bath
+      totTafm=totT
       do i=1,2*nres
          read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
       enddo
@@ -2303,6 +2401,36 @@ c-------------------------------------------------------------------------------
       enddo
       return
       end
+C---------------------------------------------------------------------------
+      subroutine read_afminp
+            implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+#ifdef MPI
+      include 'mpif.h'
+#endif
+      include 'COMMON.SETUP'
+      include 'COMMON.CONTROL'
+      include 'COMMON.CHAIN'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.SBRIDGE'
+      character*320 afmcard
+      print *, "wchodze"
+      call card_concat(afmcard)
+      call readi(afmcard,"BEG",afmbeg,0)
+      call readi(afmcard,"END",afmend,0)
+      call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
+      call reada(afmcard,"VEL",velAFMconst,0.0d0)
+      print *,'FORCE=' ,forceAFMconst
+CCCC NOW PROPERTIES FOR AFM
+       distafminit=0.0d0
+       do i=1,3
+        distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
+       enddo
+        distafminit=dsqrt(distafminit)
+        print *,'initdist',distafminit
+      return
+      end
+
 c-------------------------------------------------------------------------------
       subroutine read_dist_constr
       implicit real*8 (a-h,o-z)