correction after cherry-pick and merge in last commit
[unres.git] / source / cluster / wham / src / energy_p_new.F
index 1b69eb3..af14761 100644 (file)
@@ -107,7 +107,7 @@ C
 #ifdef SPLITELE
       etot=wsc*evdw+wscp*evdw2+welec*fact(1)*ees+wvdwpp*evdw1
      & +wang*ebe+wtor*fact(1)*etors+wscloc*escloc
-     & +wstrain*ehpb+nss*ebr+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5
+     & +wstrain*ehpb+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5
      & +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4
      & +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6
      & +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d
@@ -115,7 +115,7 @@ C
 #else
       etot=wsc*evdw+wscp*evdw2+welec*fact(1)*(ees+evdw1)
      & +wang*ebe+wtor*fact(1)*etors+wscloc*escloc
-     & +wstrain*ehpb+nss*ebr+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5
+     & +wstrain*ehpb+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5
      & +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4
      & +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6
      & +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d
@@ -152,6 +152,7 @@ C
       energia(18)=estr
       energia(19)=esccor
       energia(20)=edihcnstr
+cc      if (dyn_ss) call dyn_set_nss
 c detecting NaNQ
       i=0
 #ifdef WINPGI
@@ -218,7 +219,7 @@ cd        write (iout,*) i,g_corr5_loc(i)
      &   +wturn4*fact(3)*gel_loc_turn4(i)
      &   +wturn3*fact(2)*gel_loc_turn3(i)
      &   +wturn6*fact(5)*gel_loc_turn6(i)
-     &   +wel_loc*fact(2)*gel_loc_loc(i)+
+     &   +wel_loc*fact(2)*gel_loc_loc(i)
      &   +wsccor*fact(1)*gsccor_loc(i)
       enddo
       endif
@@ -336,7 +337,7 @@ C
       implicit real*8 (a-h,o-z)
       include 'DIMENSIONS'
       include 'sizesclu.dat'
-      include "DIMENSIONS.COMPAR"
+c      include "DIMENSIONS.COMPAR"
       parameter (accur=1.0d-10)
       include 'COMMON.GEO'
       include 'COMMON.VAR'
@@ -497,7 +498,7 @@ C
       implicit real*8 (a-h,o-z)
       include 'DIMENSIONS'
       include 'sizesclu.dat'
-      include "DIMENSIONS.COMPAR"
+c      include "DIMENSIONS.COMPAR"
       include 'COMMON.GEO'
       include 'COMMON.VAR'
       include 'COMMON.LOCAL'
@@ -586,7 +587,7 @@ C
       implicit real*8 (a-h,o-z)
       include 'DIMENSIONS'
       include 'sizesclu.dat'
-      include "DIMENSIONS.COMPAR"
+c      include "DIMENSIONS.COMPAR"
       include 'COMMON.GEO'
       include 'COMMON.VAR'
       include 'COMMON.LOCAL'
@@ -713,7 +714,7 @@ C
       implicit real*8 (a-h,o-z)
       include 'DIMENSIONS'
       include 'sizesclu.dat'
-      include "DIMENSIONS.COMPAR"
+c      include "DIMENSIONS.COMPAR"
       include 'COMMON.GEO'
       include 'COMMON.VAR'
       include 'COMMON.LOCAL'
@@ -723,6 +724,7 @@ C
       include 'COMMON.INTERACT'
       include 'COMMON.IOUNITS'
       include 'COMMON.CALC'
+      include 'COMMON.SBRIDGE'
       logical lprn
       common /srutu/icall
       integer icant
@@ -748,6 +750,12 @@ C Calculate SC interaction energy.
 C
         do iint=1,nint_gr(i)
           do j=istart(i,iint),iend(i,iint)
+            IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN
+              call dyn_ssbond_ene(i,j,evdwij)
+              evdw=evdw+evdwij
+c              if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)')
+c     &                        'evdw',i,j,evdwij,' ss'
+            ELSE
             ind=ind+1
             itypj=itype(j)
             dscj_inv=vbld_inv(j+nres)
@@ -830,6 +838,7 @@ C Calculate the radial part of the gradient
 C Calculate angular part of the gradient.
             call sc_grad
             endif
+            ENDIF    ! SSBOND
           enddo      ! j
         enddo        ! iint
       enddo          ! i
@@ -844,7 +853,7 @@ C
       implicit real*8 (a-h,o-z)
       include 'DIMENSIONS'
       include 'sizesclu.dat'
-      include "DIMENSIONS.COMPAR"
+c      include "DIMENSIONS.COMPAR"
       include 'COMMON.GEO'
       include 'COMMON.VAR'
       include 'COMMON.LOCAL'
@@ -854,6 +863,7 @@ C
       include 'COMMON.INTERACT'
       include 'COMMON.IOUNITS'
       include 'COMMON.CALC'
+      include 'COMMON.SBRIDGE'
       common /srutu/ icall
       logical lprn
       integer icant
@@ -879,6 +889,13 @@ C Calculate SC interaction energy.
 C
         do iint=1,nint_gr(i)
           do j=istart(i,iint),iend(i,iint)
+C in case of diagnostics    write (iout,*) "TU SZUKAJ",i,j,dyn_ss_mask(i),dyn_ss_mask(j)
+            IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN
+              call dyn_ssbond_ene(i,j,evdwij)
+              evdw=evdw+evdwij
+c              if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)')
+c     &                        'evdw',i,j,evdwij,' ss'
+            ELSE
             ind=ind+1
             itypj=itype(j)
             dscj_inv=vbld_inv(j+nres)
@@ -961,6 +978,7 @@ C Calculate the radial part of the gradient
 C Calculate angular part of the gradient.
             call sc_grad
             endif
+            ENDIF    ! dyn_ss
           enddo      ! j
         enddo        ! iint
       enddo          ! i
@@ -2794,16 +2812,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 'sizesclu.dat'
       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
@@ -2818,43 +2836,88 @@ 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 (.not.dyn_ss .and. i.le.nss) then
+C 15/02/13 CC dynamic SSbond - additional check
         if (ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then
           call ssbond_ene(iii,jjj,eij)
           ehpb=ehpb+2*eij
+cd          write (iout,*) "eij",eij
+        endif
+        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
@@ -2915,7 +2978,7 @@ C
       deltat12=om2-om1+2.0d0
       cosphi=om12-om1*om2
       eij=akcm*deltad*deltad+akth*(deltat1*deltat1+deltat2*deltat2)
-     &  +akct*deltad*deltat12
+     &  +akct*deltad*deltat12+ebr
      &  +v1ss*cosphi+v2ss*cosphi*cosphi+v3ss*cosphi*cosphi*cosphi
 c      write(iout,*) i,j,"rij",rij,"d0cm",d0cm," akcm",akcm," akth",akth,
 c     &  " akct",akct," deltad",deltad," deltat",deltat1,deltat2,
@@ -4207,7 +4270,7 @@ c       write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg)
       do i=1,ndih_constr
         itori=idih_constr(i)
         phii=phi(itori)
-        difi=phii-phi0(i)
+        difi=pinorm(phii-phi0(i))
         if (difi.gt.drange(i)) then
           difi=difi-drange(i)
           edihcnstr=edihcnstr+0.25d0*ftors*difi**4
@@ -4217,10 +4280,10 @@ c       write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg)
           edihcnstr=edihcnstr+0.25d0*ftors*difi**4
           gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
         endif
-!        write (iout,'(2i5,2f8.3,2e14.5)') i,itori,rad2deg*phii,
-!     &    rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg)
+c        write (iout,'(2i5,2f8.3,2e14.5)') i,itori,rad2deg*phii,
+c     &    rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg)
       enddo
-!      write (iout,*) 'edihcnstr',edihcnstr
+      write (iout,*) 'edihcnstr',edihcnstr
       return
       end
 c------------------------------------------------------------------------------
@@ -4289,24 +4352,26 @@ c       write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg)
       enddo
 ! 6/20/98 - dihedral angle constraints
       edihcnstr=0.0d0
+c      write (iout,*) "Dihedral angle restraint energy"
       do i=1,ndih_constr
-        print *,"i",i
         itori=idih_constr(i)
         phii=phi(itori)
-        difi=phii-phi0(i)
+        difi=pinorm(phii-phi0(i))
+c        write (iout,'(2i5,4f8.3,2e14.5)') i,itori,rad2deg*phii,
+c     &    rad2deg*difi,rad2deg*phi0(i),rad2deg*drange(i)
         if (difi.gt.drange(i)) then
           difi=difi-drange(i)
           edihcnstr=edihcnstr+0.25d0*ftors*difi**4
           gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
+c          write (iout,*) 0.25d0*ftors*difi**4
         else if (difi.lt.-drange(i)) then
           difi=difi+drange(i)
           edihcnstr=edihcnstr+0.25d0*ftors*difi**4
           gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
+c          write (iout,*) 0.25d0*ftors*difi**4
         endif
-!        write (iout,'(2i5,2f8.3,2e14.5)') i,itori,rad2deg*phii,
-!     &    rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg)
       enddo
-!      write (iout,*) 'edihcnstr',edihcnstr
+c      write (iout,*) 'edihcnstr',edihcnstr
       return
       end
 c----------------------------------------------------------------------------
@@ -4409,26 +4474,58 @@ C Set lprn=.true. for debugging
 c      lprn=.true.
 c      write (iout,*) "EBACK_SC_COR",iphi_start,iphi_end,nterm_sccor
       esccor=0.0D0
-      do i=iphi_start,iphi_end
+      do i=itau_start,itau_end
         esccor_ii=0.0D0
-        itori=itype(i-2)
-        itori1=itype(i-1)
+        isccori=isccortyp(itype(i-2))
+        isccori1=isccortyp(itype(i-1))
         phii=phi(i)
+cccc  Added 9 May 2012
+cc Tauangle is torsional engle depending on the value of first digit 
+c(see comment below)
+cc Omicron is flat angle depending on the value of first digit 
+c(see comment below)
+
+
+        do intertyp=1,3 !intertyp
+cc Added 09 May 2012 (Adasko)
+cc  Intertyp means interaction type of backbone mainchain correlation: 
+c   1 = SC...Ca...Ca...Ca
+c   2 = Ca...Ca...Ca...SC
+c   3 = SC...Ca...Ca...SCi
         gloci=0.0D0
-        do j=1,nterm_sccor
-          v1ij=v1sccor(j,itori,itori1)
-          v2ij=v2sccor(j,itori,itori1)
-          cosphi=dcos(j*phii)
-          sinphi=dsin(j*phii)
+        if (((intertyp.eq.3).and.((itype(i-2).eq.10).or.
+     &      (itype(i-1).eq.10).or.(itype(i-2).eq.21).or.
+     &      (itype(i-1).eq.21)))
+     &    .or. ((intertyp.eq.1).and.((itype(i-2).eq.10)
+     &     .or.(itype(i-2).eq.21)))
+     &    .or.((intertyp.eq.2).and.((itype(i-1).eq.10).or.
+     &      (itype(i-1).eq.21)))) cycle
+        if ((intertyp.eq.2).and.(i.eq.4).and.(itype(1).eq.21)) cycle
+        if ((intertyp.eq.1).and.(i.eq.nres).and.(itype(nres).eq.21))
+     & cycle
+        do j=1,nterm_sccor(isccori,isccori1)
+          v1ij=v1sccor(j,intertyp,isccori,isccori1)
+          v2ij=v2sccor(j,intertyp,isccori,isccori1)
+          cosphi=dcos(j*tauangle(intertyp,i))
+          sinphi=dsin(j*tauangle(intertyp,i))
           esccor=esccor+v1ij*cosphi+v2ij*sinphi
-          gloci=gloci+fact*j*(v2ij*cosphi-v1ij*sinphi)
+#ifdef DEBUG
+          esccor_ii=esccor_ii+v1ij*cosphi+v2ij*sinphi
+#endif
+          gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi)
         enddo
+        gloc_sc(intertyp,i-3,icg)=gloc_sc(intertyp,i-3,icg)+wsccor*gloci
+c        write (iout,*) "WTF",intertyp,i,itype(i),v1ij*cosphi+v2ij*sinphi
+c     &gloc_sc(intertyp,i-3,icg)
         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,
-     &  (v1sccor(j,itori,itori1),j=1,6),(v2sccor(j,itori,itori1),j=1,6)
-        gsccor_loc(i-3)=gloci
+     &  (v1sccor(j,intertyp,itori,itori1),j=1,6)
+     & ,(v2sccor(j,intertyp,itori,itori1),j=1,6)
+        gsccor_loc(i-3)=gsccor_loc(i-3)+gloci
+       enddo !intertyp
       enddo
+
       return
       end
 c------------------------------------------------------------------------------
@@ -4796,6 +4893,7 @@ C This subroutine calculates multi-body contributions to hydrogen-bonding
 C Set lprn=.true. for debugging
       lprn=.false.
       eturn6=0.0d0
+      ecorr6=0.0d0
 #ifdef MPL
       n_corr=0
       n_corr1=0
@@ -4972,10 +5070,10 @@ cd                write(2,*)'wcorr6',wcorr6,' wturn6',wturn6
 cd                write(2,*)'ijkl',i,j,i+1,j1 
                 if (wcorr6.gt.0.0d0 .and. (j.ne.i+4 .or. j1.ne.i+3
      &               .or. wturn6.eq.0.0d0))then
-cd                  write (iout,*) '******ecorr6: i,j,i+1,j1',i,j,i+1,j1
-                  ecorr6=ecorr6+eello6(i,j,i+1,j1,jj,kk)
-cd                write (iout,*) 'ecorr',ecorr,' ecorr5=',ecorr5,
-cd     &            'ecorr6=',ecorr6
+c                  write (iout,*) '******ecorr6: i,j,i+1,j1',i,j,i+1,j1
+c                  ecorr6=ecorr6+eello6(i,j,i+1,j1,jj,kk)
+c                write (iout,*) 'ecorr',ecorr,' ecorr5=',ecorr5,
+c     &            'ecorr6=',ecorr6, wcorr6
 cd                write (iout,'(4e15.5)') sred_geom,
 cd     &          dabs(eello4(i,j,i+1,j1,jj,kk)),
 cd     &          dabs(eello5(i,j,i+1,j1,jj,kk)),
@@ -6195,18 +6293,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))
@@ -6298,22 +6396,22 @@ c----------------------------------------------------------------------------
       include 'COMMON.GEO'
       logical swap
       double precision vv(2),pizda(2,2),auxmat(2,2),auxvec(2),
-     & auxvec1(2),auxvec2(1),auxmat1(2,2)
+     & auxvec1(2),auxvec2(2),auxmat1(2,2)
       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, 
@@ -6486,18 +6584,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 
@@ -6604,18 +6702,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