wham cluster ignores missing dat file for wdfab=0
[unres.git] / source / wham / src / dfa.F
index 922f933..381afa3 100644 (file)
@@ -1,14 +1,3 @@
-      block data dfadata
-      include 'DIMENSIONS'
-      include 'COMMON.INTERACT'
-      include 'COMMON.DFA'
-      DATA CK/1.0D0,1.58740105197D0,2.08008382305D0,2.51984209979D0/
-      DATA SCK/1.0D0,1.25992104989D0,1.44224957031D0,1.58740105197D0/
-      DATA S1/3.75D0,5.75D0,7.75D0,9.75D0/
-      DATA S2/4.25D0,6.25D0,8.25D0,10.25D0/
-      end
-c------------------------------------------------------------------------
-
       subroutine init_dfa_vars
 
       include 'DIMENSIONS'
@@ -82,9 +71,11 @@ C     read fragment informations
 C
       implicit real*8 (a-h,o-z)
       include 'DIMENSIONS'
+      include 'DIMENSIONS.ZSCOPT'
       include 'COMMON.IOUNITS'
       include 'COMMON.CHAIN'
       include 'COMMON.DFA'
+      include 'COMMON.FFIELD'
 
 
 C     NOTE THAT FILENAMES are FIXED, CURRENTLY!!
@@ -247,6 +238,7 @@ C            write(*,*) 'READ NEI:',i,j,fnei(i,j)
 C     END OF NEIGHBORING CA
 
 C     READ BETA RESTRAINT
+      if (wdfa_beta.eq.0.0) return
       open(iodfa, file='beta_dfa.dat',status='old',err=39)
       goto 40
  39   write(iout,'(a)') 'Error opening beta_dfa.dat file'
@@ -841,9 +833,9 @@ c            write(*,*) n1atom, j, dist
 
 c                  write(*,*) 'case1:',tmp_n
 
-                  t1dx=t1dx+0.0d0
-                  t1dy=t1dy+0.0d0
-                  t1dz=t1dz+0.0d0
+cc                  t1dx=t1dx+0.0d0
+cc                  t1dy=t1dy+0.0d0
+cc                  t1dz=t1dz+0.0d0
                   t2dx(j)=0.0d0
                   t2dy(j)=0.0d0
                   t2dz(j)=0.0d0
@@ -887,9 +879,9 @@ c     neighbor atoms
 
                   tmp_n=tmp_n+1.0d0
 c                  write(*,*) 'case4:',tmp_n
-                  t1dx=t1dx+0.0d0
-                  t1dy=t1dy+0.0d0
-                  t1dz=t1dz+0.0d0
+cc                  t1dx=t1dx+0.0d0
+cc                  t1dy=t1dy+0.0d0
+cc                  t1dz=t1dz+0.0d0
                   t2dx(j)=0.0d0
                   t2dy(j)=0.0d0
                   t2dz(j)=0.0d0
@@ -1046,8 +1038,8 @@ C     End of sheet variables
       double precision enebet
 
       enebet=0.0d0
-      bx=0.0d0;by=0.0d0;bz=0.0d0
-      shetfx=0.0d0;shetfy=0.0d0;shetfz=0.0d0
+c      bx=0.0d0;by=0.0d0;bz=0.0d0
+c      shetfx=0.0d0;shetfy=0.0d0;shetfz=0.0d0
 
       gdfab=0.0d0
 
@@ -2130,6 +2122,9 @@ c-------------------------------------------------------------------------------
       implicit none
       integer maxca
       parameter(maxca=800)
+      real*8 dfa_cutoff,dfa_cutoff_delta
+      parameter(dfa_cutoff=15.5d0)
+      parameter(dfa_cutoff_delta=0.5d0)
 cc******************************************************************************
 
 c      real*8 dfaexp(15001)
@@ -2163,6 +2158,9 @@ cc**********************************************************************
       real*8 uum, uup
       real*8 vbeta,vbetp,vbetm,y,y1,y2,yshe1,yshe2,yy1,yy2
 
+      real*8 shetfx(maxca),shetfy(maxca),shetfz(maxca)
+      common /shetf/  shetfx,shetfy,shetfz
+
       common /sheca/    bx,by,bz
       common /phys1/    inb,nmax,iselect
       common /kyori2/   dis
@@ -2181,6 +2179,7 @@ cc**********************************************************************
       common /sinsin/   sth
       
       real*8 r_pair_mat(maxca,maxca)
+      real*8 e_gcont,fprim_gcont,de_gcont
 ci      integer istrand(maxca,maxca)
 ci      integer istrand_p(maxca,maxca),istrand_m(maxca,maxca)
 ci      common  /shetest/ istrand,istrand_p,istrand_m
@@ -2201,6 +2200,12 @@ c
 
       do i=1,inb-7
          do j=i+4,inb-3
+
+            if (dis(i,j).lt.dfa_cutoff) then
+            call gcont(dis(i,j),dfa_cutoff-dfa_cutoff_delta,1.0D0,
+     &                     dfa_cutoff_delta,e_gcont,fprim_gcont)
+
+      
             ip=i+1
             ipp=i+2
             jp=j+1
@@ -2415,7 +2420,30 @@ c           uum=vbetam1(i,j)+vbetam2(i,j)
            vbet(i,j)=uup+uum
            vbetp=vbetp+uup
            vbetm=vbetm+uum
-           vbeta=vbeta+vbet(i,j)
+           vbeta=vbeta+vbet(i,j)*e_gcont
+
+              
+           if (dis(i,j) .ge. dfa_cutoff-2*dfa_cutoff_delta) then
+c gradient correction from gcont
+             de_gcont=vbet(i,j)*fprim_gcont/dis(i,j)
+             shetfx(i)=shetfx(i) + de_gcont*rx(i,j)
+             shetfy(i)=shetfy(i) + de_gcont*ry(i,j)
+             shetfz(i)=shetfz(i) + de_gcont*rz(i,j)
+
+             shetfx(j)=shetfx(j) - de_gcont*rx(i,j)
+             shetfy(j)=shetfy(j) - de_gcont*ry(i,j)
+             shetfz(j)=shetfz(j) - de_gcont*rz(i,j)
+
+c energy correction from gcont
+             vbet(i,j)=vbet(i,j)*e_gcont
+             vbetap(i,j)=vbetap(i,j)*e_gcont
+             vbetap1(i,j)=vbetap1(i,j)*e_gcont
+             vbetap2(i,j)=vbetap2(i,j)*e_gcont
+             vbetam(i,j)=vbetam(i,j)*e_gcont
+             vbetam1(i,j)=vbetam1(i,j)*e_gcont
+             vbetam2(i,j)=vbetam2(i,j)*e_gcont
+           endif
+
 
 ci         elseif(istrand(i,j).eq.0)then
 ci           vbet(i,j)=0
@@ -2434,7 +2462,16 @@ c           write(*,*) 'vbetp:',vbetp
 c           write(*,*) 'vbetm:',vbetm
 c           write(*,*) 'vbet(i,j):',vbet(i,j)
 c           stop
-
+            
+            else
+              vbetap(i,j)=0
+              vbetap1(i,j)=0
+              vbetap2(i,j)=0
+              vbetam(i,j)=0
+              vbetam1(i,j)=0
+              vbetam2(i,j)=0
+              vbet(i,j)=0
+            endif
         enddo
       enddo
 
@@ -2458,6 +2495,9 @@ c-------------------------------------------------------------------------------
       implicit none
       integer maxca
       parameter(maxca=800)
+      real*8 dfa_cutoff,dfa_cutoff_delta
+      parameter(dfa_cutoff=15.5d0)
+      parameter(dfa_cutoff_delta=0.5d0)
 cc**********************************************************************
       real*8 vbet(maxca,maxca)
       real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
@@ -2516,6 +2556,7 @@ c     local variables
       real*8  c3,v3,cc2,cc3,dmm3,dmm3__,c4,v4,c7,v7,cc7,c8,v8,cc8
       real*8  c9,v9,cc9,dmm9,dmm9__,c10,v10,dmm2,dmm1__,dmm2_1,dmm2_2
       real*8  dmm7,dmm8,dmm7__,dmm8_1,dmm8_2
+      real*8 e_gcont,fprim_gcont
 C--------------------------------------------------------------------------------
       do i=4,inb-4
          im3=i-3
@@ -2755,6 +2796,9 @@ c----------------------------------------------------------------------------
       implicit none
       integer maxca
       parameter(maxca=800)
+      real*8 dfa_cutoff,dfa_cutoff_delta
+      parameter(dfa_cutoff=15.5d0)
+      parameter(dfa_cutoff_delta=0.5d0)
 cc**********************************************************************
       real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
       real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
@@ -2795,11 +2839,17 @@ c     local variables
       real*8 y66x,y66y,y66z,yy6,yyy4,yyy5a,yyy5b
       real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z
       real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy8,yyy9a,yyy9b
+      real*8 e_gcont,fprim_gcont
 c********************************************************************************
       do i=3,inb-5
          imm=i-2
          im=i-1
          do j=i+2,inb-3
+
+            if (dis(imm,j).lt.dfa_cutoff) then
+            call gcont(dis(imm,j),dfa_cutoff-dfa_cutoff_delta,1.0D0,
+     &                     dfa_cutoff_delta,e_gcont,fprim_gcont)
+
             jp=j+1
             jpp=j+2
             
@@ -2923,7 +2973,8 @@ ci     &   .and.(istrand_p(imm,j)+istrand_m(imm,j)).ge.1) then
 !     $           -sy1*vbetam1(imm,j)-sy2*vbetam2(imm,j)
 !            shefz(i,5)=shefz(i,5)
 !     $           -sz1*vbetam1(imm,j)-sz2*vbetam2(imm,j)
-            
+
+             endif            
 ci          endif
 
          enddo
@@ -2936,6 +2987,9 @@ c--------------------------------------------------------------------------c
       implicit none
       integer maxca
       parameter(maxca=800)
+      real*8 dfa_cutoff,dfa_cutoff_delta
+      parameter(dfa_cutoff=15.5d0)
+      parameter(dfa_cutoff_delta=0.5d0)
 cc**********************************************************************
       real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
       real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
@@ -2976,11 +3030,17 @@ C     local variables
       real*8  yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z
       real*8  yy1010,yyy10,yy10,y10x,y10y,y10z,yyy8,yyy9a,yyy9b,yyy4
       real*8  yyy3a,yyy3b,y66z,yy6,yyy5a,yyy5b
+      real*8 e_gcont,fprim_gcont
 C********************************************************************************      
       do i=2,inb-6
          ip=i+1
          im=i-1
          do j=i+3,inb-3
+
+            if (dis(im,j).lt.dfa_cutoff) then
+            call gcont(dis(im,j),dfa_cutoff-dfa_cutoff_delta,1.0D0,
+     &                     dfa_cutoff_delta,e_gcont,fprim_gcont)
+
             jp=j+1
             jpp=j+2
 
@@ -3102,7 +3162,8 @@ ci     &    .and.(istrand_p(im,j)+istrand_m(im,j)).ge.1) then
 !     $           -sy1*vbetam1(im,j)-sy2*vbetam2(im,j)
 !            shefz(i,6)=shefz(i,6)
 !     $           -sz1*vbetam1(im,j)-sz2*vbetam2(im,j)
-          
+
+            endif          
 ci         endif
      
          enddo
@@ -3115,6 +3176,9 @@ c-----------------------------------------------------------------------
       implicit none
       integer maxca
       parameter(maxca=800)
+      real*8 dfa_cutoff,dfa_cutoff_delta
+      parameter(dfa_cutoff=15.5d0)
+      parameter(dfa_cutoff_delta=0.5d0)
 cc**********************************************************************
       real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
       real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
@@ -3155,12 +3219,18 @@ C     local variables
       real*8  yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y
       real*8  yy1010,yyy10,yy10,y10x,y10y,y10z,yyy4,yyy5a,yyy5b,yy6
       real*8  yyy9a,yyy9b,y5z,y66z,y9z,yyy8
+      real*8 e_gcont,fprim_gcont
 C********************************************************************************          
       
       do j=7,inb-1
          jm=j-1
          jmm=j-2
          do i=1,j-6
+
+            if (dis(i,jmm).lt.dfa_cutoff) then
+            call gcont(dis(i,jmm),dfa_cutoff-dfa_cutoff_delta,1.0D0,
+     &                     dfa_cutoff_delta,e_gcont,fprim_gcont)
+
             ip=i+1
             ipp=i+2
 
@@ -3284,6 +3354,7 @@ ci     &   .and.(istrand_p(i,jmm)+istrand_m(i,jmm)).ge.1) then
 !            shefz(j,11)=shefz(j,11)
 !     $           -sz1*vbetam1(i,jmm)-sz2*vbetam2(i,jmm)
       
+            endif
 ci         endif
          
          enddo
@@ -3296,6 +3367,9 @@ c-----------------------------------------------------------------------
       implicit none
       integer maxca
       parameter(maxca=800)
+      real*8 dfa_cutoff,dfa_cutoff_delta
+      parameter(dfa_cutoff=15.5d0)
+      parameter(dfa_cutoff_delta=0.5d0)
 cc**********************************************************************
       real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
       real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
@@ -3335,11 +3409,17 @@ C     local variables
       real*8 sx,sy,sz,sx1,sy1,sz1,sx2,sy2,sz2,y6x,y6y,y6z,y66x,y66y,y66z
       real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z
       real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy10a,yyy10b,yy6,yyy8
+      real*8 e_gcont,fprim_gcont
 !c*************************************************************************c      
       do j=6,inb-2
          jp=j+1
          jm=j-1
          do i=1,j-5
+
+            if (dis(i,jm).lt.dfa_cutoff) then
+            call gcont(dis(i,jm),dfa_cutoff-dfa_cutoff_delta,1.0D0,
+     &                     dfa_cutoff_delta,e_gcont,fprim_gcont)
+
             ip=i+1
             ipp=i+2
 
@@ -3456,6 +3536,8 @@ ci     &   .and.(istrand_p(i,jm)+istrand_m(i,jm)).ge.1) then
             shefz(j,12)=shefz(j,12)-sz*vbetam(i,jm)
      $           -sz1*vbetam1(i,jm)-sz2*vbetam2(i,jm)
       
+            endif
+           
 ci         endif
          
          ENDDO