corrrection in energies (no trash)
[unres4.git] / source / wham / conform_compar.f90
index 701e920..8476646 100644 (file)
 
       use calc_data
       use geometry_data, only:c,dc,dc_norm
-      use energy_data, only:itype,maxcont
+      use energy_data, only:itype,maxcont,molnum
 !      implicit none
 !      include 'DIMENSIONS'
 !      include 'DIMENSIONS.ZSCOPT'
 !      include 'COMMON.CALC'
 !      include 'COMMON.CONTPAR'
 !      include 'COMMON.LOCAL'
-      integer :: ist,ien,kkk,iti,itj,itypi,itypj,i1,i2,it1,it2
+      integer :: ist,ien,kkk,iti,itj,itypi,itypj,i1,i2,it1,it2,mnum,mnum2
       real(kind=8) :: csc !el,dist
       real(kind=8),dimension(maxcont) :: cscore,omt1,omt2,omt12,&
           ddsc,ddla,ddlb
-      integer :: ncont
+      integer :: ncont,mhum
       integer,dimension(2,maxcont) :: icont
       real(kind=8) :: u,v,a(3),b(3),dla,dlb
       logical :: lprint
       kkk=3
       if (lprint) then
       do i=1,nres
-        write (iout,110) restyp(itype(i)),i,c(1,i),c(2,i),&
+        mnum=molnum(i)
+        write (iout,110) restyp(itype(i,mnum),mnum),i,c(1,i),c(2,i),&
           c(3,i),dc(1,nres+i),dc(2,nres+i),dc(3,nres+i),&
           dc_norm(1,nres+i),dc_norm(2,nres+i),dc_norm(3,nres+i)
       enddo
       endif
   110 format (a,'(',i3,')',9f8.3)
       do i=ist,ien-kkk
-        iti=iabs(itype(i))
-        if (iti.le.0 .or. iti.gt.ntyp) cycle
+        mnum=molnum(i)
+        iti=iabs(itype(i,mnum))
+        if (iti.le.0 .or. iti.gt.ntyp_molec(mnum)) cycle
         do j=i+kkk,ien
-          itj=iabs(itype(j))
-          if (itj.le.0 .or. itj.gt.ntyp) cycle
+          mnum2=molnum(j)
+          itj=iabs(itype(j,mnum2))
+          if (itj.le.0 .or. itj.gt.ntyp_molec(mnum2)) cycle
           itypi=iti
           itypj=itj
           xj = c(1,nres+j)-c(1,nres+i)    
       if (lprint) then
         write (iout,'(a)') 'Contact map:'
         do i=1,ncont
+          mnum=molnum(i)
           i1=icont(1,i)
           i2=icont(2,i)
-          it1=itype(i1)
-          it2=itype(i2)
+          it1=itype(i1,molnum(i1))
+          it2=itype(i2,molnum(i2))
+!          print *,"CONTACT",i1,mnum,it1,it2
           write (iout,'(i3,2x,a,i4,2x,a,i4,5f8.3,3f10.5)') &
-           i,restyp(it1),i1,restyp(it2),i2,cscore(i),&
+           i,restyp(it1,mnum),i1,restyp(it2,mnum),i2,cscore(i),&
            sc_cutoff(iabs(it1),iabs(it2)),ddsc(i),ddla(i),ddlb(i),&
            omt1(i),omt2(i),omt12(i)
         enddo
 !----------------------------------------------------------------------------
       subroutine contact(lprint,ncont,icont)
 
-      use energy_data, only: nnt,nct,itype,ipot,maxcont,sigma,sigmaii
+      use energy_data, only: nnt,nct,itype,ipot,maxcont,sigma,sigmaii,molnum
 !      include 'DIMENSIONS'
 !      include 'COMMON.IOUNITS'
 !      include 'COMMON.CHAIN'
       logical :: lprint
       integer :: kkk,i,j,i1,i2,it1,it2,iti,itj
       real(kind=8) :: rcomp
+      integer :: mnum,mnum2
       ncont=0
       kkk=3
 !     print *,'nnt=',nnt,' nct=',nct
       do i=nnt+kkk,nct
-        iti=iabs(itype(i))
+        mnum=molnum(i)
+        iti=iabs(itype(i,1))
         do j=nnt,i-kkk
-          itj=iabs(itype(j))
+          mnum2=molnum(j)
+          itj=iabs(itype(j,1))
           if (ipot.ne.4) then
 !           rcomp=sigmaii(iti,itj)+1.0D0
             rcomp=facont*sigmaii(iti,itj)
       if (lprint) then
         write (iout,'(a)') 'Contact map:'
         do i=1,ncont
+           mnum=molnum(i)
           i1=icont(1,i)
           i2=icont(2,i)
-          it1=itype(i1)
-          it2=itype(i2)
+          it1=itype(i1,1)
+          it2=itype(i2,1)
           write (iout,'(i3,2x,a,i4,2x,a,i4)') &
-           i,restyp(it1),i1,restyp(it2),i2
+           i,restyp(it1,mnum),i1,restyp(it2,mnum),i2
         enddo
       endif
       return
       subroutine pept_cont(lprint,ncont,icont)
 
       use geometry_data, only:c
-      use energy_data, only:maxcont,nnt,nct,itype
+      use energy_data, only:maxcont,nnt,nct,itype,molnum
 !      implicit none
 !      include 'DIMENSIONS'
 !      include 'DIMENSIONS.ZSCOPT'
 !      include 'COMMON.FFIELD'
 !      include 'COMMON.NAMES'
       integer :: ncont,icont(2,maxcont)
-      integer :: i,j,k,kkk,i1,i2,it1,it2
+      integer :: i,j,k,kkk,i1,i2,it1,it2,mnum
       logical :: lprint
 !el      real(kind=8) :: dist
       real(kind=8) :: rcomp=5.5d0
       if (lprint) then
         write (iout,'(a)') 'PP contact map:'
         do i=1,ncont
+          mnum=molnum(i)
           i1=icont(1,i)
           i2=icont(2,i)
-          it1=itype(i1)
-          it2=itype(i2)
+          it1=itype(i1,mnum)
+          it2=itype(i2,mnum)
           write (iout,'(i3,2x,a,i4,2x,a,i4)') &
-           i,restyp(it1),i1,restyp(it2),i2
+           i,restyp(it1,mnum),i1,restyp(it2,mnum),i2
         enddo
       endif
       return
       subroutine contacts_between_fragments(lprint,is,ncont,icont,&
          ncont_interfrag,icont_interfrag)
 
-      use energy_data, only:itype,maxcont
+      use energy_data, only:itype,maxcont,molnum
 !      include 'DIMENSIONS'
 !      include 'DIMENSIONS.ZSCOPT'
 !      include 'DIMENSIONS.COMPAR'
       integer :: icont(2,maxcont),ncont_interfrag(mmaxfrag),&
         icont_interfrag(2,maxcont,mmaxfrag)
       logical :: OK1,OK2,lprint
-      integer :: is,ncont,i,j,ind,nc,k,ic1,ic2,l,i1,i2,it1,it2
+      integer :: is,ncont,i,j,ind,nc,k,ic1,ic2,l,i1,i2,it1,it2,mnum
 ! Determine the contacts that occur within a fragment and between fragments.
       do i=1,nfrag(1)
         do j=1,i
           do k=1,ncont_interfrag(ind)
             i1=icont_interfrag(1,k,ind)
             i2=icont_interfrag(2,k,ind)
-            it1=itype(i1)
-            it2=itype(i2)
+            mnum=molnum(i1)
+            it1=itype(i1,mnum)
+            it2=itype(i2,molnum(i2))
             write (iout,'(i3,2x,a,i4,2x,a,i4)') &
-              i,restyp(it1),i1,restyp(it2),i2
+              i,restyp(it1,mnum),i1,restyp(it2,molnum(i2)),i2
           enddo
         enddo
         write (iout,*)
           do k=1,ncont_interfrag(ind)
             i1=icont_interfrag(1,k,ind)
             i2=icont_interfrag(2,k,ind)
-            it1=itype(i1)
-            it2=itype(i2)
+            mnum=molnum(i1)
+            it1=itype(i1,mnum)
+            it2=itype(i2,molnum(i2))
             write (iout,'(i3,2x,a,i4,2x,a,i4)') &
-              i,restyp(it1),i1,restyp(it2),i2
+              i,restyp(it1,mnum),i1,restyp(it2,molnum(i2)),i2
           enddo
         enddo
         enddo
       subroutine elecont(lprint,ncont,icont,ist,ien)
 
       use geometry_data, only:c
-      use energy_data, only:maxcont,rpp,epp,itype,itel,vblinv,vblinv2
+      use energy_data, only:maxcont,rpp,epp,itype,itel,vblinv,vblinv2,molnum
 !      implicit none
 !      include 'DIMENSIONS'
 !      include 'DIMENSIONS.ZSCOPT'
       real(kind=8),dimension(2,2) :: ael6c,ael3c,appc,bppc
       real(kind=8) :: elcutoff=-0.3d0
       real(kind=8) :: elecutoff_14=-0.5d0
-      integer :: ncont,icont(2,maxcont)
+      integer :: ncont,icont(2,maxcont),mnum
       real(kind=8) :: econt(maxcont)
 !
 ! Load the constants of peptide bond - peptide bond interactions.
         do i=1,ncont
           i1=icont(1,i)
           i2=icont(2,i)
-          it1=itype(i1)
-          it2=itype(i2)
+          it1=itype(i1,molnum(i1))
+          it2=itype(i2,molnum(i1))
           write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)') &
-           i,restyp(it1),i1,restyp(it2),i2,econt(i)
+           i,restyp(it1,molnum(i1)),i1,restyp(it2,molnum(i1)),i2,econt(i)
         enddo
       endif
 ! For given residues keep only the contacts with the greatest energy.
         write (iout,*)
         write (iout,*) 'Electrostatic contacts after pruning: '
         do i=1,ncont
+          mnum=molnum(i)
           i1=icont(1,i)
           i2=icont(2,i)
-          it1=itype(i1)
-          it2=itype(i2)
+            mnum=molnum(i1)
+            it1=itype(i1,mnum)
+            it2=itype(i2,molnum(i2))
           write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)') &
-           i,restyp(it1),i1,restyp(it2),i2,econt(i)
+           i,restyp(it1,mnum),i1,restyp(it2,molnum(i2)),i2,econt(i)
         enddo
       endif
       return
 
       use geometry_data, only:cref,nperm
       use control_data, only:symetr
-      use energy_data, only:nnt,nct,itype
+      use energy_data, only:nnt,nct,itype,molnum
 !      implicit none
 !      include 'DIMENSIONS'
 !      include 'DIMENSIONS.ZSCOPT'
                        (cref(3,jl,kkk)-cref(3,il,kkk))**2)
             dij=dist(il,jl)
             qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2)
-            if (itype(il).ne.10 .or. itype(jl).ne.10) then
+            if (itype(il,molnum(il)).ne.10 .or. itype(jl,molnum(jl)).ne.10) then
               nl=nl+1
               d0ijCM=dsqrt( &
                      (cref(1,jl+nres,kkk)-cref(1,il+nres,kkk))**2+ &
             qq = qq+qqij+qqijCM
             if (lprn) then
               write (iout,*) "il",il," jl",jl,&
-               " itype",itype(il),itype(jl)
+               " itype",itype(il,molnum(il)),itype(jl,molnum(jl))
               write (iout,*)"d0ij",d0ij," dij",dij," d0ijCM",d0ijCM,&
                " dijCM",dijCM," qqij",qqij," qqijCM",qqijCM
             endif
                          (cref(3,jl,kkk)-cref(3,il,kkk))**2)
               dij=dist(il,jl)
               qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2)
-              if (itype(il).ne.10 .or. itype(jl).ne.10) then
+              if (itype(il,molnum(il)).ne.10 .or. itype(jl,molnum(jl)).ne.10) then
                 nl=nl+1
                 d0ijCM=dsqrt( &
                        (cref(1,jl+nres,kkk)-cref(1,il+nres,kkk))**2+ &
               qq = qq+qqij+qqijCM
               if (lprn) then
                 write (iout,*) "i",i," j",j," il",il," jl",jl,&
-                 " itype",itype(il),itype(jl)
+                 " itype",itype(il,molnum(il)),itype(jl,molnum(jl))
                 write (iout,*)"d0ij",d0ij," dij",dij," d0ijCM",d0ijCM,&
                  " dijCM",dijCM," qqij",qqij," qqijCM",qqijCM
               endif
                              (cref(3,kl,kkk)-cref(3,il,kkk))**2)
                   dij=dist(il,kl)
                   qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2)
-                  if (itype(il).ne.10 .or. itype(kl).ne.10) then
+                  if (itype(il,molnum(il)).ne.10 .or. itype(kl,molnum(kl)).ne.10) then
                     nl=nl+1
                     d0ijCM=dsqrt( &
                        (cref(1,kl+nres,kkk)-cref(1,il+nres,kkk))**2+ &
                   qq = qq+qqij+qqijCM
                   if (lprn) then
                     write (iout,*) "i",i," j",j," k",k," l",l," il",il,&
-                      " kl",kl," itype",itype(il),itype(kl)
+                      " kl",kl," itype",itype(il,molnum(il)), &
+                         itype(kl,molnum(kl))
                     write (iout,*) " d0ij",d0ij," dij",dij," d0ijCM",&
                     d0ijCM," dijCM",dijCM," qqij",qqij," qqijCM",qqijCM
                   endif
 
       use control_data, only:symetr
       use geometry_data, only:nperm,cref,c
-      use energy_data, only:itype
+      use energy_data, only:itype,molnum
 !      implicit real*8 (a-h,o-z)
 !      include 'DIMENSIONS'
 !      include 'DIMENSIONS.ZSCOPT'
       do kkk=1,nperm
        nnsup=0
        do i=1,nres
-        if (itype(i).ne.ntyp1) then
+        if (itype(i,molnum(i)).ne.ntyp1_molec(molnum(i))) then
           nnsup=nnsup+1
           do j=1,3
             cc(j,nnsup)=c(j,i)
       subroutine define_fragments
 
       use geometry_data, only:rad2deg
-      use energy_data, only:itype
+      use energy_data, only:itype,molnum
       use compare_data, only:nhfrag,nbfrag,bfrag,hfrag
 !      implicit real*8 (a-h,o-z)
 !      include 'DIMENSIONS'
 !      include 'COMMON.INTERACT'
 !      include 'COMMON.NAMES'
       integer :: nstrand,istrand(2,nres/2)
-      integer :: nhairp,ihairp(2,nres/5) 
+      integer :: nhairp,ihairp(2,nres/5),mnum
       character(len=16) :: strstr(4)=reshape((/'helix','hairpin',&
                           'strand','strand pair'/),shape(strstr))
       integer :: j,i,ii,i1,i2,i3,i4,it1,it2,it3,it4
       write (iout,*) "The following primary fragments were found:"
       write (iout,*) "Helices:",nhfrag
       do i=1,nhfrag
+        mnum=molnum(i)
         i1=ifrag(1,1,i)
         i2=ifrag(2,1,i)
-        it1=itype(i1)
-        it2=itype(i2)
+        it1=itype(i1,mnum)
+        it2=itype(i2,mnum)
         write (iout,'(i3,2x,a,i4,2x,a,i4)') &
-             i,restyp(it1),i1,restyp(it2),i2
+             i,restyp(it1,mnum),i1,restyp(it2,mnum),i2
       enddo
       write (iout,*) "Hairpins:",nhairp
       do i=nhfrag+1,nhfrag+nhairp
         i1=ifrag(1,1,i)
         i2=ifrag(2,1,i)
-        it1=itype(i1)
-        it2=itype(i2)
+        it1=itype(i1,molnum(i1))
+        it2=itype(i2,molnum(i2))
         write (iout,'(i3,2x,a,i4,2x,a,i4,2x)') &
-             i,restyp(it1),i1,restyp(it2),i2
+             i,restyp(it1,molnum(i1)),i1,restyp(it2,molnum(i2)),i2
       enddo
       write (iout,*) "Far strand pairs:",nbfrag
       do i=nhfrag+nhairp+1,nhfrag+nhairp+nbfrag
         i1=ifrag(1,1,i)
         i2=ifrag(2,1,i)
-        it1=itype(i1)
-        it2=itype(i2)
+        it1=itype(i1,molnum(i1))
+        it2=itype(i2,molnum(i1))
         i3=ifrag(1,2,i)
         i4=ifrag(2,2,i)
-        it3=itype(i3)
-        it4=itype(i4)
+        it3=itype(i3,molnum(i3))
+        it4=itype(i4,molnum(i4))
         write (iout,'(i3,2x,a,i4,2x,a,i4," and ",a,i4,2x,a,i4)') &
-             i,restyp(it1),i1,restyp(it2),i2,&
-               restyp(it3),i3,restyp(it4),i4
+             i,restyp(it1,molnum(i1)),i1,restyp(it2,molnum(i2)),i2,&
+               restyp(it3,molnum(i3)),i3,restyp(it4,molnum(i4)),i4
       enddo
       write (iout,*) "Strands:",nstrand
       do i=nhfrag+nhairp+nbfrag+1,nfrag(1)
+        mnum=molnum(i)
         i1=ifrag(1,1,i)
         i2=ifrag(2,1,i)
-        it1=itype(i1)
-        it2=itype(i2)
+        it1=itype(i1,mnum)
+        it2=itype(i2,mnum)
         write (iout,'(i3,2x,a,i4,2x,a,i4)') &
-             i,restyp(it1),i1,restyp(it2),i2
+             i,restyp(it1,mnum),i1,restyp(it2,mnum),i2
       enddo
       call imysort(nfrag(1),2,maxpiece,ifrag(1,1,1),npiece(1,1),&
         istruct(1),n_shift(1,1,1),ang_cut(1),ang_cut1(1),frac_min(1),&
         nc_fragm(1,1),nc_req_setf(1,1))
       write (iout,*) "Fragments after sorting:"
       do i=1,nfrag(1)
+        mnum=molnum(i)
         i1=ifrag(1,1,i)
         i2=ifrag(2,1,i)
-        it1=itype(i1)
-        it2=itype(i2)
+        it1=itype(i1,mnum)
+        it2=itype(i2,mnum)
         write (iout,'(i3,2x,a,i4,2x,a,i4,$)') &
-             i,restyp(it1),i1,restyp(it2),i2
+             i,restyp(it1,molnum(it1)),i1,restyp(it2,molnum(it2)),i2
         if (npiece(i,1).eq.1) then
           write (iout,'(2x,a)') strstr(istruct(i))
         else
           i1=ifrag(1,2,i)
           i2=ifrag(2,2,i)
-          it1=itype(i1)
-          it2=itype(i2)
+          it1=itype(i1,mnum)
+          it2=itype(i2,mnum)
           write (iout,'(2x,a,i4,2x,a,i4,2x,a)') &
-             restyp(it1),i1,restyp(it2),i2,strstr(istruct(i))
+             restyp(it1,molnum(it1)),i1,restyp(it2,molnum(it2)),i2,strstr(istruct(i))
         endif
       enddo
       return
       subroutine secondary2(lprint,lprint_sec,ncont,icont,isecstr)
 
       use geometry_data, only:anatemp,rad2deg,phi,nstart_sup,nend_sup
-      use energy_data, only:itype,maxcont
+      use energy_data, only:itype,maxcont,molnum
       use compare_data, only:bfrag,hfrag,nbfrag,nhfrag
       use compare, only:freeres
 !      implicit real*8 (a-h,o-z)
 !      include 'COMMON.CHAIN'
 !      include 'COMMON.NAMES'
 !      include 'COMMON.INTERACT'
-      integer :: ncont,icont(2,maxcont),isec(nres,4),nsec(nres),&
-        isecstr(nres)
+      integer :: ncont,icont(2,maxcont),isec(nres,4),nsec(nres+1),&
+        isecstr(nres+1)
       logical :: lprint,lprint_sec,not_done !el,freeres
       integer :: i,j,ii1,jj1,i1,j1,ij,k,ien,ist
-      integer :: nstrand,nbeta,nhelix,iii1,jjj1
+      integer :: nstrand,nbeta,nhelix,iii1,jjj1,mnum
       real(kind=8) :: p1,p2
 !rel      external freeres
       character(len=1) :: csec(0:2)=reshape((/'-','E','H'/),shape(csec))
         write (iout,*)
         write (iout,*) "Secondary structure"
         do i=1,nres,80
+          mnum=molnum(i)
           ist=i
           ien=min0(i+79,nres)
           write (iout,*)
           write (iout,'(8(7x,i3))') (k,k=ist+9,ien,10)
-          write (iout,'(80a1)') (onelet(itype(k)),k=ist,ien) 
+          write (iout,'(80a1)') (onelet(itype(k,mnum)),k=ist,ien) 
           write (iout,'(80a1)') (csec(isecstr(k)),k=ist,ien) 
         enddo 
         write (iout,*)