Merge branch 'bartek2' into bartek
authorBartlomiej Zaborowski <bartek.zaborowski@chem.univ.gda.pl>
Tue, 2 Oct 2012 15:51:07 +0000 (17:51 +0200)
committerBartlomiej Zaborowski <bartek.zaborowski@chem.univ.gda.pl>
Tue, 2 Oct 2012 15:51:07 +0000 (17:51 +0200)
Conflicts:
source/unres/src_CSA_DiL/energy_p_new_barrier.F
source/unres/src_CSA_DiL/parmread.F

1  2 
source/unres/src_CSA_DiL/COMMON.LOCAL
source/unres/src_CSA_DiL/energy_p_new_barrier.F
source/unres/src_CSA_DiL/gen_rand_conf.F
source/unres/src_CSA_DiL/initialize_p.F
source/unres/src_CSA_DiL/parmread.F
source/unres/src_CSA_DiL/readrtns_csa.F

@@@ -2,13 -2,13 +2,14 @@@
       &  sigc0,dsc,dsc_inv,bsc,censc,gaussc,dsc0
        integer nlob
  C Parameters of the virtual-bond-angle probability distribution
-       common /thetas/ a0thet(ntyp),athet(2,ntyp),bthet(2,ntyp),
-      &  polthet(0:3,ntyp),gthet(3,ntyp),theta0(ntyp),sig0(ntyp),
-      &  sigc0(ntyp)
+       common /thetas/ a0thet(-ntyp:ntyp),athet(2,-ntyp:ntyp),-1:1,-1:1),
+      &  bthet(2,-ntyp:ntyp,-1:1,-1:1),polthet(0:3,-ntyp:ntyp),
+      &  gthet(3,-ntyp:ntyp),theta0(-ntyp:ntyp),sig0(-ntyp:ntyp),
+      &  sigc0(-ntyp:ntyp)
  C Parameters of the side-chain probability distribution
        common /sclocal/ dsc(ntyp1),dsc_inv(ntyp1),bsc(maxlob,ntyp),
 -     &  censc(3,maxlob,ntyp),gaussc(3,3,maxlob,ntyp),dsc0(ntyp1),
 +     &  censc(3,maxlob,-ntyp:ntyp),gaussc(3,3,maxlob,-ntyp,ntyp),
 +     &  dsc0(ntyp1),
       &    nlob(ntyp1)
  C Parameters of ab initio-derived potential of virtual-bond-angle bending
        integer nthetyp,ntheterm,ntheterm2,ntheterm3,nsingle,ndouble,
  c      write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
        evdw=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)
  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
  c     print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
        evdw=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)
@@@ -1276,7 -1276,7 +1276,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
@@@ -1383,8 -1383,8 +1383,8 @@@ c     els
  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)
@@@ -1519,8 -1519,8 +1519,8 @@@ c     print *,'Entering EGB nnt=',nnt,
  c     if (icall.eq.0) lprn=.false.
        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)
          do iint=1,nint_gr(i)
            do j=istart(i,iint),iend(i,iint)
              ind=ind+1
 -            itypj=itype(j)
 +            itypj=iabs(itype(j))
  c            dscj_inv=dsc_inv(itypj)
              dscj_inv=vbld_inv(j+nres)
  c            write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
@@@ -1678,8 -1678,8 +1678,8 @@@ c     print *,'Entering EGB nnt=',nnt,
  c     if (icall.eq.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)
          do iint=1,nint_gr(i)
            do j=istart(i,iint),iend(i,iint)
              ind=ind+1
 -            itypj=itype(j)
 +            itypj=iabs(itype(j))
  c            dscj_inv=dsc_inv(itypj)
              dscj_inv=vbld_inv(j+nres)
              sig0ij=sigma(itypi,itypj)
  cd    print *,'Entering Esoft_sphere nnt=',nnt,' nct=',nct
        evdw=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)
  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
@@@ -4032,7 -4032,7 +4032,7 @@@ cd    write (iout,*) 'iatscp_s=',iatscp
          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
@@@ -4126,7 -4126,7 +4126,7 @@@ cd    write (iout,*) 'iatscp_s=',iatscp
          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
@@@ -4242,8 -4242,7 +4242,8 @@@ C iii and jjj point to the residues fo
  cd        write (iout,*) "i",i," ii",ii," iii",iii," jj",jj," jjj",jjj
  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
        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)
        dzi=dc_norm(3,nres+i)
  c      dsci_inv=dsc_inv(itypi)
        dsci_inv=vbld_inv(nres+i)
 -      itypj=itype(j)
 +      itypj=iabs(itype(j))
  c      dscj_inv=dsc_inv(itypj)
        dscj_inv=vbld_inv(nres+j)
        xj=c(1,nres+j)-xi
  c 09/18/07 AL: multimodal bond potential based on AM1 CA-SC PMF's included
  c
        do i=ibond_start,ibond_end
 -        iti=itype(i)
 +        iti=iabs(itype(i))
          if (iti.ne.10) then
            nbi=nbondterm(iti)
            if (nbi.eq.1) then
@@@ -4485,7 -4484,19 +4485,23 @@@ c     write (*,'(a,i2)') 'EBEND ICG=',i
        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)
++<<<<<<< HEAD
 +        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
++>>>>>>> bartek2
          if (i.gt.3) then
  #ifdef OSF
          phii=phi(i)
@@@ -4519,15 -4530,27 +4535,27 @@@ C dependent on the adjacent virtual-bon
  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)
          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)
  #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)
  #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)
@@@ -4883,7 -4906,7 +4911,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
@@@ -5040,11 -5063,11 +5068,11 @@@ C Compute the contribution to SC energ
  
            do j=1,nlobit
  #ifdef OSF
 -            adexp=bsc(j,it)-0.5D0*contr(j,iii)+emin
 +            adexp=bsc(j,iabs(it))-0.5D0*contr(j,iii)+emin
              if(adexp.ne.adexp) adexp=1.0
              expfac=dexp(adexp)
  #else
 -            expfac=dexp(bsc(j,it)-0.5D0*contr(j,iii)+emin)
 +            expfac=dexp(bsc(j,iabs(it))-0.5D0*contr(j,iii)+emin)
  #endif
  cd          print *,'j=',j,' expfac=',expfac
              escloc_i=escloc_i+expfac
@@@ -5126,7 -5149,7 +5154,7 @@@ C Compute the contribution to SC energ
  
        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
@@@ -5609,11 -5632,6 +5637,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...
@@@ -5713,9 -5731,9 +5741,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
@@@ -5730,7 -5748,7 +5758,7 @@@ C          [v2 cos(phi/2)+v3 sin(phi/2)
  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 (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
@@@ -5805,7 -5822,7 +5833,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
@@@ -15,8 -15,8 +15,8 @@@ cd    print *,' CG Processor',me,' maxg
        maxsi=100
  cd    write (iout,*) 'Gen_Rand_conf: nstart=',nstart
        if (nstart.lt.5) then
 -        it1=itype(2)
 -        phi(4)=gen_phi(4,itype(2),itype(3))
 +        it1=iabs(itype(2))
 +        phi(4)=gen_phi(4,iabs(itype(2)),abs(itype(3)))
  c       write(iout,*)'phi(4)=',rad2deg*phi(4)
          if (nstart.lt.3) theta(3)=gen_theta(itype(2),pi,phi(4))
  c       write(iout,*)'theta(3)=',rad2deg*theta(3) 
@@@ -54,9 -54,9 +54,9 @@@
            endif
            return1
          endif
 -      it1=itype(i-1)
 -      it2=itype(i-2)
 -      it=itype(i)
 +      it1=abs(itype(i-1))
 +      it2=abs(itype(i-2))
 +      it=abs(itype(i))
  c       print *,'Gen_Rand_Conf: i=',i,' it=',it,' it1=',it1,' it2=',it2,
  c     &    ' nit=',nit,' niter=',niter,' maxgen=',maxgen
        phi(i+1)=gen_phi(i+1,it1,it)
@@@ -132,12 -132,12 +132,12 @@@ c--------------------------------------
        include 'COMMON.FFIELD'
        data redfac /0.5D0/
        overlap=.false.
 -      iti=itype(i)
 +      iti=abs(itype(i))
        if (iti.gt.ntyp) return
  C Check for SC-SC overlaps.
  cd    print *,'nnt=',nnt,' nct=',nct
        do j=nnt,i-1
 -        itj=itype(j)
 +        itj=abs(itype(j))
          if (j.lt.i-1 .or. ipot.ne.4) then
            rcomp=sigmaii(iti,itj)
          else 
@@@ -159,7 -159,7 +159,7 @@@ C SCs
        c(j,maxres2+1)=0.5D0*(c(j,i)+c(j,i+1))
        enddo
        do j=nnt,i-2
 -      itj=itype(j)
 +      itj=abs(itype(j))
  cd      print *,'overlap, p-Sc: i=',i,' j=',j,
  cd   &         ' dist=',dist(nres+j,maxres2+1)
        if (dist(nres+j,maxres2+1).lt.4.0D0*redfac) then
@@@ -238,7 -238,8 +238,8 @@@ c     print *,'gen_theta: it=',i
        endif  
        thet_pred_mean=a0thet(it)
        do k=1,2
-         thet_pred_mean=thet_pred_mean+athet(k,it)*y(k)+bthet(k,it)*z(k)
+         thet_pred_mean=thet_pred_mean+athet(k,it,1,1)*y(k)
+      &                +bthet(k,it,1,1)*z(k)
        enddo
        sig=polthet(3,it)
        do j=2,0,-1
@@@ -779,7 -780,7 +780,7 @@@ c     overlapping residues left, or fal
  
          do ires=1,ioverlap_last 
            i=ioverlap(ires)
 -          iti=itype(i)
 +          iti=abs(itype(i))
            if (iti.ne.10) then
              nsi=0
              fail=.true.
@@@ -839,8 -840,8 +840,8 @@@ C Check for SC-SC overlaps and mark res
  c      print *,'>>overlap_sc nnt=',nnt,' nct=',nct
        ind=0
        do i=iatsc_s,iatsc_e
 -        itypi=itype(i)
 -        itypi1=itype(i+1)
 +        itypi=abs(itype(i))
 +        itypi1=abs(itype(i+1))
          xi=c(1,nres+i)
          yi=c(2,nres+i)
          zi=c(3,nres+i)
@@@ -164,8 -164,12 +164,12 @@@ c      call memmon_print_usage(
        rr0(i)=0.0D0
        a0thet(i)=0.0D0
        do j=1,2
-         athet(j,i)=0.0D0
-         bthet(j,i)=0.0D0
+          do ichir1=-1,1
+           do ichir2=-1,1
+          athet(j,i,ichir1,ichir2)=0.0D0
+          bthet(j,i,ichir1,ichir2)=0.0D0
+           enddo
+          enddo
          enddo
        do j=0,3
          polthet(j,i)=0.0D0
        do iblock=1,2
        do j=-maxtor,maxtor
          do k=1,maxterm
 -          v1(k,j,i)=0.0D0
 -          v2(k,j,i)=0.0D0
 +          v1(k,j,i,iblock)=0.0D0
 +          v2(k,j,i,iblock)=0.0D0
            enddo
          enddo
        enddo
@@@ -275,13 -279,9 +279,13 @@@ c--------------------------------------
        include 'COMMON.NAMES'
        include 'COMMON.FFIELD'
        data restyp /
 +     &'DD' ,'DPR','DLY','DAR','DHI','DAS','DGL','DSG','DGN','DSN','DTH',
 +     &'DYY','DAL','DTY','DTR','DVA','DLE','DIL','DPN','MED','DCY','ZER',
       &'CYS','MET','PHE','ILE','LEU','VAL','TRP','TYR','ALA','GLY','THR',
       &'SER','GLN','ASN','GLU','ASP','HIS','ARG','LYS','PRO','D'/
        data onelet /
 +     &'z','p','k','r','h','d','e','n','q','s','t','g',
 +     &'a','y','w','v','l','i','f','m','c','x',
       &'C','M','F','I','L','V','W','Y','A','G','T',
       &'S','Q','N','E','D','H','R','K','P','X'/
        data potname /'LJ','LJK','BP','GB','GBV'/
@@@ -103,13 -103,47 +103,47 @@@ C Read the parameters of the probabilit
  C of the virtual-bond valence angles theta
  C
        do i=1,ntyp
-         read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i),j=1,2),
-      &    (bthet(j,i),j=1,2)
+         read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),
+      &    (bthet(j,i,1,1),j=1,2)
          read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
        read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
        read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
        sigc0(i)=sigc0(i)**2
        enddo
+       do i=1,ntyp
+       athet(1,i,1,-1)=athet(1,i,1,1)
+       athet(2,i,1,-1)=athet(2,i,1,1)
+       bthet(1,i,1,-1)=-bthet(1,i,1,1)
+       bthet(2,i,1,-1)=-bthet(2,i,1,1)
+       athet(1,i,-1,1)=-athet(1,i,1,1)
+       athet(2,i,-1,1)=-athet(2,i,1,1)
+       bthet(1,i,-1,1)=bthet(1,i,1,1)
+       bthet(2,i,-1,1)=bthet(2,i,1,1)
+       enddo
+       do i=-ntyp,-1
+       a0thet(i)=a0thet(-i)
+       athet(1,i,-1,-1)=athet(1,-i,1,1)
+       athet(2,i,-1,-1)=-athet(2,-i,1,1)
+       bthet(1,i,-1,-1)=bthet(1,-i,1,1)
+       bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
+       athet(1,i,-1,1)=athet(1,-i,1,1)
+       athet(2,i,-1,1)=-athet(2,-i,1,1)
+       bthet(1,i,-1,1)=-bthet(1,-i,1,1)
+       bthet(2,i,-1,1)=bthet(2,-i,1,1)
+       athet(1,i,1,-1)=-athet(1,-i,1,1)
+       athet(2,i,1,-1)=athet(2,-i,1,1)
+       bthet(1,i,1,-1)=bthet(1,-i,1,1)
+       bthet(2,i,1,-1)=-bthet(2,-i,1,1)
+       theta0(i)=theta0(-i)
+       sig0(i)=sig0(-i)
+       sigc0(i)=sigc0(-i)
+        do j=0,3
+         polthet(j,i)=polthet(j,-i)
+        enddo
+        do j=1,3
+          gthet(j,i)=gthet(j,-i)
+        enddo
+       enddo
        close (ithep)
        if (lprint) then
        if (.not.LaTeX) then
       & '        B1    ','         B2   '        
          do i=1,ntyp
            write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
-      &        a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
+      &        a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
          enddo
          write (iout,'(/a/9x,5a/79(1h-))') 
       & 'Parameters of the expression for sigma(theta_c):',
       & '   b1*10^1    ','    b2*10^1   '        
          do i=1,ntyp
            write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
-      &        a0thet(i),(100*athet(j,i),j=1,2),(10*bthet(j,i),j=1,2)
+      &        a0thet(i),(100*athet(j,i,1,1),j=1,2),
+      $        (10*bthet(j,i,1,1),j=1,2)
          enddo
        write (iout,'(/a/9x,5a/79(1h-))') 
       & 'Parameters of the expression for sigma(theta_c):',
        bsc(1,i)=0.0D0
          read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
       &    ((blower(k,l,1),l=1,k),k=1,3)
+         censc(1,1,-i)=censc(1,1,i)
+         censc(2,1,-i)=censc(2,1,i)
+         censc(3,1,-i)=-censc(3,1,i)
        do j=2,nlob(i)
          read (irotam,*,end=112,err=112) bsc(j,i)
          read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
       &                                 ((blower(k,l,j),l=1,k),k=1,3)
+         censc(1,j,-i)=censc(1,j,i)
+         censc(2,j,-i)=censc(2,j,i)
+         censc(3,j,-i)=-censc(3,j,i)
+ C BSC is amplitude of Gaussian
          enddo
        do j=1,nlob(i)
          do k=1,3
                enddo
              gaussc(k,l,j,i)=akl
              gaussc(l,k,j,i)=akl
+               if (((k.eq.3).and.(l.ne.3))
+      &        .or.((l.eq.3).and.(k.ne.3))) then
+                 gaussc(k,l,j,-i)=-akl
+                 gaussc(l,k,j,-i)=-akl
+               else
+                 gaussc(k,l,j,-i)=akl
+                 gaussc(l,k,j,-i)=akl
+               endif
              enddo
            enddo 
        enddo
        read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
        do iblock=1,2
        do i=-ntyp,-1
++<<<<<<< HEAD
 +       itortyp(i)=-itortyp(-1)
 +      enddo
 +c      write (iout,*) 'ntortyp',ntortyp
 +      do i=0,ntortyp,ntortyp-1
 +      do j=-ntortyp,ntortyp
 +        read (itorp,*,end=113,err=113) nterm(i,j,iblock),
 +     &          nlor(i,j,iblock)
 +           nterm(-i,-j,iblock)=nterm(i,j,iblock)
 +           nlor(-i,-j,iblock)=nlor(i,j,iblock)
++=======
+        itortyp(i)=-itortyp(-i)
+       enddo
+ c      write (iout,*) 'ntortyp',ntortyp
+       do i=0,ntortyp-1
+       do j=-ntortyp,ntortyp
+         read (itorp,*,end=113,err=113) nterm(i,j,iblock),
+                  nlor(i,j,iblock)
+           nterm(-i,-j,iblock)=nterm(i,j,iblock)
+      &          nlor(i,j,iblock)
++>>>>>>> bartek2
            v0ij=0.0d0
            si=-1.0d0
          do k=1,nterm(i,j,iblock)
            read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
       &      v2(k,i,j,iblock)
              v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
++<<<<<<< HEAD
 +            v2(k,-i,-j,iblock)=-v2(k,i,j,iblock) 
++=======
+             v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
++>>>>>>> bartek2
              v0ij=v0ij+si*v1(k,i,j,iblock)
              si=-si
+ c         write(iout,*) i,j,k,iblock,nterm(i,j,iblock)
+ c         write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),
+ c      &v2(k,-i,-j,iblock),v2(k,i,j,iblock)
            enddo
          do k=1,nlor(i,j,iblock)
              read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
       &        vlor2(k,i,j),vlor3(k,i,j) 
              v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
            enddo
-           v0(i,j)=v0ij
+           v0(i,j,iblock)=v0ij
+           v0(-i,-j,iblock)=v0ij
          enddo
        enddo
+       enddo
        close (itorp)
        if (lprint) then
        write (iout,'(/a/)') 'Torsional constants:'
          do j=1,ntortyp
              write (iout,*) 'ityp',i,' jtyp',j
              write (iout,*) 'Fourier constants'
 -            do k=1,nterm(i,j)
 -            write (iout,'(2(1pe15.5))') v1(k,i,j),v2(k,i,j)
 +            do k=1,nterm(i,j,iblock)
 +            write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
 +     &         v2(k,i,j,iblock)
              enddo
              write (iout,*) 'Lorenz constants'
 -            do k=1,nlor(i,j)
 +            do k=1,nlor(i,j,iblock)
              write (iout,'(3(1pe15.5))') 
       &         vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
              enddo
              endif
              read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
       &         ntermd_2(i,j,k,iblock)
+             ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
+             ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
              read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
       &         ntermd_1(i,j,k,iblock))
              read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
              read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
       &         ntermd_1(i,j,k,iblock))
  C Matrix of D parameters for one dimesional foureir series
-             do l=1,  ntermd_1(i,j,k,iblock)
+             do l=1,ntermd_1(i,j,k,iblock)
               v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
               v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
               v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
               v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
+ c            write(iout,*) "whcodze" ,
+ c     & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
              enddo
              read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
       &         v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
@@@ -553,9 -612,10 +629,10 @@@ C Matrix of D parameters for two dimesi
          enddo!j
        enddo!i
        enddo!iblock
- cc      if (lprint) then
+       if (lprint) then
        write (iout,*) 
        write (iout,*) 'Constants for double torsionals'
+       do iblock=1,2
        do i=1,ntortyp
          do j=-ntortyp,ntortyp
            do k=-ntortyp,ntortyp
            enddo
          enddo
        enddo
- cc      endif
+       enddo
+       endif
  #endif
  C
  C 5/21/07 (AL) Read coefficients of the backbone-local sidechain-local
          endif
          B1(1,i)  = b(3)
          B1(2,i)  = b(5)
+         B1(1,-i) = b(3)
+         B1(2,-i) = -b(5)
  c        b1(1,i)=0.0d0
  c        b1(2,i)=0.0d0
          B1tilde(1,i) = b(3)
-         B1tilde(2,i) =-b(5) 
+         B1tilde(2,i) =-b(5)
+         B1tilde(1,-i) =-b(3)
+         B1tilde(2,-i) =b(5)
  c        b1tilde(1,i)=0.0d0
  c        b1tilde(2,i)=0.0d0
          B2(1,i)  = b(2)
          B2(2,i)  = b(4)
+         B2(1,-i)  =b(2)
+         B2(2,-i)  =-b(4)
  c        b2(1,i)=0.0d0
  c        b2(2,i)=0.0d0
          CC(1,1,i)= b(7)
          CC(2,2,i)=-b(7)
          CC(2,1,i)= b(9)
          CC(1,2,i)= b(9)
+         CC(1,1,-i)= b(7)
+         CC(2,2,-i)=-b(7)
+         CC(2,1,-i)=-b(9)
+         CC(1,2,-i)=-b(9)
  c        CC(1,1,i)=0.0d0
  c        CC(2,2,i)=0.0d0
  c        CC(2,1,i)=0.0d0
@@@ -656,6 -728,10 +745,10 @@@ c        CC(1,2,i)=0.0d
          Ctilde(1,2,i)=b(9)
          Ctilde(2,1,i)=-b(9)
          Ctilde(2,2,i)=b(7)
+         Ctilde(1,1,-i)=b(7)
+         Ctilde(1,2,-i)=-b(9)
+         Ctilde(2,1,-i)=b(9)
+         Ctilde(2,2,-i)=b(7)
  c        Ctilde(1,1,i)=0.0d0
  c        Ctilde(1,2,i)=0.0d0
  c        Ctilde(2,1,i)=0.0d0
@@@ -664,6 -740,10 +757,10 @@@ c        Ctilde(2,2,i)=0.0d
          DD(2,2,i)=-b(6)
          DD(2,1,i)= b(8)
          DD(1,2,i)= b(8)
+         DD(1,1,-i)= b(6)
+         DD(2,2,-i)=-b(6)
+         DD(2,1,-i)=-b(8)
+         DD(1,2,-i)=-b(8)
  c        DD(1,1,i)=0.0d0
  c        DD(2,2,i)=0.0d0
  c        DD(2,1,i)=0.0d0
@@@ -672,6 -752,10 +769,10 @@@ c        DD(1,2,i)=0.0d
          Dtilde(1,2,i)=b(8)
          Dtilde(2,1,i)=-b(8)
          Dtilde(2,2,i)=b(6)
+         Dtilde(1,1,-i)=b(6)
+         Dtilde(1,2,-i)=-b(8)
+         Dtilde(2,1,-i)=b(8)
+         Dtilde(2,2,-i)=b(6)
  c        Dtilde(1,1,i)=0.0d0
  c        Dtilde(1,2,i)=0.0d0
  c        Dtilde(2,1,i)=0.0d0
@@@ -680,6 -764,10 +781,10 @@@ c        Dtilde(2,2,i)=0.0d
          EE(2,2,i)=-b(10)+b(11)
          EE(2,1,i)= b(12)-b(13)
          EE(1,2,i)= b(12)+b(13)
+         EE(1,1,-i)= b(10)+b(11)
+         EE(2,2,-i)=-b(10)+b(11)
+         EE(2,1,-i)=-b(12)+b(13)
+         EE(1,2,-i)=-b(12)-b(13)
  c        ee(1,1,i)=1.0d0
  c        ee(2,2,i)=1.0d0
  c        ee(2,1,i)=0.0d0
@@@ -493,8 -493,8 +493,8 @@@ C Assign initial virtual bond length
            vbld_inv(i)=vblinv
          enddo
          do i=2,nres-1
 -          vbld(i+nres)=dsc(itype(i))
 -          vbld_inv(i+nres)=dsc_inv(itype(i))
 +          vbld(i+nres)=dsc(iabs(itype(i)))
 +          vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
  c          write (iout,*) "i",i," itype",itype(i),
  c     &      " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
          enddo
@@@ -509,9 -509,9 +509,9 @@@ c      print '(20i4)',(itype(i),i=1,nre
  #endif
            itel(i)=0
  #ifdef PROCOR
-         else if (itype(i+1).ne.20) then
+         else if (iabs(itype(i+1)).ne.20) then
  #else
-         else if (itype(i).ne.20) then
+         else if (iabs(itype(i)).ne.20) then
  #endif
          itel(i)=1
          else
@@@ -730,7 -730,6 +730,7 @@@ C initial geometry
           enddo
           do i=2,nres-1
            omeg(i)=-120d0*deg2rad
 +          if (itype(i).le.0) omeg(i)=-omeg(i)
           enddo
          else
            if(me.eq.king.or..not.out1file)