Merge branch 'adasko' into devel
[unres.git] / source / wham / src / parmread.F
index 618bf8a..77a255f 100644 (file)
@@ -23,7 +23,6 @@ C
       include 'COMMON.FREE'
       character*1 t1,t2,t3
       character*1 onelett(4) /"G","A","P","D"/
-      character*1 toronelet(-2:2)/"p","a","G","A","P"/
       logical lprint
       dimension blower(3,3,maxlob)
       character*800 controlcard
@@ -52,10 +51,31 @@ C Assign virtual-bond length
         key = wname(i)(:ilen(wname(i)))
         call reada(controlcard,key(:ilen(key)),ww(i),1.0d0)
       enddo
-
+      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)
+c      dyn_ss=(index(controlcard,'DYN_SS').gt.0)
       write (iout,*) "iparm",iparm," myparm",myparm
-c If reading not own parameters, skip assignment
+c      do i=1,maxres
+c        dyn_ss_mask(i)=.false.
+c      enddo
+      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(me.eq.king.or..not.out1file) then
+c       print *,'indpdb=',indpdb,' pdbref=',pdbref
+c      endif
+c If reading not own parameters, skip assignment
+cc      write(iout,*) "KURWA", ww(15)
       if (iparm.eq.myparm .or. .not.separate_parset) then
 
 c
@@ -76,10 +96,12 @@ c
       wtor=ww(13)
       wtor_d=ww(14)
       wvdwpp=ww(16)
+      wstrain=ww(15)
       wbond=ww(18)
       wsccor=ww(19)
 
       endif
+cc      write(iout,*) "KURWA", wstrain,akcm,akth,wsc,dyn_ss
 
       call card_concat(controlcard,.false.)
 
@@ -102,7 +124,7 @@ c Return if not own parameters
       call reads(controlcard,"TORDPAR",tordname_t,tordname)
       open (itordp,file=tordname_t,status='old')
       rewind(itordp)
-      call reads(controlcard,"SCCORAR",sccorname_t,sccorname)
+      call reads(controlcard,"SCCORPAR",sccorname_t,sccorname)
       open (isccor,file=sccorname_t,status='old')
       rewind(isccor)
       call reads(controlcard,"FOURIER",fouriername_t,fouriername)
@@ -299,46 +321,57 @@ C
      &  ntheterm3,nsingle,ndouble
       nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
       read (ithep,*) (ithetyp(i),i=1,ntyp1)
-      do i=1,maxthetyp
-        do j=1,maxthetyp
-          do k=1,maxthetyp
-            aa0thet(i,j,k)=0.0d0
+      do i=-ntyp1,-1
+        ithetyp(i)=-ithetyp(-i)
+      enddo
+c      write (iout,*) "tu dochodze"
+      do iblock=1,2
+      do i=-maxthetyp,maxthetyp
+        do j=-maxthetyp,maxthetyp
+          do k=-maxthetyp,maxthetyp
+            aa0thet(i,j,k,iblock)=0.0d0
             do l=1,ntheterm
-              aathet(l,i,j,k)=0.0d0
+              aathet(l,i,j,k,iblock)=0.0d0
             enddo
             do l=1,ntheterm2
               do m=1,nsingle
-                bbthet(m,l,i,j,k)=0.0d0
-                ccthet(m,l,i,j,k)=0.0d0
-                ddthet(m,l,i,j,k)=0.0d0
-                eethet(m,l,i,j,k)=0.0d0
+                bbthet(m,l,i,j,k,iblock)=0.0d0
+                ccthet(m,l,i,j,k,iblock)=0.0d0
+                ddthet(m,l,i,j,k,iblock)=0.0d0
+                eethet(m,l,i,j,k,iblock)=0.0d0
               enddo
             enddo
             do l=1,ntheterm3
               do m=1,ndouble
                 do mm=1,ndouble
-                 ffthet(mm,m,l,i,j,k)=0.0d0
-                 ggthet(mm,m,l,i,j,k)=0.0d0
+                 ffthet(mm,m,l,i,j,k,iblock)=0.0d0
+                 ggthet(mm,m,l,i,j,k,iblock)=0.0d0
                 enddo
               enddo
             enddo
           enddo
         enddo
+      enddo 
       enddo
-      do i=1,nthetyp
-        do j=1,nthetyp
-          do k=1,nthetyp
-            read (ithep,'(3a)') res1,res2,res3
-            read (ithep,*) aa0thet(i,j,k)
-            read (ithep,*)(aathet(l,i,j,k),l=1,ntheterm)
+
+      do iblock=1,2
+      do i=0,nthetyp
+        do j=-nthetyp,nthetyp
+          do k=-nthetyp,nthetyp
+            read (ithep,'(6a)') res1
+            read (ithep,*) aa0thet(i,j,k,iblock)
+            read (ithep,*)(aathet(l,i,j,k,iblock),l=1,ntheterm)
             read (ithep,*)
-     &       ((bbthet(lll,ll,i,j,k),lll=1,nsingle),
-     &        (ccthet(lll,ll,i,j,k),lll=1,nsingle),
-     &        (ddthet(lll,ll,i,j,k),lll=1,nsingle),
-     &        (eethet(lll,ll,i,j,k),lll=1,nsingle),ll=1,ntheterm2)
+     &       ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
+     &        (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
+     &        (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
+     &        (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle)
+     &        ,ll=1,ntheterm2)
             read (ithep,*)
-     &      (((ffthet(llll,lll,ll,i,j,k),ffthet(lll,llll,ll,i,j,k),
-     &         ggthet(llll,lll,ll,i,j,k),ggthet(lll,llll,ll,i,j,k),
+     &      (((ffthet(llll,lll,ll,i,j,k,iblock),
+     &      ffthet(lll,llll,ll,i,j,k,iblock),
+     &         ggthet(llll,lll,ll,i,j,k,iblock)
+     &        ,ggthet(lll,llll,ll,i,j,k,iblock),
      &         llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
           enddo
         enddo
@@ -350,17 +383,53 @@ C
       do i=1,nthetyp
         do j=1,nthetyp
           do l=1,ntheterm
-            aathet(l,i,j,nthetyp+1)=aathet(l,i,j,1)
-            aathet(l,nthetyp+1,i,j)=aathet(l,1,i,j)
+            aathet(l,i,j,nthetyp+1,iblock)=0.0d0
+            aathet(l,nthetyp+1,i,j,iblock)=0.0d0
           enddo
-          aa0thet(i,j,nthetyp+1)=aa0thet(i,j,1)
-          aa0thet(nthetyp+1,i,j)=aa0thet(1,i,j)
+          aa0thet(i,j,nthetyp+1,iblock)=0.0d0
+          aa0thet(nthetyp+1,i,j,iblock)=0.0d0
         enddo
         do l=1,ntheterm
-          aathet(l,nthetyp+1,i,nthetyp+1)=aathet(l,1,i,1)
+          aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
         enddo
-        aa0thet(nthetyp+1,i,nthetyp+1)=aa0thet(1,i,1)
+        aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
+      enddo
       enddo
+C Substitution for D aminoacids from symmetry.
+      do iblock=1,2
+      do i=-nthetyp,0
+        do j=-nthetyp,nthetyp
+          do k=-nthetyp,nthetyp
+           aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
+           do l=1,ntheterm
+           aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
+           enddo
+           do ll=1,ntheterm2
+            do lll=1,nsingle
+            bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
+            ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
+            ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
+            eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
+            enddo
+          enddo
+          do ll=1,ntheterm3
+           do lll=2,ndouble
+            do llll=1,lll-1
+            ffthet(llll,lll,ll,i,j,k,iblock)=
+     &      ffthet(llll,lll,ll,-i,-j,-k,iblock)
+            ffthet(lll,llll,ll,i,j,k,iblock)=
+     &      ffthet(lll,llll,ll,-i,-j,-k,iblock)
+            ggthet(llll,lll,ll,i,j,k,iblock)=
+     &      -ggthet(llll,lll,ll,-i,-j,-k,iblock)
+            ggthet(lll,llll,ll,i,j,k,iblock)=
+     &      -ggthet(lll,llll,ll,-i,-j,-k,iblock)
+            enddo !ll
+           enddo  !lll  
+          enddo   !llll
+         enddo    !k
+        enddo     !j
+       enddo      !i
+      enddo       !iblock
 C
 C Control printout of the coefficients of virtual-bond-angle potentials
 C
@@ -372,16 +441,16 @@ C
               write (iout,'(//4a)')
      &         'Type ',onelett(i),onelett(j),onelett(k)
               write (iout,'(//a,10x,a)') " l","a[l]"
-              write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k)
+              write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
               write (iout,'(i2,1pe15.5)')
-     &           (l,aathet(l,i,j,k),l=1,ntheterm)
+     &           (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
             do l=1,ntheterm2
               write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
      &          "b",l,"c",l,"d",l,"e",l
               do m=1,nsingle
                 write (iout,'(i2,4(1pe15.5))') m,
-     &          bbthet(m,l,i,j,k),ccthet(m,l,i,j,k),
-     &          ddthet(m,l,i,j,k),eethet(m,l,i,j,k)
+     &          bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
+     &          ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
               enddo
             enddo
             do l=1,ntheterm3
@@ -390,8 +459,10 @@ C
               do m=2,ndouble
                 do n=1,m-1
                   write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
-     &              ffthet(n,m,l,i,j,k),ffthet(m,n,l,i,j,k),
-     &              ggthet(n,m,l,i,j,k),ggthet(m,n,l,i,j,k)
+     &              ffthet(n,m,l,i,j,k,iblock),
+     &              ffthet(m,n,l,i,j,k,iblock),
+     &              ggthet(n,m,l,i,j,k,iblock),
+     &              ggthet(m,n,l,i,j,k,iblock)
                 enddo
               enddo
             enddo
@@ -596,8 +667,8 @@ C
         do j=-ntortyp+1,ntortyp-1
           do k=-ntortyp+1,ntortyp-1
             read (itordp,'(3a1)') t1,t2,t3
-            if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) 
-     &        .or. t3.ne.toronelet(k)) then
+            if (t1.ne.onelett(i) .or. t2.ne.onelett(j) 
+     &        .or. t3.ne.onelett(k)) then
               write (iout,*) "Error in double torsional parameter file",
      &         i,j,k,t1,t2,t3
                stop "Error in double torsional parameter file"
@@ -683,6 +754,10 @@ CCC
 C
       read (isccor,*) nsccortyp
       read (isccor,*) (isccortyp(i),i=1,ntyp)
+      do i=-ntyp,-1
+        isccortyp(i)=-isccortyp(-i)
+       enddo
+       iscprol=isccortyp(20)
 c      write (iout,*) 'ntortyp',ntortyp
       maxinter=3
 cc maxinter is maximum interaction sites
@@ -691,12 +766,52 @@ cc maxinter is maximum interaction sites
         do j=1,nsccortyp
           read (isccor,*) nterm_sccor(i,j),nlor_sccor(i,j)
           v0ijsccor=0.0d0
+          v0ijsccor1=0.0d0
+          v0ijsccor2=0.0d0
+          v0ijsccor3=0.0d0
           si=-1.0d0
-
+          nterm_sccor(-i,j)=nterm_sccor(i,j)
+          nterm_sccor(-i,-j)=nterm_sccor(i,j)
+          nterm_sccor(i,-j)=nterm_sccor(i,j)
           do k=1,nterm_sccor(i,j)
             read (isccor,*) kk,v1sccor(k,l,i,j)
      &    ,v2sccor(k,l,i,j)
+             if (j.eq.iscprol) then
+             if (i.eq.isccortyp(10)) then
+             v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
+             v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
+             else
+              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
+     &                        +v2sccor(k,l,i,j)*dsqrt(0.75d0)
+             v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
+     &                        +v1sccor(k,l,i,j)*dsqrt(0.75d0)
+             v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
+             v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
+             v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
+             v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
+             endif
+            else
+             if (i.eq.isccortyp(10)) then
+             v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
+             v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
+             else
+               if (j.eq.isccortyp(10)) then
+             v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
+             v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
+               else
+             v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
+             v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
+             v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
+             v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
+             v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
+             v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
+                endif
+               endif
+             endif
             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
+            v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
+            v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
+            v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
             si=-si
           enddo
           do k=1,nlor_sccor(i,j)
@@ -705,7 +820,10 @@ cc maxinter is maximum interaction sites
             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
      &(1+vlor3sccor(k,i,j)**2)
           enddo
-          v0sccor(i,j)=v0ijsccor
+          v0sccor(l,i,j)=v0ijsccor
+          v0sccor(l,-i,j)=v0ijsccor1
+          v0sccor(l,i,-j)=v0ijsccor2
+          v0sccor(l,-i,-j)=v0ijsccor3   
         enddo
       enddo
       enddo
@@ -1021,8 +1139,27 @@ c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
         enddo
       enddo
 C
-C Define the SC-p interaction constants
+C Define the SC-p interaction constants and SS bond potentials
 C
+      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
+      else
+        ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
+      endif
+       write (iout,*) "Parameters of the SS-bond potential:"
+       write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
+     & " AKCT",akct
+       write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
+       write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
+       write (iout,*)" HT",Ht
+
 #ifdef OLDSCP
       do i=1,20
 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates 
@@ -1071,7 +1208,7 @@ C
 C
 C Define the constants of the disulfide bridge
 C
-      ebr=-5.50D0
+c      ebr=-5.50D0
 c
 c Old arbitrary potential - commented out.
 c
@@ -1082,21 +1219,21 @@ 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
 
-      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      if (lprint) 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
       return
       end