Dzialajaczy dla D aminokwasow whama(mono) psuedosymetria dla GABowego
[unres.git] / source / wham / src / energy_p_new.F
index 06bdd76..fe9f889 100644 (file)
@@ -368,8 +368,8 @@ cd    print *,'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
       evdw=0.0D0
       evdw_t=0.0d0
       do i=iatsc_s,iatsc_e
-        itypi=itype(i)
-        itypi1=itype(i+1)
+        itypi=iabs(itype(i))
+        itypi1=iabs(itype(i+1))
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
@@ -382,7 +382,7 @@ C
 cd        write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
 cd   &                  'iend=',iend(i,iint)
           do j=istart(i,iint),iend(i,iint)
-            itypj=itype(j)
+            itypj=iabs(itype(j))
             xj=c(1,nres+j)-xi
             yj=c(2,nres+j)-yi
             zj=c(3,nres+j)-zi
@@ -539,8 +539,8 @@ c     print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
       evdw=0.0D0
       evdw_t=0.0d0
       do i=iatsc_s,iatsc_e
-        itypi=itype(i)
-        itypi1=itype(i+1)
+        itypi=iabs(itype(i))
+        itypi1=iabs(itype(i+1))
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
@@ -549,7 +549,7 @@ C Calculate SC interaction energy.
 C
         do iint=1,nint_gr(i)
           do j=istart(i,iint),iend(i,iint)
-            itypj=itype(j)
+            itypj=iabs(itype(j))
             xj=c(1,nres+j)-xi
             yj=c(2,nres+j)-yi
             zj=c(3,nres+j)-zi
@@ -650,8 +650,8 @@ c     else
 c     endif
       ind=0
       do i=iatsc_s,iatsc_e
-        itypi=itype(i)
-        itypi1=itype(i+1)
+        itypi=iabs(itype(i))
+        itypi1=iabs(itype(i+1))
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
@@ -665,7 +665,7 @@ C
         do iint=1,nint_gr(i)
           do j=istart(i,iint),iend(i,iint)
             ind=ind+1
-            itypj=itype(j)
+            itypj=iabs(itype(j))
             dscj_inv=vbld_inv(j+nres)
             chi1=chi(itypi,itypj)
             chi2=chi(itypj,itypi)
@@ -786,8 +786,8 @@ c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
 c      if (icall.gt.0) lprn=.true.
       ind=0
       do i=iatsc_s,iatsc_e
-        itypi=itype(i)
-        itypi1=itype(i+1)
+        itypi=iabs(itype(i))
+        itypi1=iabs(itype(i+1))
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
@@ -801,7 +801,7 @@ C
         do iint=1,nint_gr(i)
           do j=istart(i,iint),iend(i,iint)
             ind=ind+1
-            itypj=itype(j)
+            itypj=iabs(itype(j))
             dscj_inv=vbld_inv(j+nres)
             sig0ij=sigma(itypi,itypj)
             chi1=chi(itypi,itypj)
@@ -931,8 +931,8 @@ c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
 c      if (icall.gt.0) lprn=.true.
       ind=0
       do i=iatsc_s,iatsc_e
-        itypi=itype(i)
-        itypi1=itype(i+1)
+        itypi=iabs(itype(i))
+        itypi1=iabs(itype(i+1))
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
@@ -946,7 +946,7 @@ C
         do iint=1,nint_gr(i)
           do j=istart(i,iint),iend(i,iint)
             ind=ind+1
-            itypj=itype(j)
+            itypj=iabs(itype(j))
             dscj_inv=vbld_inv(j+nres)
             sig0ij=sigma(itypi,itypj)
             r0ij=r0(itypi,itypj)
@@ -2785,7 +2785,7 @@ c     &   " iscp",(iscpstart(i,j),iscpend(i,j),j=1,nscp_gr(i))
         do iint=1,nscp_gr(i)
 
         do j=iscpstart(i,iint),iscpend(i,iint)
-          itypj=itype(j)
+          itypj=iabs(itype(j))
 C Uncomment following three lines for SC-p interactions
 c         xj=c(1,nres+j)-xi
 c         yj=c(2,nres+j)-yi
@@ -2869,16 +2869,16 @@ C Evaluate bridge-strain energy and its gradient in virtual-bond and SC vectors.
 C
       implicit real*8 (a-h,o-z)
       include 'DIMENSIONS'
-      include 'DIMENSIONS.ZSCOPT'
       include 'COMMON.SBRIDGE'
       include 'COMMON.CHAIN'
       include 'COMMON.DERIV'
       include 'COMMON.VAR'
       include 'COMMON.INTERACT'
+      include 'COMMON.IOUNITS'
       dimension ggg(3)
       ehpb=0.0D0
-cd    print *,'edis: nhpb=',nhpb,' fbr=',fbr
-cd    print *,'link_start=',link_start,' link_end=',link_end
+cd      write(iout,*)'edis: nhpb=',nhpb,' fbr=',fbr
+cd      write(iout,*)'link_start=',link_start,' link_end=',link_end
       if (link_end.eq.0) return
       do i=link_start,link_end
 C If ihpb(i) and jhpb(i) > NRES, this is a SC-SC distance, otherwise a
@@ -2893,43 +2893,86 @@ C iii and jjj point to the residues for which the distance is assigned.
           iii=ii
           jjj=jj
         endif
+c        write (iout,*) "i",i," ii",ii," iii",iii," jj",jj," jjj",jjj,
+c     &    dhpb(i),dhpb1(i),forcon(i)
 C 24/11/03 AL: SS bridges handled separately because of introducing a specific
 C    distance and angle dependent SS bond potential.
-        if (ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then
+        if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and.
+     & iabs(itype(jjj)).eq.1) then
           call ssbond_ene(iii,jjj,eij)
           ehpb=ehpb+2*eij
+cd          write (iout,*) "eij",eij
+        else if (ii.gt.nres .and. jj.gt.nres) then
+c Restraints from contact prediction
+          dd=dist(ii,jj)
+          if (dhpb1(i).gt.0.0d0) then
+            ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+            fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd
+c            write (iout,*) "beta nmr",
+c     &        dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+          else
+            dd=dist(ii,jj)
+            rdis=dd-dhpb(i)
+C Get the force constant corresponding to this distance.
+            waga=forcon(i)
+C Calculate the contribution to energy.
+            ehpb=ehpb+waga*rdis*rdis
+c            write (iout,*) "beta reg",dd,waga*rdis*rdis
+C
+C Evaluate gradient.
+C
+            fac=waga*rdis/dd
+          endif  
+          do j=1,3
+            ggg(j)=fac*(c(j,jj)-c(j,ii))
+          enddo
+          do j=1,3
+            ghpbx(j,iii)=ghpbx(j,iii)-ggg(j)
+            ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j)
+          enddo
+          do k=1,3
+            ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k)
+            ghpbc(k,iii)=ghpbc(k,iii)-ggg(k)
+          enddo
         else
 C Calculate the distance between the two points and its difference from the
 C target distance.
-        dd=dist(ii,jj)
-        rdis=dd-dhpb(i)
+          dd=dist(ii,jj)
+          if (dhpb1(i).gt.0.0d0) then
+            ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+            fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd
+c            write (iout,*) "alph nmr",
+c     &        dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+          else
+            rdis=dd-dhpb(i)
 C Get the force constant corresponding to this distance.
-        waga=forcon(i)
+            waga=forcon(i)
 C Calculate the contribution to energy.
-        ehpb=ehpb+waga*rdis*rdis
+            ehpb=ehpb+waga*rdis*rdis
+c            write (iout,*) "alpha reg",dd,waga*rdis*rdis
 C
 C Evaluate gradient.
 C
-        fac=waga*rdis/dd
+            fac=waga*rdis/dd
+          endif
 cd      print *,'i=',i,' ii=',ii,' jj=',jj,' dhpb=',dhpb(i),' dd=',dd,
 cd   &   ' waga=',waga,' fac=',fac
-        do j=1,3
-          ggg(j)=fac*(c(j,jj)-c(j,ii))
-        enddo
+            do j=1,3
+              ggg(j)=fac*(c(j,jj)-c(j,ii))
+            enddo
 cd      print '(i3,3(1pe14.5))',i,(ggg(j),j=1,3)
 C If this is a SC-SC distance, we need to calculate the contributions to the
 C Cartesian gradient in the SC vectors (ghpbx).
-        if (iii.lt.ii) then
+          if (iii.lt.ii) then
           do j=1,3
             ghpbx(j,iii)=ghpbx(j,iii)-ggg(j)
             ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j)
           enddo
-        endif
-        do j=iii,jjj-1
+          endif
           do k=1,3
-            ghpbc(k,j)=ghpbc(k,j)+ggg(k)
+            ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k)
+            ghpbc(k,iii)=ghpbc(k,iii)-ggg(k)
           enddo
-        enddo
         endif
       enddo
       ehpb=0.5D0*ehpb
@@ -2955,7 +2998,7 @@ C
       include 'COMMON.VAR'
       include 'COMMON.IOUNITS'
       double precision erij(3),dcosom1(3),dcosom2(3),gg(3)
-      itypi=itype(i)
+      itypi=iabs(itype(i))
       xi=c(1,nres+i)
       yi=c(2,nres+i)
       zi=c(3,nres+i)
@@ -2963,7 +3006,7 @@ C
       dyi=dc_norm(2,nres+i)
       dzi=dc_norm(3,nres+i)
       dsci_inv=dsc_inv(itypi)
-      itypj=itype(j)
+      itypj=iabs(itype(j))
       dscj_inv=dsc_inv(itypj)
       xj=c(1,nres+j)-xi
       yj=c(2,nres+j)-yi
@@ -3053,7 +3096,7 @@ c
 c 09/18/07 AL: multimodal bond potential based on AM1 CA-SC PMF's included
 c
       do i=nnt,nct
-        iti=itype(i)
+        iti=iabs(itype(i))
         if (iti.ne.10) then
           nbi=nbondterm(iti)
           if (nbi.eq.1) then
@@ -3133,6 +3176,18 @@ c      write (iout,*) 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=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
 c        if (i.gt.ithet_start .and. 
 c     &     (itel(i-1).eq.0 .or. itel(i-2).eq.0)) goto 1215
 c        if (i.gt.3 .and. (i.le.4 .or. itel(i-3).ne.0)) then
@@ -3188,8 +3243,12 @@ 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
 c        write (iout,*) "thet_pred_mean",thet_pred_mean
@@ -3197,8 +3256,16 @@ c        write (iout,*) "thet_pred_mean",thet_pred_mean
         thet_pred_mean=thet_pred_mean*ss+a0thet(it)
 c        write (iout,*) "thet_pred_mean",thet_pred_mean
 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)
@@ -3373,7 +3440,7 @@ c      write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1)
         dephii=0.0d0
         dephii1=0.0d0
         theti2=0.5d0*theta(i)
-        ityp2=ithetyp(itype(i-1))
+        ityp2=ithetyp(iabs(itype(i-1)))
         do k=1,nntheterm
           coskt(k)=dcos(k*theti2)
           sinkt(k)=dsin(k*theti2)
@@ -3385,7 +3452,7 @@ c      write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1)
 #else
           phii=phi(i)
 #endif
-          ityp1=ithetyp(itype(i-2))
+          ityp1=ithetyp(iabs(itype(i-2)))
           do k=1,nsingle
             cosph1(k)=dcos(k*phii)
             sinph1(k)=dsin(k*phii)
@@ -3406,7 +3473,7 @@ c      write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1)
 #else
           phii1=phi(i+1)
 #endif
-          ityp3=ithetyp(itype(i))
+          ityp3=ithetyp(iabs(itype(i)))
           do k=1,nsingle
             cosph2(k)=dcos(k*phii1)
             sinph2(k)=dsin(k*phii1)
@@ -3560,7 +3627,7 @@ c     write (iout,'(a)') 'ESC'
       do i=loc_start,loc_end
         it=itype(i)
         if (it.eq.10) goto 1
-        nlobit=nlob(it)
+        nlobit=nlob(iabs(it))
 c       print *,'i=',i,' it=',it,' nlobit=',nlobit
 c       write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad
         theti=theta(i+1)-pipol
@@ -3715,7 +3782,7 @@ C Compute the contribution to SC energy and derivatives
         do iii=-1,1
 
           do j=1,nlobit
-            expfac=dexp(bsc(j,it)-0.5D0*contr(j,iii)+emin)
+            expfac=dexp(bsc(j,iabs(it))-0.5D0*contr(j,iii)+emin)
 cd          print *,'j=',j,' expfac=',expfac
             escloc_i=escloc_i+expfac
             do k=1,3
@@ -3796,7 +3863,7 @@ C Compute the contribution to SC energy and derivatives
 
       dersc12=0.0d0
       do j=1,nlobit
-        expfac=dexp(bsc(j,it)-0.5D0*contr(j)+emin)
+        expfac=dexp(bsc(j,iabs(it))-0.5D0*contr(j)+emin)
         escloc_i=escloc_i+expfac
         do k=1,2
           dersc(k)=dersc(k)+Ax(k,j)*expfac
@@ -3859,7 +3926,7 @@ C
         cosfac=dsqrt(cosfac2)
         sinfac2=0.5d0/(1.0d0-costtab(i+1))
         sinfac=dsqrt(sinfac2)
-        it=itype(i)
+        it=iabs(itype(i))
         if (it.eq.10) goto 1
 c
 C  Compute the axes of tghe local cartesian coordinates system; store in
@@ -3899,7 +3966,7 @@ c
         do j = 1,3
           xx = xx + x_prime(j)*dc_norm(j,i+nres)
           yy = yy + y_prime(j)*dc_norm(j,i+nres)
-          zz = zz + z_prime(j)*dc_norm(j,i+nres)
+          zz = zz + dsign(1.0,itype(i))*z_prime(j)*dc_norm(j,i+nres)
         enddo
 
         xxtab(i)=xx
@@ -3909,7 +3976,7 @@ C
 C Compute the energy of the ith side cbain
 C
 c        write (2,*) "xx",xx," yy",yy," zz",zz
-        it=itype(i)
+        it=iabs(itype(i))
         do j = 1,65
           x(j) = sc_parmin(j,it) 
         enddo
@@ -3917,7 +3984,7 @@ c        write (2,*) "xx",xx," yy",yy," zz",zz
 Cc diagnostics - remove later
         xx1 = dcos(alph(2))
         yy1 = dsin(alph(2))*dcos(omeg(2))
-        zz1 = -dsin(alph(2))*dsin(omeg(2))
+        zz1 = -dsign(1.0,itype(i))*dsin(alph(2))*dsin(omeg(2))
         write(2,'(3f8.1,3f9.3,1x,3f9.3)') 
      &    alph(2)*rad2deg,omeg(2)*rad2deg,theta(3)*rad2deg,xx,yy,zz,
      &    xx1,yy1,zz1
@@ -4325,14 +4392,19 @@ c      lprn=.true.
       etors=0.0D0
       do i=iphi_start,iphi_end
         if (itel(i-2).eq.0 .or. itel(i-1).eq.0) goto 1215
+         if (iabs(itype(i)).eq.20) then
+         iblock=2
+         else
+         iblock=1
+         endif
         itori=itortyp(itype(i-2))
         itori1=itortyp(itype(i-1))
         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
@@ -4345,7 +4417,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)
           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 (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,1),j=1,6),(v2(j,itori,itori1,1),j=1,6)
         gloc(i-3,icg)=gloc(i-3,icg)+wtor*fact*gloci
 c       write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg)
  1215   continue
@@ -4421,16 +4493,26 @@ c     lprn=.true.
         itori=itortyp(itype(i-2))
         itori1=itortyp(itype(i-1))
         itori2=itortyp(itype(i))
+c        iblock=1
+c        if (iabs(itype(i+1)).eq.20) iblock=2
         phii=phi(i)
         phii1=phi(i+1)
         gloci1=0.0D0
         gloci2=0.0D0
+        iblock=1
+        if (iabs(itype(i+1)).eq.20) iblock=2
 C Regular cosine and sine terms
-        do j=1,ntermd_1(itori,itori1,itori2)
-          v1cij=v1c(1,j,itori,itori1,itori2)
-          v1sij=v1s(1,j,itori,itori1,itori2)
-          v2cij=v1c(2,j,itori,itori1,itori2)
-          v2sij=v1s(2,j,itori,itori1,itori2)
+c c       do j=1,ntermd_1(itori,itori1,itori2,iblock)
+c          v1cij=v1c(1,j,itori,itori1,itori2,iblock)
+c          v1sij=v1s(1,j,itori,itori1,itori2,iblock)
+c          v2cij=v1c(2,j,itori,itori1,itori2,iblock)
+c          v2sij=v1s(2,j,itori,itori1,itori2,iblock)
+       do j=1,ntermd_1(itori,itori1,itori2,iblock)
+          v1cij=v1c(1,j,itori,itori1,itori2,iblock)
+          v1sij=v1s(1,j,itori,itori1,itori2,iblock)
+          v2cij=v1c(2,j,itori,itori1,itori2,iblock)
+          v2sij=v1s(2,j,itori,itori1,itori2,iblock)
+
           cosphi1=dcos(j*phii)
           sinphi1=dsin(j*phii)
           cosphi2=dcos(j*phii1)
@@ -4440,12 +4522,12 @@ C Regular cosine and sine terms
           gloci1=gloci1+j*(v1sij*cosphi1-v1cij*sinphi1)
           gloci2=gloci2+j*(v2sij*cosphi2-v2cij*sinphi2)
         enddo
-        do k=2,ntermd_2(itori,itori1,itori2)
+        do k=2,ntermd_2(itori,itori1,itori2,iblock)
           do l=1,k-1
-            v1cdij = v2c(k,l,itori,itori1,itori2)
-            v2cdij = v2c(l,k,itori,itori1,itori2)
-            v1sdij = v2s(k,l,itori,itori1,itori2)
-            v2sdij = v2s(l,k,itori,itori1,itori2)
+            v1cdij = v2c(k,l,itori,itori1,itori2,iblock)
+            v2cdij = v2c(l,k,itori,itori1,itori2,iblock)
+            v1sdij = v2s(k,l,itori,itori1,itori2,iblock)
+            v2sdij = v2s(l,k,itori,itori1,itori2,iblock)
             cosphi1p2=dcos(l*phii+(k-l)*phii1)
             cosphi1m2=dcos(l*phii-(k-l)*phii1)
             sinphi1p2=dsin(l*phii+(k-l)*phii1)
@@ -4496,8 +4578,9 @@ c      write (iout,*) "EBACK_SC_COR",iphi_start,iphi_end,nterm_sccor
       esccor=0.0D0
       do i=itau_start,itau_end
         esccor_ii=0.0D0
-        isccori=isccortyp(itype(i-2))
-        isccori1=isccortyp(itype(i-1))
+CC TU NIE MOZE BYC IABS DO ZMIANY JAK BEDA GOTOWE POTENCJALY AM1
+        isccori=isccortyp(iabs(itype(i-2)))
+        isccori1=isccortyp(iabs(itype(i-1)))
         phii=phi(i)
 cccc  Added 9 May 2012
 cc Tauangle is torsional engle depending on the value of first digit 
@@ -6311,18 +6394,18 @@ c--------------------------------------------------------------------------
       logical lprn
       common /kutas/ lprn
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
-C                                              
-C      Parallel       Antiparallel
-C                                             
-C          o             o         
-C         /l\           /j\       
-C        /   \         /   \      
-C       /| o |         | o |\     
-C     \ j|/k\|  /   \  |/k\|l /   
-C      \ /   \ /     \ /   \ /    
-C       o     o       o     o                
-C       i             i                     
-C
+C                                                                              C
+C      Parallel       Antiparallel                                             C
+C                                                                              C
+C          o             o                                                     C
+C         /l\           /j\                                                    C 
+C        /   \         /   \                                                   C
+C       /| o |         | o |\                                                  C
+C     \ j|/k\|  /   \  |/k\|l /                                                C
+C      \ /   \ /     \ /   \ /                                                 C
+C       o     o       o     o                                                  C
+C       i             i                                                        C
+C                                                                              C
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
       itk=itortyp(itype(k))
       s1= scalar2(AEAb1(1,2,imat),CUgb2(1,i))
@@ -6418,18 +6501,18 @@ c----------------------------------------------------------------------------
       logical lprn
       common /kutas/ lprn
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
-C                                              
-C      Parallel       Antiparallel
-C                                             
-C          o             o         
-C     \   /l\           /j\   /   
-C      \ /   \         /   \ /    
-C       o| o |         | o |o     
-C     \ j|/k\|      \  |/k\|l     
-C      \ /   \       \ /   \      
-C       o             o                      
-C       i             i                     
-C
+C                                                                              C 
+C      Parallel       Antiparallel                                             C
+C                                                                              C
+C          o             o                                                     C
+C     \   /l\           /j\   /                                                C
+C      \ /   \         /   \ /                                                 C
+C       o| o |         | o |o                                                  C
+C     \ j|/k\|      \  |/k\|l                                                  C
+C      \ /   \       \ /   \                                                   C
+C       o             o                                                        C
+C       i             i                                                        C
+C                                                                              C
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 cd      write (2,*) 'eello6_graph2: i,',i,' j',j,' k',k,' l',l
 C AL 7/4/01 s1 would occur in the sixth-order moment, 
@@ -6602,18 +6685,18 @@ c----------------------------------------------------------------------------
       double precision vv(2),pizda(2,2),auxmat(2,2),auxvec(2)
       logical swap
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
-C                                              
-C      Parallel       Antiparallel
-C                                             
-C          o             o         
-C         /l\   /   \   /j\       
-C        /   \ /     \ /   \      
-C       /| o |o       o| o |\     
-C       j|/k\|  /      |/k\|l /   
-C        /   \ /       /   \ /    
-C       /     o       /     o                
-C       i             i                     
-C
+C                                                                              C
+C      Parallel       Antiparallel                                             C
+C                                                                              C
+C          o             o                                                     C
+C         /l\   /   \   /j\                                                    C
+C        /   \ /     \ /   \                                                   C
+C       /| o |o       o| o |\                                                  C
+C       j|/k\|  /      |/k\|l /                                                C
+C        /   \ /       /   \ /                                                 C
+C       /     o       /     o                                                  C
+C       i             i                                                        C
+C                                                                              C
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 C
 C 4/7/01 AL Component s1 was removed, because it pertains to the respective 
@@ -6720,18 +6803,18 @@ c----------------------------------------------------------------------------
      & auxvec1(2),auxmat1(2,2)
       logical swap
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
-C                                              
-C      Parallel       Antiparallel
-C                                             
-C          o             o         
-C         /l\   /   \   /j\       
-C        /   \ /     \ /   \      
-C       /| o |o       o| o |\     
-C     \ j|/k\|      \  |/k\|l     
-C      \ /   \       \ /   \      
-C       o     \       o     \                
-C       i             i                     
-C
+C                                                                              C
+C      Parallel       Antiparallel                                             C
+C                                                                              C
+C          o             o                                                     C 
+C         /l\   /   \   /j\                                                    C
+C        /   \ /     \ /   \                                                   C
+C       /| o |o       o| o |\                                                  C
+C     \ j|/k\|      \  |/k\|l                                                  C
+C      \ /   \       \ /   \                                                   C
+C       o     \       o     \                                                  C
+C       i             i                                                        C
+C                                                                              C
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 C
 C 4/7/01 AL Component s1 was removed, because it pertains to the respective