unres src_MD-M from branch Multichain
[unres.git] / source / unres / src_MD / energy_p_new_barrier.F
index cb97c0a..b76681e 100644 (file)
@@ -99,6 +99,12 @@ c      if (modecalc.eq.12.or.modecalc.eq.14) then
 c        call int_from_cart1(.false.)
 c      endif
 #endif     
+#ifndef DFA
+      edfadis=0.0d0
+      edfator=0.0d0
+      edfanei=0.0d0
+      edfabet=0.0d0
+#endif
 #ifdef TIMING
 #ifdef MPI
       time00=MPI_Wtime()
@@ -132,6 +138,7 @@ C
 C Calculate electrostatic (H-bonding) energy of the main chain.
 C
   107 continue
+#ifdef DFA
 C     BARTEK for dfa test!
       if (wdfa_dist.gt.0) then 
         call edfad(edfadis)
@@ -156,6 +163,7 @@ c      print*, 'edfan is finished!', edfanei
       else
         edfabet=0
       endif
+#endif
 c      print*, 'edfab is finished!', edfabet
 cmc
 cmc Sep-06: egb takes care of dynamic ss bonds too
@@ -543,6 +551,7 @@ cMS$ATTRIBUTES C ::  proc_proc
       include 'COMMON.TIME1'
       include 'COMMON.MAXGRAD'
       include 'COMMON.SCCOR'
+      include 'COMMON.MD'
 #ifdef TIMING
 #ifdef MPI
       time01=MPI_Wtime()
@@ -813,6 +822,14 @@ c      enddo
 #endif
         enddo
       enddo 
+      if (constr_homology.gt.0) then
+        do i=1,nct
+          do j=1,3
+            gradc(j,i,icg)=gradc(j,i,icg)+duscdiff(j,i)
+            gradx(j,i,icg)=gradx(j,i,icg)+duscdiffx(j,i)
+          enddo
+        enddo
+      endif
 #ifdef DEBUG
       write (iout,*) "gloc before adding corr"
       do i=1,4*nres
@@ -4862,6 +4879,8 @@ C
      & sinph1ph2(maxdouble,maxdouble)
       logical lprn /.false./, lprn1 /.false./
       etheta=0.0D0
+c      write (iout,*) "EBEND ithet_start",ithet_start,
+c     &     " ithet_end",ithet_end
       do i=ithet_start,ithet_end
         if ((itype(i-1).eq.ntyp1).or.(itype(i-2).eq.ntyp1).or.
      &(itype(i).eq.ntyp1)) cycle
@@ -5851,7 +5870,7 @@ C Proline-Proline pair is a special case...
 c------------------------------------------------------------------------------
 c LICZENIE WIEZOW Z ROWNANIA ENERGII MODELLERA
       subroutine e_modeller(ehomology_constr)
-      ehomology_constr=0.0
+      ehomology_constr=0.0d0
       write (iout,*) "!!!!!UWAGA, JESTEM W DZIWNEJ PETLI, TEST!!!!!"
       return
       end
@@ -6019,12 +6038,25 @@ c     write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
          dij=dist(i,j)
 c        write (iout,*) "dij(",i,j,") =",dij
          do k=1,constr_homology
+c           write(iout,*) ii,k,i,j,l_homo(k,ii),dij,odl(k,ii)
+           if(.not.l_homo(k,ii)) cycle
            distance(k)=odl(k,ii)-dij
 c          write (iout,*) "distance(",k,") =",distance(k)
+c
+c          For Gaussian-type Urestr
+c
            distancek(k)=0.5d0*distance(k)**2*sigma_odl(k,ii) ! waga_dist rmvd from Gaussian argument
 c          write (iout,*) "sigma_odl(",k,ii,") =",sigma_odl(k,ii)
 c          write (iout,*) "distancek(",k,") =",distancek(k)
 c          distancek(k)=0.5d0*waga_dist*distance(k)**2*sigma_odl(k,ii)
+c
+c          For Lorentzian-type Urestr
+c
+           if (waga_dist.lt.0.0d0) then
+              sigma_odlir(k,ii)=dsqrt(1/sigma_odl(k,ii))
+              distancek(k)=distance(k)**2/(sigma_odlir(k,ii)*
+     &                     (distance(k)**2+sigma_odlir(k,ii)**2))
+           endif
          enddo
          
          min_odl=minval(distancek)
@@ -6040,8 +6072,19 @@ c        write (iout,* )"min_odl",min_odl
 c Nie wiem po co to liczycie jeszcze raz!
 c            odleg3=-waga_dist(iset)*((distance(i,j,k)**2)/ 
 c     &              (2*(sigma_odl(i,j,k))**2))
+           if(.not.l_homo(k,ii)) cycle
+           if (waga_dist.ge.0.0d0) then
+c
+c          For Gaussian-type Urestr
+c
             godl(k)=dexp(-distancek(k)+min_odl)
             odleg2=odleg2+godl(k)
+c
+c          For Lorentzian-type Urestr
+c
+           else
+            odleg2=odleg2+distancek(k)
+           endif
 
 ccc       write(iout,779) i,j,k, "odleg2=",odleg2, "odleg3=", odleg3,
 ccc     & "dEXP(odleg3)=", dEXP(odleg3),"distance(i,j,k)^2=",
@@ -6055,16 +6098,42 @@ c        write (iout,*) "ii i j",ii,i,j," odleg2",odleg2 ! sum of exps
          write (iout,*) "godl",(godl(k),k=1,constr_homology) ! exponents
          write (iout,*) "ii i j",ii,i,j," odleg2",odleg2 ! sum of exps
 #endif
-         odleg=odleg-dLOG(odleg2/constr_homology)+min_odl
+           if (waga_dist.ge.0.0d0) then
+c
+c          For Gaussian-type Urestr
+c
+              odleg=odleg-dLOG(odleg2/constr_homology)+min_odl
+c
+c          For Lorentzian-type Urestr
+c
+           else
+              odleg=odleg+odleg2/constr_homology
+           endif
+c
 c        write (iout,*) "odleg",odleg ! sum of -ln-s
 c Gradient
-         sum_godl=odleg2
+c
+c          For Gaussian-type Urestr
+c
+         if (waga_dist.ge.0.0d0) sum_godl=odleg2
          sum_sgodl=0.0d0
          do k=1,constr_homology
 c            godl=dexp(((-(distance(i,j,k)**2)/(2*(sigma_odl(i,j,k))**2))
 c     &           *waga_dist)+min_odl
 c          sgodl=-godl(k)*distance(k)*sigma_odl(k,ii)*waga_dist
+c
+         if(.not.l_homo(k,ii)) cycle
+         if (waga_dist.ge.0.0d0) then
+c          For Gaussian-type Urestr
+c
            sgodl=-godl(k)*distance(k)*sigma_odl(k,ii) ! waga_dist rmvd
+c
+c          For Lorentzian-type Urestr
+c
+         else
+           sgodl=-2*sigma_odlir(k,ii)*(distance(k)/(distance(k)**2+
+     &           sigma_odlir(k,ii)**2)**2)
+         endif
            sum_sgodl=sum_sgodl+sgodl
 
 c            sgodl2=sgodl2+sgodl
@@ -6072,12 +6141,22 @@ c      write(iout,*) i, j, k, distance(i,j,k), "W GRADIENCIE1"
 c      write(iout,*) "constr_homology=",constr_homology
 c      write(iout,*) i, j, k, "TEST K"
          enddo
-
-       if (homol_nset.gt.1)then
-         grad_odl3=waga_dist1(iset)*sum_sgodl/(sum_godl*dij)
-       else
-         grad_odl3=waga_dist*sum_sgodl/(sum_godl*dij)
-       endif
+         if (waga_dist.ge.0.0d0) then
+c
+c          For Gaussian-type Urestr
+c
+            grad_odl3=waga_homology(iset)*waga_dist
+     &                *sum_sgodl/(sum_godl*dij)
+c
+c          For Lorentzian-type Urestr
+c
+         else
+c Original grad expr modified by analogy w Gaussian-type Urestr grad
+c           grad_odl3=-waga_homology(iset)*waga_dist*sum_sgodl
+            grad_odl3=-waga_homology(iset)*waga_dist*
+     &                sum_sgodl/(constr_homology*dij)
+         endif
+c
 c        grad_odl3=sum_sgodl/(sum_godl*dij)
 
 
@@ -6175,11 +6254,7 @@ c         sgdih=-gdih(k)*dih_diff(k)*sigma_dih(k,i)*waga_angle
           sum_sgdih=sum_sgdih+sgdih
         enddo
 c       grad_dih3=sum_sgdih/sum_gdih
-        if (homol_nset.gt.1)then
-         grad_dih3=waga_angle1(iset)*sum_sgdih/sum_gdih
-        else
-         grad_dih3=waga_angle*sum_sgdih/sum_gdih
-        endif
+        grad_dih3=waga_homology(iset)*waga_angle*sum_sgdih/sum_gdih
 
 c      write(iout,*)i,k,gdih,sgdih,beta(i+1,i+2,i+3,i+4),grad_dih3
 ccc      write(iout,747) "GRAD_KAT_1", i, nphi, icg, grad_dih3,
 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
+        gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)
+     &      +sum_sgtheta/sum_gtheta*waga_theta
+     &               *waga_homology(iset)
+c        dutheta(i-2)=sum_sgtheta/sum_gtheta*waga_theta
+c     &               *waga_homology(iset)
 c       dutheta(i)=sum_sgtheta/sum_gtheta
 c
 c       Uconst_back=Uconst_back+waga_theta*utheta(i) ! waga_theta added as weight
@@ -6365,7 +6441,7 @@ c         sum_sguscdiff=sum_sguscdiff+sum_guscdiff
 c
 c
 c        New implementation
-         sum_guscdiff = waga_d*sum_guscdiff
+         sum_guscdiff = waga_homology(iset)*waga_d*sum_guscdiff
          do jik=1,3
             duscdiff(jik,i-1)=duscdiff(jik,i-1)+
      &      sum_guscdiff*(dXX_C1tab(jik,i)*dxx+
 c Addition of energy of theta angle and SC local geom over constr_homologs ref strs
 c
 c     ehomology_constr=odleg+kat
-      if (homol_nset.gt.1)then
-       ehomology_constr=waga_dist1(iset)*odleg+waga_angle1(iset)*kat
-     &              +waga_theta*Eval
-     &              +waga_d*Erot     
+c
+c     For Lorentzian-type Urestr
+c
+
+      if (waga_dist.ge.0.0d0) then
+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     write (iout,*) "ehomology_constr=",ehomology_constr
       else
-       ehomology_constr=waga_dist*odleg+waga_angle*kat+waga_theta*Eval
-     &              +waga_d*Erot
+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     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
 c
 c FP 01/15 end