Adam's changes to wham and cluster following previous commit
[unres.git] / source / wham / src / energy_p_new.F
index 816e38e..81abea5 100644 (file)
@@ -2,6 +2,7 @@
       implicit real*8 (a-h,o-z)
       include 'DIMENSIONS'
       include 'DIMENSIONS.ZSCOPT'
+      include 'DIMENSIONS.FREE'
 
 #ifndef ISNAN
       external proc_proc
@@ -133,7 +134,7 @@ c      write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t
      & +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4
      & +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6
      & +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d
-     & +wbond*estr+wsccor*fact(1)*esccor+ehomology_constr
+     & +wbond*estr+wsccor*fact(1)*esccor!+ehomology_constr
      & +wdfa_dist*edfadis+wdfa_tor*edfator+wdfa_nei*edfanei
      & +wdfa_beta*edfabet
 #else
@@ -144,7 +145,7 @@ c      write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t
      & +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4
      & +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6
      & +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d
-     & +wbond*estr+wsccor*fact(1)*esccor+ehomology_constr
+     & +wbond*estr+wsccor*fact(1)*esccor!+ehomology_constr
      & +wdfa_dist*edfadis+wdfa_tor*edfator+wdfa_nei*edfanei
      & +wdfa_beta*edfabet
 #endif
@@ -1820,6 +1821,7 @@ C
       implicit real*8 (a-h,o-z)
       include 'DIMENSIONS'
       include 'DIMENSIONS.ZSCOPT'
+      include 'DIMENSIONS.FREE'
       include 'COMMON.CONTROL'
       include 'COMMON.IOUNITS'
       include 'COMMON.GEO'
@@ -3149,7 +3151,7 @@ c MODELLER restraint function
       implicit real*8 (a-h,o-z)
       include 'DIMENSIONS'
       include 'DIMENSIONS.ZSCOPT'
-
+      include 'DIMENSIONS.FREE'
       integer nnn, i, j, k, ki, irec, l
       integer katy, odleglosci, test7
       real*8 odleg, odleg2, odleg3, kat, kat2, kat3, gdih(max_template)
@@ -3180,7 +3182,7 @@ c
       include 'COMMON.SETUP'
       include 'COMMON.NAMES'
 
-      do i=1,19
+      do i=1,max_template
         distancek(i)=9999999.9
       enddo
 
@@ -3510,8 +3512,6 @@ c
 c        sgtheta=-gtheta(k)*theta_diff(k)*sigma_theta(k,i)*waga_theta ! right functional form?
           sum_sgtheta=sum_sgtheta+sgtheta ! cum variable
         enddo
-c       grad_theta3=sum_sgtheta/sum_gtheta 1/*theta(i)? s. line below
-c       grad_theta3=sum_sgtheta/sum_gtheta
 c
 c       Final value of gradient using same var as in Econstr_back
         dutheta(i-2)=sum_sgtheta/sum_gtheta*waga_theta
 c
 c          For Gaussian-type Urestr
 c
-        ehomology_constr=(waga_dist*odleg+waga_angle*kat+
-     &              waga_theta*Eval+waga_d*Erot)*waga_homology(iset)
+c        ehomology_constr=(waga_dist*odleg+waga_angle*kat+
+c     &              waga_theta*Eval+waga_d*Erot)*waga_homology(iset)
+        ehomology_constr=waga_dist*odleg+waga_angle*kat+
+     &              waga_theta*Eval+waga_d*Erot
 c     write (iout,*) "ehomology_constr=",ehomology_constr
       else
 c
 c          For Lorentzian-type Urestr
 c  
-        ehomology_constr=(-waga_dist*odleg+waga_angle*kat+
-     &              waga_theta*Eval+waga_d*Erot)*waga_homology(iset)
+c        ehomology_constr=(-waga_dist*odleg+waga_angle*kat+
+c     &              waga_theta*Eval+waga_d*Erot)*waga_homology(iset)
+        ehomology_constr=-waga_dist*odleg+waga_angle*kat+
+     &              waga_theta*Eval+waga_d*Erot
 c     write (iout,*) "ehomology_constr=",ehomology_constr
       endif
-c     write (iout,*) "odleg",odleg," kat",kat," Uconst_back",Uconst_back
-c     write (iout,*) "ehomology_constr",ehomology_constr
-c     ehomology_constr=odleg+kat+Uconst_back
+#ifdef DEBUG
+      write (iout,*) "odleg",waga_dist,odleg," kat",waga_angle,kat,
+     & "Eval",waga_theta,eval,
+     &   "Erot",waga_d,Erot
+      write (iout,*) "ehomology_constr",ehomology_constr
+#endif
       return
 
   748 format(a8,f12.3,a6,f12.3,a7,f12.3)
@@ -3739,6 +3746,7 @@ c
       implicit real*8 (a-h,o-z)
       include 'DIMENSIONS'
       include 'DIMENSIONS.ZSCOPT'
+      include 'DIMENSIONS.FREE'
       include 'COMMON.LOCAL'
       include 'COMMON.GEO'
       include 'COMMON.INTERACT'
@@ -4065,6 +4073,7 @@ C
       implicit real*8 (a-h,o-z)
       include 'DIMENSIONS'
       include 'DIMENSIONS.ZSCOPT'
+      include 'DIMENSIONS.FREE'
       include 'COMMON.LOCAL'
       include 'COMMON.GEO'
       include 'COMMON.INTERACT'
@@ -4129,7 +4138,8 @@ c      write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1)
           enddo
         else
           phii1=0.0d0
-          ityp3=nthetyp+1
+c          ityp3=nthetyp+1
+          ityp3=ithetyp(itype(i))
           do k=1,nsingle
             cosph2(k)=0.0d0
             sinph2(k)=0.0d0
@@ -4545,6 +4555,7 @@ C
       implicit real*8 (a-h,o-z)
       include 'DIMENSIONS'
       include 'DIMENSIONS.ZSCOPT'
+      include 'DIMENSIONS.FREE'
       include 'COMMON.GEO'
       include 'COMMON.LOCAL'
       include 'COMMON.VAR'
@@ -5196,6 +5207,7 @@ c        amino-acid residues.
       implicit real*8 (a-h,o-z)
       include 'DIMENSIONS'
       include 'DIMENSIONS.ZSCOPT'
+      include 'DIMENSIONS.FREE'
       include 'COMMON.VAR'
       include 'COMMON.GEO'
       include 'COMMON.LOCAL'
@@ -5249,6 +5261,9 @@ c   3 = SC...Ca...Ca...SCi
           cosphi=dcos(j*tauangle(intertyp,i))
           sinphi=dsin(j*tauangle(intertyp,i))
           esccor=esccor+v1ij*cosphi+v2ij*sinphi
+#ifdef DEBUG
+          esccor_ii=esccor_ii+v1ij*cosphi+v2ij*sinphi
+#endif
           gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi)
         enddo
         gloc_sc(intertyp,i-3,icg)=gloc_sc(intertyp,i-3,icg)+wsccor*gloci
@@ -5261,6 +5276,9 @@ c     &gloc_sc(intertyp,i-3,icg)
      & ,(v2sccor(j,intertyp,itori,itori1),j=1,6)
         gsccor_loc(i-3)=gsccor_loc(i-3)+gloci
        enddo !intertyp
+#ifdef DEBUG
+       write (iout,*) "i",i,(tauangle(j,i),j=1,3),esccor_ii
+#endif
       enddo
 c        do i=1,nres
 c        write (iout,*) "W@T@F",  gloc_sc(1,i,icg),gloc(i,icg)
@@ -5582,6 +5600,10 @@ c     &         ' jj=',jj,' kk=',kk
 C Contacts I-J and (I+1)-(J+1) or (I+1)-(J-1) occur simultaneously. 
 C The system gains extra energy.
               ecorr=ecorr+ehbcorr(i,j,i+1,j1,jj,kk,0.72D0,0.32D0)
+#ifdef DEBUG
+              write (iout,*) "ecorr",i,j,i+1,j1,
+     &               ehbcorr(i,j,i+1,j1,jj,kk,0.72D0,0.32D0)
+#endif
               n_corr=n_corr+1
             else if (j1.eq.j) then
 C Contacts I-J and I-(J+1) occur simultaneously. 
@@ -5869,11 +5891,11 @@ cd    ees0pkl=0.0D0
 cd    ees0pij=1.0D0
 cd    ees0mkl=0.0D0
 cd    ees0mij=1.0D0
-c     write (iout,*)'Contacts have occurred for peptide groups',i,j,
-c    &   ' and',k,l
-c     write (iout,*)'Contacts have occurred for peptide groups',
-c    &  i,j,' fcont:',eij,' eij',' eesij',ees0pij,ees0mij,' and ',k,l
-c    & ,' fcont ',ekl,' eeskl',ees0pkl,ees0mkl,' ees=',ees
+cd      write (iout,*)'Contacts have occurred for peptide groups',i,j,
+cd     &   ' and',k,l
+cd      write (iout,*)'Contacts have occurred for peptide groups',
+cd     &  i,j,' fcont:',eij,' eij',' eesij',ees0pij,ees0mij,' and ',k,l
+cd     & ,' fcont ',ekl,' eeskl',ees0pkl,ees0mkl,' ees=',ees
 C Calculate the multi-body contribution to energy.
       ecorr=ecorr+ekont*ees
       if (calc_grad) then