working gradient for cations

[unres4.git] / source / unres / compare.F90

diff --git a/source/unres/compare.F90 b/source/unres/compare.F90
index 34f8a2a..45ed3e5 100644 (file)
--- a/source/unres/compare.F90
+++ b/source/unres/compare.F90
@@ -41,14 +41,25 @@
       logical :: lprint
 !el local variables
       real(kind=8) :: co,rcomp
-      integer :: kkk,i,j,i1,i2,it1,it2,iti,itj
+      integer :: kkk,i,j,i1,i2,it1,it2,iti,itj,inum,jnum
 
       ncont=0
       kkk=3
       do i=nnt+kkk,nct
-        iti=iabs(itype(i,1))
+        iti=iabs(itype(i,molnum(i)))
+        if (molnum(i).lt.3) then
+                inum=i+nres
+        else
+                inum=i
+        endif
+
         do j=nnt,i-kkk
-          itj=iabs(itype(j,1))
+          itj=iabs(itype(j,molnum(i)))
+        if (molnum(j).lt.3) then
+                jnum=j+nres
+        else
+                jnum=j
+        endif
           if (ipot.ne.4) then
 !           rcomp=sigmaii(iti,itj)+1.0D0
             rcomp=facont*sigmaii(iti,itj)
@@ -58,7 +69,7 @@
           endif
 !         rcomp=6.5D0
 !         print *,'rcomp=',rcomp,' dist=',dist(nres+i,nres+j)
-         if (dist(nres+i,nres+j).lt.rcomp) then
+         if (dist(inum,jnum).lt.rcomp) then
             ncont=ncont+1
             icont(1,ncont)=i
             icont(2,ncont)=j
@@ -300,6 +311,7 @@
       ncont=0
       ees=0.0
       evdw=0.0
+      print *, "nntt,nct",nnt,nct-2
       do 1 i=nnt,nct-2
         if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) goto 1
         xi=c(1,i)
@@ -492,6 +504,8 @@
       enddo
 
       call elecont(lprint,ncont,icont)
+      print *,"after elecont"
+      if (nres_molec(1).eq.0) return
 
 ! finding parallel beta
 !d      write (iout,*) '------- looking for parallel beta -----------'
@@ -742,7 +756,7 @@
        write(12,'(a37)') "DefPropRes 'coil' '! (helix | sheet)'"
        write(12,'(a20)') "XMacStand ribbon.mac"
          
-        
+      if (nres_molec(1).eq.0) return
        write(iout,*) 'UNRES seq:'
        do j=1,nbfrag
         write(iout,*) 'beta ',(bfrag(i,j),i=1,4)
@@ -815,6 +829,7 @@
 !     &             obr,non_conv)
 !        rms=dsqrt(rms)
         call rmsd(rms)
+!        print *,"before contact"
 !elte(iout,*) "rms_nacc before contact"
         call contact(.false.,ncont,icont,co)
         frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
@@ -877,6 +892,7 @@
 !      else
 !      do kkk=1,nperm
       iatom=0
+      print *,nz_start,nz_end,nstart_seq-nstart_sup
       do i=nz_start,nz_end
         iatom=iatom+1
         iti=itype(i,1)
@@ -4180,7 +4196,7 @@
       do i=nnt,nct
         if (itype(i,1).ne.10) then
 !d        print *,'i=',i,' itype=',itype(i,1),' theta=',theta(i+1)  
-          call gen_side(itype(i,1),theta(i+1),alph(i),omeg(i),fail)
+          call gen_side(itype(i,1),theta(i+1),alph(i),omeg(i),fail,1)
         endif
       enddo
       call chainbuild
@@ -4196,7 +4212,7 @@
         alph0=alph(ind_sc)
         omeg0=omeg(ind_sc)
         call gen_side(itype(ind_sc,1),theta(ind_sc+1),alph(ind_sc),&
-             omeg(ind_sc),fail)
+             omeg(ind_sc),fail,1)
         call chainbuild
         call etotal(energia)
 !d      write (iout,'(a,i5,a,i4,2(a,f8.3),2(a,1pe14.5))')