X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD%2Fenergy_p_new_barrier.F;h=168a00020bf741025afa0eab7ece4da0eea8c416;hb=bee3b66ba871f85ef86d5328a081dcf2e974f6f5;hp=5be1ec503517f6370f37803f3614d3ebeb77cc71;hpb=e8aa82ee9b741d0ccebec2934f21bd1f204cdee7;p=unres.git

diff --git a/source/unres/src_MD/energy_p_new_barrier.F b/source/unres/src_MD/energy_p_new_barrier.F
index 5be1ec5..168a000 100644
--- a/source/unres/src_MD/energy_p_new_barrier.F
+++ b/source/unres/src_MD/energy_p_new_barrier.F
@@ -263,7 +263,7 @@ cd    print *,'nterm=',nterm
        edihcnstr=0
       endif
 
-      if (constr_homology.ge.1) then
+      if (constr_homology.ge.1.and.waga_homology(iset).ne.0d0) then
         call e_modeller(ehomology_constr)
 c        print *,'iset=',iset,'me=',me,ehomology_constr,
 c     &  'Processor',fg_rank,' CG group',kolor,
@@ -822,7 +822,7 @@ c      enddo
 #endif
         enddo
       enddo 
-      if (constr_homology.gt.0) then
+      if (constr_homology.gt.0.and.waga_homology(iset).ne.0d0) then
         do i=1,nct
           do j=1,3
             gradc(j,i,icg)=gradc(j,i,icg)+duscdiff(j,i)
@@ -1781,9 +1781,10 @@ c     &        " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
      &        evdwij
             endif
 
-            if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') 
-     &                        'evdw',i,j,evdwij
-
+            if (energy_dec) then
+              write (iout,'(a6,2i5,0pf7.3)') 'evdw',i,j,evdwij
+              call flush(iout)
+            endif
 C Calculate gradient components.
             e1=e1*eps1*eps2rt**2*eps3rt**2
             fac=-expon*(e1+evdwij)*rij_shift
@@ -4584,6 +4585,8 @@ c
       do i=ibondp_start,ibondp_end
         diff = vbld(i)-vbldp0
 c        write (iout,*) i,vbld(i),vbldp0,diff,AKP*diff*diff
+        if (energy_dec)    write (iout,'(a7,i5,4f7.3)') 
+     &     "estr bb",i,vbld(i),vbldp0,diff,AKP*diff*diff
         estr=estr+diff*diff
         do j=1,3
           gradb(j,i-1)=AKP*diff*dc(j,i-1)/vbld(i)
@@ -4602,6 +4605,12 @@ c
             diff=vbld(i+nres)-vbldsc0(1,iti)
 c            write (iout,*) i,iti,vbld(i+nres),vbldsc0(1,iti),diff,
 c     &      AKSC(1,iti),AKSC(1,iti)*diff*diff
+            if (energy_dec)  then
+              write (iout,*) 
+     &         "estr sc",i,iti,vbld(i+nres),vbldsc0(1,iti),diff,
+     &         AKSC(1,iti),AKSC(1,iti)*diff*diff
+              call flush(iout)
+            endif
             estr=estr+0.5d0*AKSC(1,iti)*diff*diff
             do j=1,3
               gradbx(j,i)=AKSC(1,iti)*diff*dc(j,i+nres)/vbld(i+nres)
@@ -4894,7 +4903,7 @@ c     &     " ithet_end",ithet_end
           sinkt(k)=dsin(k*theti2)
         enddo
 C        if (i.gt.3) then
-         if (i.gt.3 .and. itype(i-3).ne.ntyp1) then
+        if (i.gt.3 .and. itype(imax0(i-3,1)).ne.ntyp1) then
 #ifdef OSF
           phii=phi(i)
           if (phii.ne.phii) phii=150.0
@@ -5041,6 +5050,8 @@ c        lprn1=.true.
      &   phii1*rad2deg,ethetai
 c        lprn1=.false.
         etheta=etheta+ethetai
+        if (energy_dec) write (iout,'(a6,i5,0pf7.3)')
+     &      'ebend',i,ethetai
         if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*dephii
         if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*dephii1
         gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wang*dethetai
@@ -6018,7 +6029,7 @@ c
 c
 
 
-      do i=1,19
+      do i=1,max_template
         distancek(i)=9999999.9
       enddo
 
@@ -6038,6 +6049,8 @@ 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
@@ -6070,6 +6083,7 @@ 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
@@ -6119,6 +6133,7 @@ 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
@@ -6578,12 +6593,14 @@ C 6/23/01 Compute double torsional energy
       include 'COMMON.IOUNITS'
       include 'COMMON.FFIELD'
       include 'COMMON.TORCNSTR'
+      include 'COMMON.CONTROL'
       logical lprn
 C Set lprn=.true. for debugging
       lprn=.false.
 c     lprn=.true.
       etors_d=0.0D0
       do i=iphid_start,iphid_end
+        etors_d_ii=0.0D0
         itori=itortyp(itype(i-2))
         itori1=itortyp(itype(i-1))
         itori2=itortyp(itype(i))
@@ -6602,6 +6619,8 @@ c     lprn=.true.
           sinphi2=dsin(j*phii1)
           etors_d=etors_d+v1cij*cosphi1+v1sij*sinphi1+
      &     v2cij*cosphi2+v2sij*sinphi2
+          if (energy_dec) etors_d_ii=etors_d_ii+
+     &     v1cij*cosphi1+v1sij*sinphi1+v2cij*cosphi2+v2sij*sinphi2
           gloci1=gloci1+j*(v1sij*cosphi1-v1cij*sinphi1)
           gloci2=gloci2+j*(v2sij*cosphi2-v2cij*sinphi2)
         enddo
@@ -6617,12 +6636,17 @@ c     lprn=.true.
             sinphi1m2=dsin(l*phii-(k-l)*phii1)
             etors_d=etors_d+v1cdij*cosphi1p2+v2cdij*cosphi1m2+
      &        v1sdij*sinphi1p2+v2sdij*sinphi1m2
+            if (energy_dec) etors_d_ii=etors_d_ii+
+     &        v1cdij*cosphi1p2+v2cdij*cosphi1m2+
+     &        v1sdij*sinphi1p2+v2sdij*sinphi1m2
             gloci1=gloci1+l*(v1sdij*cosphi1p2+v2sdij*cosphi1m2
      &        -v1cdij*sinphi1p2-v2cdij*sinphi1m2)
             gloci2=gloci2+(k-l)*(v1sdij*cosphi1p2-v2sdij*cosphi1m2
      &        -v1cdij*sinphi1p2+v2cdij*sinphi1m2) 
           enddo
         enddo
+        if (energy_dec) write (iout,'(a6,i5,0pf7.3)')
+     &        'etor_d',i,etors_d_ii
         gloc(i-3,icg)=gloc(i-3,icg)+wtor_d*gloci1
         gloc(i-2,icg)=gloc(i-2,icg)+wtor_d*gloci2
 c        write (iout,*) "gloci", gloc(i-3,icg)