Lorentzian-type distance restraint energy implemented

diff --git a/source/unres/src_MD/COMMON.MD b/source/unres/src_MD/COMMON.MD
index 24fc1b3..bc0014b 100644 (file)
--- a/source/unres/src_MD/COMMON.MD
+++ b/source/unres/src_MD/COMMON.MD
@@ -15,7 +15,8 @@
      & vtot(MAXRES2),Gvec(maxres2,maxres2),Geigen(maxres2)
 
        real*8 odl(max_template,maxdim),sigma_odl(max_template,maxdim),
-     &    dih(max_template,maxres),sigma_dih(max_template,maxres)
+     &    dih(max_template,maxres),sigma_dih(max_template,maxres),
+     &    sigma_odlir(max_template,maxdim)
 c
 c    Specification of new variables used in  subroutine e_modeller
 c    modified by FP (Nov.,2014)
@@ -62,9 +63,9 @@ c
      & lim_odl,lim_dih,ires_homo,jres_homo,link_start_homo,
      & link_end_homo,idihconstr_start_homo,idihconstr_end_homo,
 c
-c    FP (30/10/2014)
+c    FP (30/10/2014,04/03/2015)
 c
-     & xxtpl,yytpl,zztpl,thetatpl,sigma_theta,sigma_d
+     & xxtpl,yytpl,zztpl,thetatpl,sigma_theta,sigma_d,sigma_odlir
 c
       common /qmeas/ qfrag,qpair,qinfrag,qinpair,wfrag,wpair,eq_time,
      & Ucdfrag,Ucdpair,dUdconst,dUdxconst,dqwol,dxqwol,Uconst,

diff --git a/source/unres/src_MD/cinfo.f b/source/unres/src_MD/cinfo.f
index 81455fd..a688eb6 100644 (file)
--- a/source/unres/src_MD/cinfo.f
+++ b/source/unres/src_MD/cinfo.f
@@ -1,10 +1,10 @@
 C DO NOT EDIT THIS FILE - IT HAS BEEN GENERATED BY COMPINFO.C
-C 0 40376 39
+C 0 40376 53
       subroutine cinfo
       include 'COMMON.IOUNITS'
       write(iout,*)'++++ Compile info ++++'
-      write(iout,*)'Version 0.40376 build 39'
-      write(iout,*)'compiled Tue Mar  3 12:57:02 2015'
+      write(iout,*)'Version 0.40376 build 53'
+      write(iout,*)'compiled Fri Mar  6 10:04:56 2015'
       write(iout,*)'compiled by felipe@piasek4'
       write(iout,*)'OS name:    Linux '
       write(iout,*)'OS release: 3.2.0-70-generic '

diff --git a/source/unres/src_MD/energy_p_new_barrier.F b/source/unres/src_MD/energy_p_new_barrier.F
index 15dfccd..d07d135 100644 (file)
--- a/source/unres/src_MD/energy_p_new_barrier.F
+++ b/source/unres/src_MD/energy_p_new_barrier.F
@@ -5851,7 +5851,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
@@ -6021,10 +6021,21 @@ c        write (iout,*) "dij(",i,j,") =",dij
          do k=1,constr_homology
            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 +6051,18 @@ 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 (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 +6076,41 @@ 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 (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,9 +6118,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
-
-         grad_odl3=waga_homology(iset)*waga_dist
-     &            *sum_sgodl/(sum_godl*dij)
+         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)
 
 
@@ -6447,8 +6506,25 @@ c
 c Addition of energy of theta angle and SC local geom over constr_homologs ref strs
 c
 c     ehomology_constr=odleg+kat
-       ehomology_constr=(waga_dist*odleg+waga_angle*kat+waga_theta*Eval
-     &              +waga_d*Erot)*waga_homology(iset)
+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
+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

diff --git a/source/unres/src_MD/readrtns.F b/source/unres/src_MD/readrtns.F
index 636c6a4..b731cf2 100644 (file)
--- a/source/unres/src_MD/readrtns.F
+++ b/source/unres/src_MD/readrtns.F
@@ -2841,7 +2841,7 @@ c    &       sigma_odl_temp(i+nnt-1,j+nnt-1,k)
 c         enddo
 c 1401   continue
 c         close (ientin)
-        if (waga_dist.gt.0.0d0) then
+        if (waga_dist.ne.0.0d0) then
           ii=0
           do i = nnt,nct-2 ! right? without parallel.
             do j=i+2,nct ! right?
@@ -3006,7 +3006,7 @@ c              read (ientin,*) sigma_d(k,i) ! 1st variant
         endif
         close(ientin)
       enddo
-      if (waga_dist.gt.0.0d0) lim_odl=ii
+      if (waga_dist.ne.0.0d0) lim_odl=ii
       if (constr_homology.gt.0) call homology_partition
       if (constr_homology.gt.0) call init_int_table
 cd      write (iout,*) "homology_partition: lim_theta= ",lim_theta,