correction in wham and UNRES for lipid and correlation
[unres.git] / source / wham / src-M / parmread.F
index 914c090..8db78b9 100644 (file)
@@ -21,6 +21,8 @@ C
       include 'COMMON.SCCOR'
       include 'COMMON.SCROT'
       include 'COMMON.FREE'
+      include 'COMMON.SHIELD'
+      include 'COMMON.CONTROL'
       character*1 t1,t2,t3
       character*1 onelett(4) /"G","A","P","D"/
       character*1 toronelet(-2:2) /"p","a","G","A","P"/
@@ -35,6 +37,8 @@ C
       character*16 key
       integer iparm
       double precision ip,mp
+      character*6 res1
+C      write (iout,*) "KURWA"
 C
 C Body
 C
@@ -55,6 +59,66 @@ C Assign virtual-bond length
 
       write (iout,*) "iparm",iparm," myparm",myparm
 c If reading not own parameters, skip assignment
+      call reada(controlcard,"D0CM",d0cm,3.78d0)
+      call reada(controlcard,"AKCM",akcm,15.1d0)
+      call reada(controlcard,"AKTH",akth,11.0d0)
+      call reada(controlcard,"AKCT",akct,12.0d0)
+      call reada(controlcard,"V1SS",v1ss,-1.08d0)
+      call reada(controlcard,"V2SS",v2ss,7.61d0)
+      call reada(controlcard,"V3SS",v3ss,13.7d0)
+      call reada(controlcard,"EBR",ebr,-5.50D0)
+      call reada(controlcard,"DTRISS",dtriss,1.0D0)
+      call reada(controlcard,"ATRISS",atriss,0.3D0)
+      call reada(controlcard,"BTRISS",btriss,0.02D0)
+      call reada(controlcard,"CTRISS",ctriss,1.0D0)
+      dyn_ss=(index(controlcard,'DYN_SS').gt.0)
+      write(iout,*) "ATRISS",atriss
+      write(iout,*) "BTRISS",btriss
+      write(iout,*) "CTRISS",ctriss
+      write(iout,*) "DTRISS",dtriss
+
+C      do i=1,maxres
+C        dyn_ss_mask(i)=.false.
+C      enddo
+C      ebr=-12.0D0
+c
+c Old arbitrary potential - commented out.
+c
+c      dbr= 4.20D0
+c      fbr= 3.30D0
+c
+c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
+c energy surface of diethyl disulfide.
+c A. Liwo and U. Kozlowska, 11/24/03
+c
+      D0CM = 3.78d0
+      AKCM = 15.1d0
+      AKTH = 11.0d0
+      AKCT = 12.0d0
+      V1SS =-1.08d0
+      V2SS = 7.61d0
+      V3SS = 13.7d0
+
+      do i=1,maxres-1
+        do j=i+1,maxres
+          dyn_ssbond_ij(i,j)=1.0d300
+        enddo
+      enddo
+      call reada(controlcard,"HT",Ht,0.0D0)
+C      if (dyn_ss) then
+C        ss_depth=ebr/wsc-0.25*eps(1,1)
+C        write(iout,*) HT,wsc,eps(1,1),'KURWA'
+C        Ht=Ht/wsc-0.25*eps(1,1)
+       
+C        akcm=akcm*whpb/wsc
+C        akth=akth*whpb/wsc
+C        akct=akct*whpb/wsc
+C        v1ss=v1ss*whpb/wsc
+C        v2ss=v2ss*whpb/wsc
+C        v3ss=v3ss*whpb/wsc
+C      else
+C        ss_depth=ebr/whpb-0.25*eps(1,1)*wsc/whpb
+C      endif
 
       if (iparm.eq.myparm .or. .not.separate_parset) then
 
@@ -78,7 +142,10 @@ c
       wvdwpp=ww(16)
       wbond=ww(18)
       wsccor=ww(19)
+      whpb=ww(15)
+      wstrain=ww(15)
       wliptran=ww(22)
+      wshield=ww(25)
       endif
 
       call card_concat(controlcard,.false.)
@@ -300,7 +367,7 @@ C
 C Read the parameters of Utheta determined from ab initio surfaces
 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
 C
-c      write (iout,*) "tu dochodze"
+      write (iout,*) "tu dochodze"
       read (ithep,*) nthetyp,ntheterm,ntheterm2,
      &  ntheterm3,nsingle,ndouble
       nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
@@ -308,7 +375,7 @@ c      write (iout,*) "tu dochodze"
       do i=-ntyp1,-1
         ithetyp(i)=-ithetyp(-i)
       enddo
-c      write (iout,*) "tu dochodze"
+      write (iout,*) "tu dochodze"
       do iblock=1,2
       do i=-maxthetyp,maxthetyp
         do j=-maxthetyp,maxthetyp
@@ -337,11 +404,13 @@ c      write (iout,*) "tu dochodze"
         enddo
       enddo
       enddo
+C      write (iout,*) "KURWA1"
       do iblock=1,2
       do i=0,nthetyp
         do j=-nthetyp,nthetyp
           do k=-nthetyp,nthetyp
             read (ithep,'(6a)') res1
+            write(iout,*) res1,i,j,k
             read (ithep,*) aa0thet(i,j,k,iblock)
             read (ithep,*)(aathet(l,i,j,k,iblock),l=1,ntheterm)
             read (ithep,*)
@@ -359,6 +428,7 @@ c      write (iout,*) "tu dochodze"
           enddo
         enddo
       enddo
+C       write(iout,*) "KURWA1.1"
 C
 C For dummy ends assign glycine-type coefficients of theta-only terms; the
 C coefficients of theta-and-gamma-dependent terms are zero.
@@ -378,6 +448,7 @@ C
         aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
       enddo
       enddo
+C       write(iout,*) "KURWA1.5"
 C Substitution for D aminoacids from symmetry.
       do iblock=1,2
       do i=-nthetyp,0
@@ -456,7 +527,7 @@ C
       call flush(iout)
       endif
 #endif
-
+C      write(iout,*) 'KURWA2'
 #ifdef CRYST_SC
 C
 C Read the parameters of the probability distribution/energy expression
@@ -564,6 +635,7 @@ C
       enddo
 #endif
       close(irotam)
+C      write (iout,*) 'KURWAKURWA'
 #ifdef CRYST_TOR
 C
 C Read torsional parameters in old format
@@ -1040,7 +1112,8 @@ C---------------------- GB or BP potential -----------------------------
       read (isidep,*)(alp(i),i=1,ntyp)
       do i=1,ntyp
        read (isidep,*)(epslip(i,j),j=i,ntyp)
-C       print *,"WARNING!!"
+C       write(iout,*) "WARNING!!",i,ntyp
+       write(iout,*) "epslip", i, (epslip(i,j),j=i,ntyp)
 C       do j=1,ntyp
 C       epslip(i,j)=epslip(i,j)+0.05d0
 C       enddo
@@ -1211,7 +1284,7 @@ C
 C
 C Define the constants of the disulfide bridge
 C
-      ebr=-5.50D0
+C      ebr=-12.0D0
 c
 c Old arbitrary potential - commented out.
 c
@@ -1222,21 +1295,53 @@ c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
 c energy surface of diethyl disulfide.
 c A. Liwo and U. Kozlowska, 11/24/03
 c
-      D0CM = 3.78d0
-      AKCM = 15.1d0
-      AKTH = 11.0d0
-      AKCT = 12.0d0
-      V1SS =-1.08d0
-      V2SS = 7.61d0
-      V3SS = 13.7d0
+C      D0CM = 3.78d0
+C      AKCM = 15.1d0
+C      AKTH = 11.0d0
+C      AKCT = 12.0d0
+C      V1SS =-1.08d0
+C      V2SS = 7.61d0
+C      V3SS = 13.7d0
+      write (iout,*) dyn_ss,'dyndyn'
+      if (dyn_ss) then
+        ss_depth=ebr/wsc-0.25*eps(1,1)
+C        write(iout,*) akcm,whpb,wsc,'KURWA'
+        Ht=Ht/wsc-0.25*eps(1,1)
 
-      if (lprint) then
+        akcm=akcm*whpb/wsc
+        akth=akth*whpb/wsc
+        akct=akct*whpb/wsc
+        v1ss=v1ss*whpb/wsc
+        v2ss=v2ss*whpb/wsc
+        v3ss=v3ss*whpb/wsc
+      else
+        ss_depth=ebr/whpb-0.25*eps(1,1)*wsc/whpb
+      endif
+
+C      if (lprint) then
       write (iout,'(/a)') "Disulfide bridge parameters:"
       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
      & ' v3ss:',v3ss
-      endif
+C      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 
       return
       end