projects / unres.git / blobdiff
? search:
summary | shortlog | log | commit | commitdiff | tree
raw | inline | side by side
Merge branch 'lipid' into AFM
[unres.git] / source / unres / src_MD-M / parmread.F
diff --git a/source/unres/src_MD-M/parmread.F b/source/unres/src_MD-M/parmread.F
index 183a917..ae4d710 100644 (file)
--- a/source/unres/src_MD-M/parmread.F
+++ b/source/unres/src_MD-M/parmread.F
@@ -26,6 +26,8 @@ C
       include 'COMMON.SBRIDGE'
       include 'COMMON.MD'
       include 'COMMON.SETUP'
+      include 'COMMON.CONTROL'
+      include 'COMMON.SHIELD'
       character*1 t1,t2,t3
       character*1 onelett(4) /"G","A","P","D"/
       character*1 toronelet(-2:2) /"p","a","G","A","P"/
@@ -73,6 +75,7 @@ c
 #else
       read (ibond,*) junk,vbldp0,vbldpdum,akp,rjunk,mp,ip,pstok
       do i=1,ntyp
+      print *,i
         read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
      &   j=1,nbondterm(i)),msc(i),isc(i),restok(i)
         dsc(i) = vbldsc0(1,i)
@@ -944,6 +947,16 @@ c        B2(1,i)  = b(2)
 c        B2(2,i)  = b(4)
 c        B2(1,-i)  =b(2)
 c        B2(2,-i)  =-b(4)
+        B1tilde(1,i) = b(3,i)
+        B1tilde(2,i) =-b(5,i)
+C        B1tilde(1,-i) =-b(3,i)
+C        B1tilde(2,-i) =b(5,i)
+        b1tilde(1,i)=0.0d0
+        b1tilde(2,i)=0.0d0
+        B2(1,i)  = b(2,i)
+        B2(2,i)  = b(4,i)
+C        B2(1,-i)  =b(2,i)
+C        B2(2,-i)  =-b(4,i)
 
 c        b2(1,i)=0.0d0
 c        b2(2,i)=0.0d0
@@ -1082,7 +1095,7 @@ C
      & ', exponents are ',expon,2*expon 
       goto (10,20,30,30,40) ipot
 C----------------------- LJ potential ---------------------------------
-   10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
+   10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
      &   (sigma0(i),i=1,ntyp)
       if (lprint) then
        write (iout,'(/a/)') 'Parameters of the LJ potential:'
@@ -1094,7 +1107,7 @@ C----------------------- LJ potential ---------------------------------
       endif
       goto 50
 C----------------------- LJK potential --------------------------------
-   20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
+   20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
       if (lprint) then
        write (iout,'(/a/)') 'Parameters of the LJK potential:'
@@ -1108,7 +1121,7 @@ C----------------------- LJK potential --------------------------------
       goto 50
 C---------------------- GB or BP potential -----------------------------
    30 do i=1,ntyp
-       read (isidep,*,end=116,err=116)(eps(i,j),j=i,ntyp)
+       read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
       enddo
       read (isidep,*,end=116,err=116)(sigma0(i),i=1,ntyp)
       read (isidep,*,end=116,err=116)(sigii(i),i=1,ntyp)
@@ -1141,7 +1154,7 @@ C For the GB potential convert sigma'**2 into chi'
       endif
       goto 50
 C--------------------- GBV potential -----------------------------------
-   40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
+   40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
       if (lprint) then
@@ -1290,7 +1303,7 @@ c      lprint=.false.
 C
 C Define the constants of the disulfide bridge
 C
-      ebr=-5.50D0
+C      ebr=-12.00D0
 c
 c Old arbitrary potential - commented out.
 c
@@ -1301,13 +1314,13 @@ 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
 c      akcm=0.0d0
 c      akth=0.0d0
 c      akct=0.0d0
@@ -1315,14 +1328,33 @@ c      v1ss=0.0d0
 c      v2ss=0.0d0
 c      v3ss=0.0d0
       
-      if(me.eq.king) 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      if(me.eq.king) then
+C      write (iout,'(/a)') "Disulfide bridge parameters:"
+C      write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
+C      write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
+C      write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
+C      write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
+C     &  ' v3ss:',v3ss
+C      endif
+C set the variables used for shielding effect
+C      write (iout,*) "SHIELD MODE",shield_mode
+C      if (shield_mode.gt.0) then
+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
+C      VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
+C      VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
+C      write (iout,*) VSolvSphere,VSolvSphere_div
+C long axis of side chain 
+C      do i=1,ntyp
+C      long_r_sidechain(i)=vbldsc0(1,i)
+C      short_r_sidechain(i)=sigma0(i)
+C      enddo
+C lets set the buffor value
+C      buff_shield=1.0d0
+C      endif
       return
   111 write (iout,*) "Error reading bending energy parameters."
       goto 999
@@ -1394,6 +1426,22 @@ c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
 #else
       call getenv(var,val)
 #endif
-
+C set the variables used for shielding effect
+C      if (shield_mode.gt.0) then
+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
+C      VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
+C      VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
+C long axis of side chain 
+C      do i=1,ntyp
+C      long_r_sidechain(i)=vbldsc0(1,i)
+C      short_r_sidechain(i)=sigma0(i)
+C      enddo
+C lets set the buffor value
+C      buff_shield=1.0d0
+C      endif
       return
       end
UNRESPACK 3.x