adding ebend_nucl to UCGM+some further reading
[unres4.git] / source / unres / geometry.f90
index 0f3febc..6895472 100644 (file)
@@ -1,4 +1,4 @@
-      module geometry
+             module geometry
 !-----------------------------------------------------------------------------
       use io_units
       use names
@@ -83,6 +83,7 @@
       nres2=2*nres
 ! Set lprn=.true. for debugging
       lprn = .false.
+      print *,"I ENTER CHAINBUILD"
 !
 ! Define the origin and orientation of the coordinate system and locate the
 ! first three CA's and SC(2).
         be1=rad2deg*beta(nres+i,i,nres2+2,i+1)
         alfai=0.0D0
         if (i.gt.2) alfai=rad2deg*alpha(i-2,i-1,i)
-        write (iout,1212) restyp(itype(i)),i,dist(i-1,i),&
+        write (iout,1212) restyp(itype(i,1),1),i,dist(i-1,i),&
         alfai,be,dist(nres+i,i),rad2deg*alpha(nres+i,i,nres2+2),be1
       enddo   
  1212 format (a3,'(',i3,')',2(f10.5,2f10.2))
       real(kind=8) :: alphi,omegi,theta2
       real(kind=8) :: dsci,dsci_inv,sinalphi,cosalphi,cosomegi,sinomegi
       real(kind=8) :: xp,yp,zp,cost2,sint2,rj
-!      dsci=dsc(itype(i))
-!      dsci_inv=dsc_inv(itype(i))
+!      dsci=dsc(itype(i,1))
+!      dsci_inv=dsc_inv(itype(i,1))
       dsci=vbld(i+nres)
       dsci_inv=vbld_inv(i+nres)
 #ifdef OSF
       xx(1)= xp*cost2+yp*sint2
       xx(2)=-xp*sint2+yp*cost2
       xx(3)= zp
-!d    print '(a3,i3,3f10.5,5x,3f10.5)',restyp(itype(i)),i,
+!d    print '(a3,i3,3f10.5,5x,3f10.5)',restyp(itype(i,1)),i,
 !d   &   xp,yp,zp,(xx(k),k=1,3)
       do j=1,3
         xloc(j,i)=xx(j)
         be=0.0D0
         if (i.gt.2) then
         if (i.le.nres) phi(i+1)=beta(i-2,i-1,i,i+1)
-        if ((itype(i).ne.10).and.(itype(i-1).ne.10)) then
+        if ((itype(i,1).ne.10).and.(itype(i-1,1).ne.10)) then
          tauangle(3,i+1)=beta(i+nres-1,i-1,i,i+nres)
         endif
-        if (itype(i-1).ne.10) then
+        if (itype(i-1,1).ne.10) then
          tauangle(1,i+1)=beta(i-1+nres,i-1,i,i+1)
          omicron(1,i)=alpha(i-2,i-1,i-1+nres)
          omicron(2,i)=alpha(i-1+nres,i-1,i)
         endif
-        if (itype(i).ne.10) then
+        if (itype(i,1).ne.10) then
          tauangle(2,i+1)=beta(i-2,i-1,i,i+nres)
         endif
         endif
         vbld(i)=dist(i-1,i)
         vbld_inv(i)=1.0d0/vbld(i)
         vbld(nres+i)=dist(nres+i,i)
-        if (itype(i).ne.10) then
+        if (itype(i,1).ne.10) then
           vbld_inv(nres+i)=1.0d0/vbld(nres+i)
         else
           vbld_inv(nres+i)=0.0d0
       enddo
       if (lprn) then
       do i=2,nres
-       write (iout,1212) restyp(itype(i)),i,vbld(i),&
+       write (iout,1212) restyp(itype(i,1),1),i,vbld(i),&
        rad2deg*theta(i),rad2deg*phi(i),vbld(nres+i),&
        rad2deg*alph(i),rad2deg*omeg(i)
       enddo
       print *,'dv=',dv
       do 10 it=1,1 
         if (it.eq.10) goto 10 
-        open (20,file=restyp(it)//'_distr.sdc',status='unknown')
+        open (20,file=restyp(it,1)//'_distr.sdc',status='unknown')
         call gen_side(it,90.0D0 * deg2rad,al,om,fail)
         close (20)
         goto 10
-        open (20,file=restyp(it)//'_distr1.sdc',status='unknown')
+        open (20,file=restyp(it,1)//'_distr1.sdc',status='unknown')
         do i=0,90
           do j=0,72
             prob(j,i)=0.0D0
       maxsi=100
 !d    write (iout,*) 'Gen_Rand_conf: nstart=',nstart
       if (nstart.lt.5) then
-        it1=iabs(itype(2))
-        phi(4)=gen_phi(4,iabs(itype(2)),iabs(itype(3)))
+        it1=iabs(itype(2,1))
+        phi(4)=gen_phi(4,iabs(itype(2,1)),iabs(itype(3,1)))
 !       write(iout,*)'phi(4)=',rad2deg*phi(4)
-        if (nstart.lt.3) theta(3)=gen_theta(iabs(itype(2)),pi,phi(4))
+        if (nstart.lt.3) theta(3)=gen_theta(iabs(itype(2,1)),pi,phi(4))
 !       write(iout,*)'theta(3)=',rad2deg*theta(3) 
         if (it1.ne.10) then
           nsi=0
           endif
           return 1
         endif
-       it1=iabs(itype(i-1))
-       it2=iabs(itype(i-2))
-       it=iabs(itype(i))
+       it1=iabs(itype(i-1,1))
+       it2=iabs(itype(i-2,1))
+       it=iabs(itype(i,1))
 !       print *,'Gen_Rand_Conf: i=',i,' it=',it,' it1=',it1,' it2=',it2,
 !    &    ' nit=',nit,' niter=',niter,' maxgen=',maxgen
        phi(i+1)=gen_phi(i+1,it1,it)
       nres2=2*nres
       data redfac /0.5D0/
       overlap=.false.
-      iti=iabs(itype(i))
+      iti=iabs(itype(i,1))
       if (iti.gt.ntyp) return
 ! Check for SC-SC overlaps.
 !d    print *,'nnt=',nnt,' nct=',nct
       do j=nnt,i-1
-        itj=iabs(itype(j))
+        itj=iabs(itype(j,1))
         if (j.lt.i-1 .or. ipot.ne.4) then
           rcomp=sigmaii(iti,itj)
         else 
        c(j,nres2+3)=0.5D0*(c(j,i)+c(j,i+1))
       enddo
       do j=nnt,i-2
-       itj=iabs(itype(j))
+       itj=iabs(itype(j,1))
 !d      print *,'overlap, p-Sc: i=',i,' j=',j,
 !d   &         ' dist=',dist(nres+j,maxres2+1)
        if (dist(nres+j,nres2+3).lt.4.0D0*redfac) then
 
         do ires=1,ioverlap_last 
           i=ioverlap(ires)
-          iti=iabs(itype(i))
+          iti=iabs(itype(i,1))
           if (iti.ne.10) then
             nsi=0
             fail=.true.
 !      print *,'>>overlap_sc nnt=',nnt,' nct=',nct
       ind=0
       do i=iatsc_s,iatsc_e
-        itypi=iabs(itype(i))
-        itypi1=iabs(itype(i+1))
+        itypi=iabs(itype(i,1))
+        itypi1=iabs(itype(i+1,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=iabs(itype(j))
+            itypj=iabs(itype(j,1))
             dscj_inv=dsc_inv(itypj)
             sig0ij=sigma(itypi,itypj)
             chi1=chi(itypi,itypj)
        endif
       endif
       do i=1,nres-1
+       if (molnum(i).ne.1) cycle
 !in wham      do i=1,nres
-        iti=itype(i)
-        if (dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0) then
+        iti=itype(i,1)
+        if ((dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0).and.&
+       (iti.ne.ntyp1  .and. itype(i+1,1).ne.ntyp1)) then
           write (iout,'(a,i4)') 'Bad Cartesians for residue',i
 !test          stop
         endif
       enddo
 !el -----
 !#ifdef WHAM_RUN
-!      if (itype(1).eq.ntyp1) then
+!      if (itype(1,1).eq.ntyp1) then
 !        do j=1,3
 !          c(j,1)=c(j,2)+(c(j,3)-c(j,4))
 !        enddo
 !      endif
-!      if (itype(nres).eq.ntyp1) then
+!      if (itype(nres,1).eq.ntyp1) then
 !        do j=1,3
 !          c(j,nres)=c(j,nres-1)+(c(j,nres-2)-c(j,nres-3))
 !        enddo
 !      endif
 !#endif
 !      if (unres_pdb) then
-!        if (itype(1).eq.21) then
+!        if (itype(1,1).eq.21) then
 !          theta(3)=90.0d0*deg2rad
 !          phi(4)=180.0d0*deg2rad
 !          vbld(2)=3.8d0
 !          vbld_inv(2)=1.0d0/vbld(2)
 !        endif
-!        if (itype(nres).eq.21) then
+!        if (itype(nres,1).eq.21) then
 !          theta(nres)=90.0d0*deg2rad
 !          phi(nres)=180.0d0*deg2rad
 !          vbld(nres)=3.8d0
            +(c(j,i+1)-c(j,i))*vbld_inv(i+1))
 ! in wham            c(j,maxres2)=0.5D0*(c(j,i-1)+c(j,i+1)
           enddo
-          iti=itype(i)
+          iti=itype(i,1)
           di=dist(i,nres+i)
 !#ifndef WHAM_RUN
 ! 10/03/12 Adam: Correction for zero SC-SC bond length
-          if (itype(i).ne.10 .and. itype(i).ne.ntyp1 .and. di.eq.0.0d0) &
-           di=dsc(itype(i))
+          
+          if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1 .and. di.eq.0.0d0) &
+           di=dsc(itype(i,molnum(i)))
           vbld(i+nres)=di
-          if (itype(i).ne.10) then
+          if (itype(i,1).ne.10) then
             vbld_inv(i+nres)=1.0d0/di
           else
             vbld_inv(i+nres)=0.0d0
             alph(i)=alpha(nres+i,i,nres2+2)
             omeg(i)=beta(nres+i,i,nres2+2,i+1)
           endif
+          if (iti.ne.0) then
           if(me.eq.king.or..not.out1file)then
            if (lprn) &
-           write (iout,'(a3,i4,7f10.3)') restyp(iti),i,vbld(i),&
+           write (iout,'(a3,i4,7f10.3)') restyp(iti,1),i,vbld(i),&
            rad2deg*theta(i),rad2deg*phi(i),dsc(iti),vbld(nres+i),&
            rad2deg*alph(i),rad2deg*omeg(i)
           endif
+          else
+          if(me.eq.king.or..not.out1file)then
+           if (lprn) &
+           write (iout,'(a3,i4,7f10.3)') restyp(iti,1),i,vbld(i),&
+           rad2deg*theta(i),rad2deg*phi(i),dsc(iti+1),vbld(nres+i),&
+           rad2deg*alph(i),rad2deg*omeg(i)
+          endif
+          endif
         enddo
       else if (lprn) then
         do i=2,nres
-          iti=itype(i)
+          iti=itype(i,1)
           if(me.eq.king.or..not.out1file) &
-           write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),&
+           write (iout,'(a3,i4,7f10.3)') restyp(iti,1),i,dist(i,i-1),&
            rad2deg*theta(i),rad2deg*phi(i)
         enddo
       endif
         enddo
       enddo
       do i=2,nres-1
-        if (itype(i).ne.10) then
+        if (itype(i,1).ne.10) then
           do j=1,3
             dc_norm(j,i+nres)=vbld_inv(i+nres)*(c(j,i+nres)-c(j,i))
           enddo
         cosfac=dsqrt(cosfac2)
         sinfac2=0.5d0/(1.0d0-costtab(i+1))
         sinfac=dsqrt(sinfac2)
-        it=itype(i)
+        it=itype(i,1)
 
         if ((it.ne.10).and.(it.ne.ntyp1)) then
 !el        if (it.ne.10) then
       enddo
       if (lprn) then
         do i=2,nres
-          iti=itype(i)
+          iti=itype(i,1)
           if(me.eq.king.or..not.out1file) &
-           write (iout,'(a3,i4,3f10.5)') restyp(iti),i,xxref(i),&
+           write (iout,'(a3,i4,3f10.5)') restyp(iti,1),i,xxref(i),&
             yyref(i),zzref(i)
         enddo
       endif
       integer :: i,j,ires,nscat
       real(kind=8),dimension(3,20) :: sccor
       real(kind=8) :: sccmj
+!        print *,"I am in sccenter",ires,nscat
       do j=1,3
         sccmj=0.0D0
         do i=1,nscat
-          sccmj=sccmj+sccor(j,i) 
+          sccmj=sccmj+sccor(j,i)
+          print *,"insccent", ires,sccor(j,i) 
         enddo
         dc(j,ires)=sccmj/nscat
       enddo
       do i=1,nres-1
        vbld(i+1)=vbl
        vbld_inv(i+1)=1.0d0/vbld(i+1)
-       vbld(i+1+nres)=dsc(itype(i+1))
-       vbld_inv(i+1+nres)=dsc_inv(itype(i+1))
+       vbld(i+1+nres)=dsc(itype(i+1,1))
+       vbld_inv(i+1+nres)=dsc_inv(itype(i+1,1))
 !       print *,vbld(i+1),vbld(i+1+nres)
       enddo
       return
        do j=1,3
          gcart(j,1)=gcart(j,1)+gloc(1,icg)*dphi(j,1,4) &
            +gloc(nres-2,icg)*dtheta(j,1,3)      
-         if(itype(2).ne.10) then
+         if(itype(2,1).ne.10) then
           gcart(j,1)=gcart(j,1)+gloc(ialph(2,1),icg)*dalpha(j,1,2)+ &
           gloc(ialph(2,1)+nside,icg)*domega(j,1,2)             
         endif
        do j=1,3
          gcart(j,2)=gcart(j,2)+gloc(1,icg)*dphi(j,2,4)+ &
         gloc(nres-2,icg)*dtheta(j,2,3)+gloc(nres-1,icg)*dtheta(j,1,4)
-        if(itype(2).ne.10) then
+        if(itype(2,1).ne.10) then
           gcart(j,2)=gcart(j,2)+gloc(ialph(2,1),icg)*dalpha(j,2,2)+ &
           gloc(ialph(2,1)+nside,icg)*domega(j,2,2)
         endif
-               if(itype(3).ne.10) then
+               if(itype(3,1).ne.10) then
          gcart(j,2)=gcart(j,2)+gloc(ialph(3,1),icg)*dalpha(j,1,3)+ &
           gloc(ialph(3,1)+nside,icg)*domega(j,1,3)
         endif
            gcart(j,3)=gcart(j,3)+gloc(1,icg)*dphi(j,3,4)+gloc(2,icg)* &
            dphi(j,2,5)+gloc(nres-1,icg)*dtheta(j,2,4)+gloc(nres,icg)* &
            dtheta(j,1,5)
-         if(itype(3).ne.10) then
+         if(itype(3,1).ne.10) then
           gcart(j,3)=gcart(j,3)+gloc(ialph(3,1),icg)* &
           dalpha(j,2,3)+gloc(ialph(3,1)+nside,icg)*domega(j,2,3)
          endif
-        if(itype(4).ne.10) then
+        if(itype(4,1).ne.10) then
           gcart(j,3)=gcart(j,3)+gloc(ialph(4,1),icg)* &
           dalpha(j,1,4)+gloc(ialph(4,1)+nside,icg)*domega(j,1,4)
          endif
           +gloc(i-1,icg)*dphi(j,2,i+2)+ &
           gloc(i,icg)*dphi(j,1,i+3)+gloc(nres+i-4,icg)*dtheta(j,2,i+1)+ &
           gloc(nres+i-3,icg)*dtheta(j,1,i+2)
-          if(itype(i).ne.10) then
+          if(itype(i,1).ne.10) then
            gcart(j,i)=gcart(j,i)+gloc(ialph(i,1),icg)*dalpha(j,2,i)+ &
            gloc(ialph(i,1)+nside,icg)*domega(j,2,i)
           endif
-          if(itype(i+1).ne.10) then
+          if(itype(i+1,1).ne.10) then
            gcart(j,i)=gcart(j,i)+gloc(ialph(i+1,1),icg)*dalpha(j,1,i+1) &
            +gloc(ialph(i+1,1)+nside,icg)*domega(j,1,i+1)
           endif
           dphi(j,3,nres-1)+gloc(nres-3,icg)*dphi(j,2,nres) &
            +gloc(2*nres-6,icg)* &
            dtheta(j,2,nres-1)+gloc(2*nres-5,icg)*dtheta(j,1,nres)
-          if(itype(nres-2).ne.10) then
+          if(itype(nres-2,1).ne.10) then
               gcart(j,nres-2)=gcart(j,nres-2)+gloc(ialph(nres-2,1),icg)* &
              dalpha(j,2,nres-2)+gloc(ialph(nres-2,1)+nside,icg)* &
               domega(j,2,nres-2)
           endif
-          if(itype(nres-1).ne.10) then
+          if(itype(nres-1,1).ne.10) then
              gcart(j,nres-2)=gcart(j,nres-2)+gloc(ialph(nres-1,1),icg)* &
             dalpha(j,1,nres-1)+gloc(ialph(nres-1,1)+nside,icg)* &
              domega(j,1,nres-1)
        do j=1,3
         gcart(j,nres-1)=gcart(j,nres-1)+gloc(nres-3,icg)*dphi(j,3,nres)+ &
        gloc(2*nres-5,icg)*dtheta(j,2,nres)
-        if(itype(nres-1).ne.10) then
+        if(itype(nres-1,1).ne.10) then
           gcart(j,nres-1)=gcart(j,nres-1)+gloc(ialph(nres-1,1),icg)* &
          dalpha(j,2,nres-1)+gloc(ialph(nres-1,1)+nside,icg)* &
           domega(j,2,nres-1)
         enddo
 !   The side-chain vector derivatives
         do i=2,nres-1
-         if(itype(i).ne.10 .and. itype(i).ne.ntyp1) then       
+         if(itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then   
             do j=1,3   
               gxcart(j,i)=gxcart(j,i)+gloc(ialph(i,1),icg)*dalpha(j,3,i) &
               +gloc(ialph(i,1)+nside,icg)*domega(j,3,i)
 !          write (iout,*) "poczotkoawy",i,gloc_sc(1,i,icg)
 !       enddo
        if (nres.lt.2) return
-       if ((nres.lt.3).and.(itype(1).eq.10)) return
-       if ((itype(1).ne.10).and.(itype(1).ne.ntyp1)) then
+       if ((nres.lt.3).and.(itype(1,1).eq.10)) return
+       if ((itype(1,1).ne.10).and.(itype(1,1).ne.ntyp1)) then
         do j=1,3
 !c Derviative was calculated for oposite vector of side chain therefore
 ! there is "-" sign before gloc_sc
            dtauangle(j,1,1,3)
          gcart(j,1)=gcart(j,1)+gloc_sc(1,0,icg)* &
            dtauangle(j,1,2,3)
-          if ((itype(2).ne.10).and.(itype(2).ne.ntyp1)) then
+          if ((itype(2,1).ne.10).and.(itype(2,1).ne.ntyp1)) then
          gxcart(j,1)= gxcart(j,1) &
                      -gloc_sc(3,0,icg)*dtauangle(j,3,1,3)
          gcart(j,1)=gcart(j,1)+gloc_sc(3,0,icg)* &
           endif
        enddo
        endif
-         if ((nres.ge.3).and.(itype(3).ne.10).and.(itype(3).ne.ntyp1)) &
+         if ((nres.ge.3).and.(itype(3,1).ne.10).and.(itype(3,1).ne.ntyp1)) &
       then
          do j=1,3
          gcart(j,1)=gcart(j,1)+gloc_sc(2,1,icg)*dtauangle(j,2,1,4)
 
 !     Calculating the remainder of dE/ddc2
        do j=1,3
-         if((itype(2).ne.10).and.(itype(2).ne.ntyp1)) then
-           if (itype(1).ne.10) gxcart(j,2)=gxcart(j,2)+ &
+         if((itype(2,1).ne.10).and.(itype(2,1).ne.ntyp1)) then
+           if (itype(1,1).ne.10) gxcart(j,2)=gxcart(j,2)+ &
                                gloc_sc(3,0,icg)*dtauangle(j,3,3,3)
-        if ((itype(3).ne.10).and.(nres.ge.3).and.(itype(3).ne.ntyp1)) &
+        if ((itype(3,1).ne.10).and.(nres.ge.3).and.(itype(3,1).ne.ntyp1)) &
          then
            gxcart(j,2)=gxcart(j,2)-gloc_sc(3,1,icg)*dtauangle(j,3,1,4)
 !c                  the   - above is due to different vector direction
 !           write(iout,*) gloc_sc(1,1,icg),dtauangle(j,1,1,4),"gx"
           endif
          endif
-         if ((itype(1).ne.10).and.(itype(1).ne.ntyp1)) then
+         if ((itype(1,1).ne.10).and.(itype(1,1).ne.ntyp1)) then
           gcart(j,2)=gcart(j,2)+gloc_sc(1,0,icg)*dtauangle(j,1,3,3)
 !           write(iout,*)  gloc_sc(1,0,icg),dtauangle(j,1,3,3)
         endif
-         if ((itype(3).ne.10).and.(nres.ge.3)) then
+         if ((itype(3,1).ne.10).and.(nres.ge.3)) then
           gcart(j,2)=gcart(j,2)+gloc_sc(2,1,icg)*dtauangle(j,2,2,4)
 !           write(iout,*) gloc_sc(2,1,icg),dtauangle(j,2,2,4)
          endif
-         if ((itype(4).ne.10).and.(nres.ge.4)) then
+         if ((itype(4,1).ne.10).and.(nres.ge.4)) then
           gcart(j,2)=gcart(j,2)+gloc_sc(2,2,icg)*dtauangle(j,2,1,5)
 !           write(iout,*) gloc_sc(2,2,icg),dtauangle(j,2,1,5)
          endif
 
-!      write(iout,*) gcart(j,2),itype(2),itype(1),itype(3), "gcart2"
+!      write(iout,*) gcart(j,2),itype(2,1),itype(1,1),itype(3,1), "gcart2"
        enddo
 !    If there are more than five residues
       if(nres.ge.5) then                        
         do i=3,nres-2
          do j=1,3
 !          write(iout,*) "before", gcart(j,i)
-          if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
+          if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
           gxcart(j,i)=gxcart(j,i)+gloc_sc(2,i-2,icg) &
           *dtauangle(j,2,3,i+1) &
           -gloc_sc(1,i-1,icg)*dtauangle(j,1,1,i+2)
           *dtauangle(j,1,2,i+2)
 !                   write(iout,*) "new",j,i,
 !     &  gcart(j,i),gloc_sc(1,i-1,icg),dtauangle(j,1,2,i+2)
-          if (itype(i-1).ne.10) then
+          if (itype(i-1,1).ne.10) then
            gxcart(j,i)=gxcart(j,i)+gloc_sc(3,i-2,icg) &
       *dtauangle(j,3,3,i+1)
           endif
-          if (itype(i+1).ne.10) then
+          if (itype(i+1,1).ne.10) then
            gxcart(j,i)=gxcart(j,i)-gloc_sc(3,i-1,icg) &
       *dtauangle(j,3,1,i+2)
            gcart(j,i)=gcart(j,i)+gloc_sc(3,i-1,icg) &
       *dtauangle(j,3,2,i+2)
           endif
           endif
-          if (itype(i-1).ne.10) then
+          if (itype(i-1,1).ne.10) then
            gcart(j,i)=gcart(j,i)+gloc_sc(1,i-2,icg)* &
            dtauangle(j,1,3,i+1)
           endif
-          if (itype(i+1).ne.10) then
+          if (itype(i+1,1).ne.10) then
            gcart(j,i)=gcart(j,i)+gloc_sc(2,i-1,icg)* &
            dtauangle(j,2,2,i+2)
 !          write(iout,*) "numer",i,gloc_sc(2,i-1,icg),
 !     &    dtauangle(j,2,2,i+2)
           endif
-          if (itype(i+2).ne.10) then
+          if (itype(i+2,1).ne.10) then
            gcart(j,i)=gcart(j,i)+gloc_sc(2,i,icg)* &
            dtauangle(j,2,1,i+3)
           endif
 !  Setting dE/ddnres-1       
       if(nres.ge.4) then
          do j=1,3
-         if ((itype(nres-1).ne.10).and.(itype(nres-1).ne.ntyp1)) then
+         if ((itype(nres-1,1).ne.10).and.(itype(nres-1,1).ne.ntyp1)) then
          gxcart(j,nres-1)=gxcart(j,nres-1)+gloc_sc(2,nres-3,icg) &
           *dtauangle(j,2,3,nres)
 !          write (iout,*) "gxcart(nres-1)", gloc_sc(2,nres-3,icg),
 !     &     dtauangle(j,2,3,nres), gxcart(j,nres-1)
-         if (itype(nres-2).ne.10) then
+         if (itype(nres-2,1).ne.10) then
         gxcart(j,nres-1)=gxcart(j,nres-1)+gloc_sc(3,nres-3,icg) &
           *dtauangle(j,3,3,nres)
           endif
-         if ((itype(nres).ne.10).and.(itype(nres).ne.ntyp1)) then
+         if ((itype(nres,1).ne.10).and.(itype(nres,1).ne.ntyp1)) then
         gxcart(j,nres-1)=gxcart(j,nres-1)-gloc_sc(3,nres-2,icg) &
           *dtauangle(j,3,1,nres+1)
         gcart(j,nres-1)=gcart(j,nres-1)+gloc_sc(3,nres-2,icg) &
           *dtauangle(j,3,2,nres+1)
           endif
          endif
-         if ((itype(nres-2).ne.10).and.(itype(nres-2).ne.ntyp1)) then
+         if ((itype(nres-2,1).ne.10).and.(itype(nres-2,1).ne.ntyp1)) then
             gcart(j,nres-1)=gcart(j,nres-1)+gloc_sc(1,nres-3,icg)* &
          dtauangle(j,1,3,nres)
          endif
-          if ((itype(nres).ne.10).and.(itype(nres).ne.ntyp1)) then
+          if ((itype(nres,1).ne.10).and.(itype(nres,1).ne.ntyp1)) then
             gcart(j,nres-1)=gcart(j,nres-1)+gloc_sc(2,nres-2,icg)* &
            dtauangle(j,2,2,nres+1)
 !           write (iout,*) "gcart(nres-1)", gloc_sc(2,nres-2,icg),
-!     &     dtauangle(j,2,2,nres+1), itype(nres-1),itype(nres)
+!     &     dtauangle(j,2,2,nres+1), itype(nres-1,1),itype(nres,1)
            endif
          enddo
       endif
 !  Settind dE/ddnres       
-       if ((nres.ge.3).and.(itype(nres).ne.10).and. &
-          (itype(nres).ne.ntyp1))then
+       if ((nres.ge.3).and.(itype(nres,1).ne.10).and. &
+          (itype(nres,1).ne.ntyp1))then
        do j=1,3
         gxcart(j,nres)=gxcart(j,nres)+gloc_sc(3,nres-2,icg) &
        *dtauangle(j,3,3,nres+1)+gloc_sc(2,nres-2,icg) &
 
       write (iout,100)
       do i=1,nres
-        write (iout,110) restyp(itype(i)),i,c(1,i),c(2,i),&
+        write (iout,110) restyp(itype(i,1),1),i,c(1,i),c(2,i),&
           c(3,i),c(1,nres+i),c(2,nres+i),c(3,nres+i)
       enddo
   100 format (//'              alpha-carbon coordinates       ',&