X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;ds=sidebyside;f=source%2Funres%2Fsrc_CSA_DiL%2Fenergy_p_new_barrier.F;h=c1e8ad3025e2678f4ce802dca2c8faf38ad2f461;hb=f690e8b70bab14132839afebf080d4a28363b226;hp=bb1c8a14cd33f6001b73370f91e6aaf63caa1db7;hpb=112d5a40499d065fa16ade60c2eb3f815f136ca4;p=unres.git

diff --git a/source/unres/src_CSA_DiL/energy_p_new_barrier.F b/source/unres/src_CSA_DiL/energy_p_new_barrier.F
index bb1c8a1..c1e8ad3 100644
--- a/source/unres/src_CSA_DiL/energy_p_new_barrier.F
+++ b/source/unres/src_CSA_DiL/energy_p_new_barrier.F
@@ -4485,7 +4485,19 @@ c     write (*,'(a,i2)') 'EBEND ICG=',icg
       do i=ithet_start,ithet_end
 C Zero the energy function and its derivative at 0 or pi.
         call splinthet(theta(i),0.5d0*delta,ss,ssd)
-        it=iabs(itype(i-1))
+        it=itype(i-1)
+        ichir1=isign(1,itype(i-2))
+        ichir2=isign(1,itype(i))
+        if (itype(i-2).eq.10) ichir1=isign(1,itype(i-1))
+        if (itype(i).eq.10) ichir2=isign(1,itype(i-1))
+        if (itype(i-1).eq.10) then
+         itype1=isign(10,itype(i-2))
+         ichir11=isign(1,itype(i-2))
+         ichir12=isign(1,itype(i-2))
+         itype2=isign(10,itype(i))
+         ichir21=isign(1,itype(i))
+         ichir22=isign(1,itype(i))
+        endif
         if (i.gt.3) then
 #ifdef OSF
 	  phii=phi(i)
@@ -4519,15 +4531,27 @@ C dependent on the adjacent virtual-bond-valence angles (gamma1 & gamma2).
 C In following comments this theta will be referred to as t_c.
         thet_pred_mean=0.0d0
         do k=1,2
-          athetk=athet(k,it)
-          bthetk=bthet(k,it)
+          athetk=athet(k,it,ichir1,ichir2)
+          bthetk=bthet(k,it,ichir1,ichir2)
+        if (it.eq.10) then
+           athetk=athet(k,itype1,ichir11,ichir12)
+           bthetk=bthet(k,itype2,ichir21,ichir22)
+        endif
           thet_pred_mean=thet_pred_mean+athetk*y(k)+bthetk*z(k)
         enddo
         dthett=thet_pred_mean*ssd
         thet_pred_mean=thet_pred_mean*ss+a0thet(it)
 C Derivatives of the "mean" values in gamma1 and gamma2.
-        dthetg1=(-athet(1,it)*y(2)+athet(2,it)*y(1))*ss
-        dthetg2=(-bthet(1,it)*z(2)+bthet(2,it)*z(1))*ss
+        dthetg1=(-athet(1,it,ichir1,ichir2)*y(2)
+     &+athet(2,it,ichir1,ichir2)*y(1))*ss
+        dthetg2=(-bthet(1,it,ichir1,ichir2)*z(2)
+     &          +bthet(2,it,ichir1,ichir2)*z(1))*ss
+        if (it.eq.10) then
+      dthetg1=(-athet(1,itype1,ichir11,ichir12)*y(2)
+     &+athet(2,itype1,ichir11,ichir12)*y(1))*ss
+        dthetg2=(-bthet(1,itype2,ichir21,ichir22)*z(2)
+     &         +bthet(2,itype2,ichir21,ichir22)*z(1))*ss
+        endif
         if (theta(i).gt.pi-delta) then
           call theteng(pi-delta,thet_pred_mean,theta0(it),f0,fprim0,
      &         E_tc0)
@@ -5609,6 +5633,11 @@ c      lprn=.true.
       etors_ii=0.0D0
 	itori=itortyp(itype(i-2))
 	itori1=itortyp(itype(i-1))
+        if (iabs(itype(i)).eq.20) then
+        iblock=2
+        else
+        iblock=1
+        endif
         phii=phi(i)
         gloci=0.0D0
 C Proline-Proline pair is a special case...
@@ -5708,9 +5737,9 @@ c     lprn=.true.
         phii=phi(i)
         gloci=0.0D0
 C Regular cosine and sine terms
-        do j=1,nterm(itori,itori1)
-          v1ij=v1(j,itori,itori1)
-          v2ij=v2(j,itori,itori1)
+        do j=1,nterm(itori,itori1,iblock)
+          v1ij=v1(j,itori,itori1,iblock)
+          v2ij=v2(j,itori,itori1,iblock)
           cosphi=dcos(j*phii)
           sinphi=dsin(j*phii)
           etors=etors+v1ij*cosphi+v2ij*sinphi
@@ -5725,7 +5754,7 @@ C          [v2 cos(phi/2)+v3 sin(phi/2)]^2 + 1
 C
         cosphi=dcos(0.5d0*phii)
         sinphi=dsin(0.5d0*phii)
-        do j=1,nlor(itori,itori1)
+        do j=1,nlor(itori,itori1,iblock)
           vl1ij=vlor1(j,itori,itori1)
           vl2ij=vlor2(j,itori,itori1)
           vl3ij=vlor3(j,itori,itori1)
@@ -5738,13 +5767,14 @@ C
           gloci=gloci+vl1ij*(vl3ij*cosphi-vl2ij*sinphi)*pom
         enddo
 C Subtract the constant term
-        etors=etors-v0(itori,itori1)
+        etors=etors-v0(itori,itori1,iblock)
           if (energy_dec) write (iout,'(a6,i5,0pf7.3)')
-     &         'etor',i,etors_ii-v0(itori,itori1)
+     &         'etor',i,etors_ii-v0(itori,itori1,iblock)
         if (lprn)
      &  write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)')
      &  restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,itori,itori1,
-     &  (v1(j,itori,itori1),j=1,6),(v2(j,itori,itori1),j=1,6)
+     &  (v1(j,itori,itori1,iblock),j=1,6),
+     &  (v2(j,itori,itori1,iblock),j=1,6)
         gloc(i-3,icg)=gloc(i-3,icg)+wtor*gloci
 c       write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg)
       enddo
@@ -5799,7 +5829,7 @@ c     lprn=.true.
         itori1=itortyp(itype(i-1))
         itori2=itortyp(itype(i))
         iblock=1
-        if (iabs(itype(i+1).eq.20)) iblock=2
+        if (iabs(itype(i+1)).eq.20) iblock=2
         phii=phi(i)
         phii1=phi(i+1)
         gloci1=0.0D0