Merge branch 'lipid' into AFM
[unres.git] / source / unres / src_MD-M / readrtns_CSA.F
index 889e5a5..ebb59da 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
@@ -98,6 +99,10 @@ c      print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
 C Set up the time limit (caution! The time must be input in minutes!)
       read_cart=index(controlcard,'READ_CART').gt.0
       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
+C this variable with_theta_constr is the variable which allow to read and execute the
+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
       unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
@@ -141,6 +146,15 @@ C Set up the time limit (caution! The time must be input in minutes!)
       selfguide=(index(controlcard,'SELFGUIDE'))
       print *,'AFMlog',AFMlog,selfguide,"KUPA"
       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,'MAXGEN',maxgen,10000)
       call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
       call readi(controlcard,"KDIAG",kdiag,0)
@@ -244,8 +258,24 @@ C      endif
       write(iout,*) "bordlipbot=",bordlipbot
       write(iout,*) "bufliptop=",bufliptop
       write(iout,*) "buflipbot=",buflipbot
-
-
+      write (iout,*) "SHIELD MODE",shield_mode
+      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,7 +601,7 @@ C
       integer rescode
       double precision x(maxvar)
       character*256 pdbfile
-      character*320 weightcard
+      character*400 weightcard
       character*80 weightcard_t,ucase
       dimension itype_pdb(maxres)
       common /pizda/ itype_pdb
@@ -713,6 +743,14 @@ C 12/1/95 Added weight for the multi-body term WCORR
       call reada(weightcard,"V2SS",v2ss,7.61d0)
       call reada(weightcard,"V3SS",v3ss,13.7d0)
       call reada(weightcard,"EBR",ebr,-5.50D0)
+      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,*) "ATRISS=", atriss
+      write (iout,*) "BTRISS=", btriss
+      write (iout,*) "CTRISS=", ctriss
+      write (iout,*) "DTRISS=", dtriss
       dyn_ss=(index(weightcard,'DYN_SS').gt.0)
       do i=1,maxres
         dyn_ss_mask(i)=.false.
@@ -733,7 +771,11 @@ C 12/1/95 Added weight for the multi-body term WCORR
         v2ss=v2ss*wstrain/wsc
         v3ss=v3ss*wstrain/wsc
       else
-        ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
+        if (wstrain.ne.0.0) then
+         ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
+        else
+          ss_depth=0.0
+        endif
       endif
 
       if(me.eq.king.or..not.out1file) then
@@ -755,9 +797,11 @@ C 12/1/95 Added weight for the multi-body term WCORR
   33    write (iout,'(a)') 'Error opening PDB file.'
         stop
   34    continue
-c        print *,'Begin reading pdb data'
+c        write (iout,*) 'Begin reading pdb data'
+c        call flush(iout)
         call readpdb
-c        print *,'Finished reading pdb data'
+c        write (iout,*) 'Finished reading pdb data'
+c        call flush(iout)
         if(me.eq.king.or..not.out1file)
      &   write (iout,'(a,i3,a,i3)')'nsup=',nsup,
      &   ' nstart_sup=',nstart_sup
@@ -847,27 +891,78 @@ C 8/13/98 Set limits to generating the dihedral angles
       enddo
       read (inp,*) ndih_constr
       if (ndih_constr.gt.0) then
-        read (inp,*) ftors
-        read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
+C        read (inp,*) ftors
+        read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
+     &  i=1,ndih_constr)
         if(me.eq.king.or..not.out1file)then
          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
         endif
         do i=1,ndih_constr
           phi0(i)=deg2rad*phi0(i)
           drange(i)=deg2rad*drange(i)
         enddo
-        if(me.eq.king.or..not.out1file)
-     &   write (iout,*) 'FTORS',ftors
+C        if(me.eq.king.or..not.out1file)
+C     &   write (iout,*) 'FTORS',ftors
         do i=1,ndih_constr
           ii = idih_constr(i)
           phibound(1,ii) = phi0(i)-drange(i)
           phibound(2,ii) = phi0(i)+drange(i)
         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 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
+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 
+        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
+        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
+C
+C      with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
+C      write (iout,*) "with_dihed_constr ",with_dihed_constr
       nnt=1
 #ifdef MPI
       if (me.eq.king) then
@@ -954,7 +1049,9 @@ czscore          call geom_to_var(nvar,coord_exp_zs(1,1))
           enddo
           call contact(.true.,ncont_ref,icont_ref,co)
         endif
-c        write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
+        endif
+        print *, "A TU"
+        write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
         call flush(iout)
         if (constr_dist.gt.0) call read_dist_constr
         write (iout,*) "After read_dist_constr nhpb",nhpb
@@ -975,7 +1072,7 @@ c        write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
      &     restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
         enddo
         endif
-      endif
+C      endif
       if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
      &    .and. modecalc.ne.8 .and. modecalc.ne.9 .and. 
      &    modecalc.ne.10) then
@@ -1035,6 +1132,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.'
@@ -1147,6 +1245,7 @@ cd    write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
      &  write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)') 
      &  'Processor',myrank,': end reading molecular data.'
 #endif
+      print *,"A TU?"
       return
       end
 c--------------------------------------------------------------------------
@@ -2303,7 +2402,8 @@ c-------------------------------------------------------------------------------
       integer ifrag_(2,100),ipair_(2,100)
       double precision wfrag_(100),wpair_(100)
       character*500 controlcard
-c      write (iout,*) "Calling read_dist_constr"
+      print *, "WCHODZE"
+      write (iout,*) "Calling read_dist_constr"
 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
 c      call flush(iout)
       call card_concat(controlcard)
@@ -2397,11 +2497,30 @@ c            write (iout,*) "j",j," k",k
         enddo
         endif
       enddo 
+      print *,ndist_
       do i=1,ndist_
-        read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
+        if (constr_dist.eq.11) then
+        read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
+     &     ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
+        fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
+        else
+C        print *,"in else"
+        read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
+     &     ibecarb(i),forcon(nhpb+1)
+        endif
         if (forcon(nhpb+1).gt.0.0d0) then
           nhpb=nhpb+1
-          dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
+          if (ibecarb(i).gt.0) then
+            ihpb(i)=ihpb(i)+nres
+            jhpb(i)=jhpb(i)+nres
+          endif
+          if (dhpb(nhpb).eq.0.0d0)
+     &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
+        endif
+C        read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
+C        if (forcon(nhpb+1).gt.0.0d0) then
+C          nhpb=nhpb+1
+C          dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
 #ifdef MPI
           if (.not.out1file .or. me.eq.king)
      &    write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
@@ -2410,7 +2529,7 @@ c            write (iout,*) "j",j," k",k
           write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 #endif
-        endif
+
       enddo
       call flush(iout)
       return