source/wham/conform_compar.F90

   1       module conform_compar
   2 !-----------------------------------------------------------------------------
   3       use names
   4       use io_units
   5       use geometry_data, only:nres
   6       use math, only:pinorm
   7       use geometry, only:dist
   8       use regularize_, only:fitsq
   9 !
  10       use wham_data
  11 #ifndef CLUSTER
  12       use w_compar_data
  13 #endif
  14 #ifdef MPI
  15       use MPI_data
  16 !      include "COMMON.MPI"
  17 #endif
  18       implicit none
  19 !-----------------------------------------------------------------------------
  20 !
  21 !
  22 !-----------------------------------------------------------------------------
  23       contains
  24 #ifndef CLUSTER
  25 !-----------------------------------------------------------------------------
  26 ! conf_compar.F
  27 !-----------------------------------------------------------------------------
  28       subroutine conf_compar(jcon)
  29 !      implicit real*8 (a-h,o-z)
  30       use energy_data, only:icont,ncont,nnt,nct,maxcont!,&
  31 !                      nsccont_frag_ref,isccont_frag_ref
  32 #ifdef MPI
  33       include "mpif.h"
  34 #endif
  35 !      include 'DIMENSIONS'
  36 !      include 'DIMENSIONS.ZSCOPT'
  37 !      include 'DIMENSIONS.COMPAR'
  38 !      include 'DIMENSIONS.FREE'
  39 !      include 'COMMON.CONTROL'
  40 !      include 'COMMON.IOUNITS'
  41 !      include 'COMMON.COMPAR'
  42 !      include 'COMMON.CHAIN'
  43 !      include 'COMMON.INTERACT'
  44 !      include 'COMMON.VAR'
  45 !      include 'COMMON.PEPTCONT'
  46 !      include 'COMMON.CONTACTS1'
  47 !      include 'COMMON.HEADER'
  48 !      include 'COMMON.FREE'
  49 !      include 'COMMON.ENERGIES'
  50 !#ifdef MPI
  51 !      include 'COMMON.MPI'
  52 !#endif
  53 !      integer ilen
  54 !      external ilen
  55       logical :: lprn,print_class
  56       integer :: ncont_frag(mmaxfrag),&
  57        icont_frag(2,maxcont,mmaxfrag),ncontsc,&
  58        icontsc(1,maxcont),nsccont_frag(mmaxfrag),&
  59        isccont_frag(2,maxcont,mmaxfrag)
  60       integer :: isecstr(nres)
  61       integer :: itemp(maxfrag)
  62       character(len=4) :: liczba
  63       real(kind=8) :: Epot,rms
  64       integer :: jcon,i,j,ind,ncnat,nsec_match,ishift,ishif1,ishif2,&
  65                  nc_match,ncon_match,iclass_rms,ishifft_rms,ishiff,ishif
  66       integer :: k,kk,iclass_con,iscor,ik,ishifft_con,idig,iex,im
  67 !      print *,"Enter conf_compar",jcon
  68        lprn=.false.
  69        print_class=.false.
  70         write (iout,*) "before anything"
  71         call flush(iout)
  72
  73 !      call angnorm12(rmsang)
  74 ! Level 1: check secondary and supersecondary structure
  75 !      call elecont(lprn,ncont,icont,nnt,nct)
  76       if (lprn) then
  77         write (iout,*) "elecont finished"
  78         call flush(iout)
  79       endif
  80       call secondary2(lprn,.false.,ncont,icont,isecstr)
  81       if (lprn) then
  82         write (iout,*) "secondary2 finished"
  83         call flush(iout)
  84       endif
  85       call contact(lprn,ncontsc,icontsc,nnt,nct)
  86       if (lprn) then
  87          write(iout,*) "Assigning electrostatic contacts"
  88          call flush(iout)
  89       endif
  90       call contacts_between_fragments(lprn,3,ncont,icont,ncont_frag,&
  91          icont_frag)
  92       if (lprn) then
  93         write(iout,*) "Assigning sidechain contacts"
  94         call flush(iout)
  95       endif
  96       call contacts_between_fragments(lprn,3,ncontsc,icontsc,&
  97          nsccont_frag,isccont_frag)
  98       if (lprn) then
  99         write(iout,*) "--> After contacts_between_fragments"
 100         call flush(iout)
 101       endif
 102       do i=1,nlevel
 103         do j=1,isnfrag(nlevel+1)
 104           iclass(j,i)=0
 105         enddo
 106       enddo
 107       do j=1,nfrag(1)
 108         ind = icant(j,j)
 109         if (lprn) then
 110           write (iout,'(80(1h=))')
 111           write (iout,*) "Level",1," fragment",j
 112           write (iout,'(80(1h=))')
 113         endif
 114         call flush(iout)
 115         rmsfrag(j,1)=rmscalc(0,1,j,jcon,lprn)
 116 ! Compare electrostatic contacts in the current conf with that in the native
 117 ! structure.
 118         if (lprn) write (iout,*) &
 119           "Comparing electrostatic contact map and local structure"
 120         call flush(iout)
 121         ncnat=ncont_frag_ref(ind)
 122 !        write (iout,*) "before match_contact:",nc_fragm(j,1),
 123 !     &   nc_req_setf(j,1)
 124 !        call flush(iout)
 125         call match_secondary(j,isecstr,nsec_match,lprn)
 126         if (lprn) write (iout,*) "Fragment",j," nsec_match",&
 127           nsec_match," length",len_frag(j,1)," min_len",&
 128           frac_sec*len_frag(j,1)
 129         if (nsec_match.lt.frac_sec*len_frag(j,1)) then
 130           iclass(j,1)=0
 131           if (lprn) write (iout,*) "Fragment",j,&
 132             " has incorrect secondary structure"
 133         else
 134           iclass(j,1)=1
 135           if (lprn) write (iout,*) "Fragment",j,&
 136             " has correct secondary structure"
 137         endif
 138         if (ielecont(j,1).gt.0) then
 139           call match_contact(ishif1,ishif2,nc_match,ncon_match,&
 140             ncont_frag_ref(ind),icont_frag_ref(1,1,ind),&
 141             ncont_frag(ind),icont_frag(1,1,ind),&
 142             j,n_shift(1,j,1),n_shift(2,j,1),nc_fragm(j,1),&
 143             nc_req_setf(j,1),istruct(j),.true.,lprn)
 144         else if (isccont(j,1).gt.0) then
 145           call match_contact(ishif1,ishif2,nc_match,ncon_match,&
 146             nsccont_frag_ref(ind),isccont_frag_ref(1,1,ind),&
 147             nsccont_frag(ind),isccont_frag(1,1,ind),&
 148             j,n_shift(1,j,1),n_shift(2,j,1),nc_fragm(j,1),&
 149             nc_req_setf(j,1),istruct(j),.true.,lprn)
 150         else if (iloc(j).gt.0) then
 151 !          write (iout,*) "n_shif",n_shift(1,j,1),n_shift(2,j,1)
 152           call match_contact(ishif1,ishif2,nc_match,ncon_match,&
 153             0,icont_frag_ref(1,1,ind),&
 154             ncont_frag(ind),icont_frag(1,1,ind),&
 155             j,n_shift(1,j,1),n_shift(2,j,1),nc_fragm(j,1),&
 156             0,istruct(j),.true.,lprn)
 157 !          write (iout,*) "n_shif",n_shift(1,j,1),n_shift(2,j,1)
 158         else
 159           ishif=0
 160           nc_match=1
 161         endif
 162         if (lprn) write (iout,*) "ishif1",ishif1," ishif2",ishif2
 163         ishif=ishif1
 164         qfrag(j,1)=qwolynes(1,j)
 165         if (iabs(ishif2).gt.iabs(ishif1)) ishif=ishif2
 166         if (lprn) write (iout,*) "ishift",ishif," nc_match",nc_match
 167 !        write (iout,*) "j",j," ishif",ishif," rms",rmsfrag(j,1)
 168         if (irms(j,1).gt.0) then
 169           if (rmsfrag(j,1).le.rmscutfrag(1,j,1)) then
 170             iclass_rms=2
 171             ishifft_rms=0
 172           else
 173             ishiff=0
 174             rms=1.0d2
 175             iclass_rms=0
 176             do while (rms.gt.rmscutfrag(1,j,1) .and. &
 177                ishiff.lt.n_shift(1,j,1))
 178               ishiff=ishiff+1
 179               rms=rmscalc(-ishiff,1,j,jcon,lprn)
 180 !              write(iout,*)"jcon,i,j,ishiff",jcon,i,j,-ishiff,
 181 !     &          " rms",rms," rmscut",rmscutfrag(1,j,1)
 182               if (lprn) write (iout,*) "rms",rmsfrag(j,1)
 183               if (rms.gt.rmscutfrag(1,j,1)) then
 184                 rms=rmscalc(ishiff,1,j,jcon,lprn)
 185 !                write (iout,*) "jcon,1,j,ishiff",jcon,1,j,ishiff,
 186 !     &           " rms",rms
 187               endif
 188               if (lprn) write (iout,*) "rms",rmsfrag(j,1)
 189             enddo
 190 !            write (iout,*) "After loop: rms",rms,
 191 !     &        " rmscut",rmscutfrag(1,j,1)
 192 !            write (iout,*) "iclass_rms",iclass_rms
 193             if (rms.le.rmscutfrag(1,j,1)) then
 194               ishifft_rms=ishiff
 195               rmsfrag(j,1)=rms
 196               iclass_rms=1
 197             endif
 198 !            write (iout,*) "iclass_rms",iclass_rms
 199           endif
 200 !          write (iout,*) "ishif",ishif
 201           if (iabs(ishifft_rms).gt.iabs(ishif)) ishif=ishifft_rms
 202         else
 203           iclass_rms=1
 204         endif
 205 !        write (iout,*) "ishif",ishif," iclass",iclass(j,1),
 206 !     &    " iclass_rms",iclass_rms
 207         if (nc_match.gt.0 .and. iclass_rms.gt.0) then
 208           if (ishif.eq.0) then
 209             iclass(j,1)=iclass(j,1)+6
 210           else
 211             iclass(j,1)=iclass(j,1)+2
 212           endif
 213         endif
 214         ncont_nat(1,j,1)=nc_match
 215         ncont_nat(2,j,1)=ncon_match
 216         ishifft(j,1)=ishif
 217 !        write (iout,*) "iclass",iclass(j,1)
 218       enddo
 219 ! Next levels: Check arrangements of elementary fragments.
 220       do i=2,nlevel
 221         do j=1,nfrag(i)
 222         if (i .eq. 2) ind = icant(ipiece(1,j,i),ipiece(2,j,i))
 223         if (lprn) then
 224             write (iout,'(80(1h=))')
 225             write (iout,*) "Level",i," fragment",j
 226             write (iout,'(80(1h=))')
 227         endif
 228 ! If an elementary fragment doesn't exist, don't check higher hierarchy levels.
 229         do k=1,npiece(j,i)
 230           ik=ipiece(k,j,i)
 231           if (iclass(ik,1).eq.0) then
 232             iclass(j,i)=0
 233             goto 12
 234           endif
 235         enddo
 236         if (i.eq.2 .and. ielecont(j,i).gt.0) then
 237           iclass_con=0
 238           ishifft_con=0
 239           if (lprn) write (iout,*) &
 240            "Comparing electrostatic contact map: fragments",&
 241             ipiece(1,j,i),ipiece(2,j,i)," ind",ind
 242           call match_contact(ishif1,ishif2,nc_match,ncon_match,&
 243            ncont_frag_ref(ind),icont_frag_ref(1,1,ind),&
 244            ncont_frag(ind),icont_frag(1,1,ind),&
 245            j,n_shift(1,j,i),n_shift(2,j,i),nc_fragm(j,i),&
 246            nc_req_setf(j,i),2,.false.,lprn)
 247           ishif=ishif1
 248           if (iabs(ishif2).gt.iabs(ishif1)) ishif=ishif2
 249           if (nc_match.gt.0) then
 250             if (ishif.eq.0) then
 251               iclass_con=2
 252             else
 253               iclass_con=1
 254             endif
 255           endif
 256           ncont_nat(1,j,i)=nc_match
 257           ncont_nat(2,j,i)=ncon_match
 258           ishifft_con=ishif
 259         else if (i.eq.2 .and. isccont(j,i).gt.0) then
 260           iclass_con=0
 261           ishifft_con=0
 262           if (lprn) write (iout,*) &
 263            "Comparing sidechain contact map: fragments",&
 264            ipiece(1,j,i),ipiece(2,j,i)," ind",ind
 265           call match_contact(ishif1,ishif2,nc_match,ncon_match,&
 266            nsccont_frag_ref(ind),isccont_frag_ref(1,1,ind),&
 267            nsccont_frag(ind),isccont_frag(1,1,ind),&
 268            j,n_shift(1,j,i),n_shift(2,j,i),nc_fragm(j,i),&
 269            nc_req_setf(j,i),2,.false.,lprn)
 270           ishif=ishif1
 271           if (iabs(ishif2).gt.iabs(ishif1)) ishif=ishif2
 272           if (nc_match.gt.0) then
 273             if (ishif.eq.0) then
 274               iclass_con=2
 275             else
 276               iclass_con=1
 277             endif
 278           endif
 279           ncont_nat(1,j,i)=nc_match
 280           ncont_nat(2,j,i)=ncon_match
 281           ishifft_con=ishif
 282         else if (i.eq.2) then
 283           iclass_con=2
 284           ishifft_con=0
 285         endif
 286         if (i.eq.2) qfrag(j,2)=qwolynes(2,j)
 287         if (lprn) write (iout,*) &
 288           "Comparing rms: fragments",&
 289            (ipiece(k,j,i),k=1,npiece(j,i))
 290         rmsfrag(j,i)=rmscalc(0,i,j,jcon,lprn)
 291         if (irms(j,i).gt.0) then
 292           iclass_rms=0
 293           ishifft_rms=0
 294           if (lprn) write (iout,*) "rms",rmsfrag(j,i)
 295 !          write (iout,*) "i",i," j",j," rmsfrag",rmsfrag(j,i),
 296 !     &     " rmscutfrag",rmscutfrag(1,j,i)
 297           if (rmsfrag(j,i).le.rmscutfrag(1,j,i)) then
 298             iclass_rms=2
 299             ishifft_rms=0
 300           else
 301             ishif=0
 302             rms=1.0d2
 303             do while (rms.gt.rmscutfrag(1,j,i) .and. &
 304                ishif.lt.n_shift(1,j,i))
 305               ishif=ishif+1
 306               rms=rmscalc(-ishif,i,j,jcon,lprn)
 307 !              print *,"jcon,i,j,ishif",jcon,i,j,-ishif," rms",rms
 308               if (lprn) write (iout,*) "rms",rmsfrag(j,i)
 309               if (rms.gt.rmscutfrag(1,j,i)) then
 310                 rms=rmscalc(ishif,i,j,jcon,lprn)
 311 !                print *,"jcon,i,j,ishif",jcon,i,j,ishif," rms",rms
 312               endif
 313               if (lprn) write (iout,*) "rms",rms
 314             enddo
 315             if (rms.le.rmscutfrag(1,j,i)) then
 316               ishifft_rms=ishif
 317               rmsfrag(j,i)=rms
 318               iclass_rms=1
 319             endif
 320           endif
 321         endif
 322         if (irms(j,i).eq.0 .and. ielecont(j,i).eq.0 .and. &
 323           isccont(j,i).eq.0 ) then
 324           write (iout,*) "Error: no measure of comparison specified:",&
 325             " level",i," part",j
 326           stop
 327         endif
 328         if (lprn) &
 329         write (iout,*) "iclass_con",iclass_con," iclass_rms",iclass_rms
 330         if (i.eq.2) then
 331           iclass(j,i) = min0(iclass_con,iclass_rms)
 332           if (iabs(ishifft_rms).gt.iabs(ishifft_con)) then
 333             ishifft(j,i)=ishifft_rms
 334           else
 335             ishifft(j,i)=ishifft_con
 336           endif
 337         else if (i.gt.2) then
 338           iclass(j,i) = iclass_rms
 339           ishifft(j,i)= ishifft_rms
 340         endif
 341    12   continue
 342         enddo
 343       enddo
 344       rms_nat=rmsnat(jcon)
 345       qnat=qwolynes(0,0)
 346 ! Compute the structural class
 347       iscor=0
 348       IF (.NOT. BINARY) THEN
 349       do i=1,nlevel
 350         IF (I.EQ.1) THEN
 351         do j=1,nfrag(i)
 352           itemp(j)=iclass(j,i)
 353         enddo
 354         do kk=-1,1
 355           do j=1,nfrag(i)
 356             idig = 2*isnfrag(nlevel+1)-2*isnfrag(i)-kk*nfrag(i)-j
 357             iex = 2**idig
 358             im=mod(itemp(j),2)
 359             itemp(j)=itemp(j)/2
 360 !            write (iout,*) "i",i," j",j," idig",idig," iex",iex,
 361 !     &        " iclass",iclass(j,i)," im",im
 362             iscor=iscor+im*iex
 363           enddo
 364         enddo
 365         ELSE
 366         do j=1,nfrag(i)
 367           idig = 2*isnfrag(nlevel+1)-2*isnfrag(i)-j
 368           iex = 2**idig
 369           if (iclass(j,i).gt.0) then
 370             im=1
 371           else
 372             im=0
 373           endif
 374 !          write (iout,*) "i",i," j",j," idig",idig," iex",iex,
 375 !     &      " iclass",iclass(j,i)," im",im
 376           iscor=iscor+im*iex
 377         enddo
 378         do j=1,nfrag(i)
 379           idig = 2*isnfrag(nlevel+1)-2*isnfrag(i)-nfrag(i)-j
 380           iex = 2**idig
 381           if (iclass(j,i).gt.1) then
 382             im=1
 383           else
 384             im=0
 385           endif
 386 !          write (iout,*) "i",i," j",j," idig",idig," iex",iex,
 387 !     &      " iclass",iclass(j,i)," im",im
 388           iscor=iscor+im*iex
 389         enddo
 390         ENDIF
 391       enddo
 392       iscore=iscor
 393       ENDIF
 394       if (print_class) then
 395 #ifdef MPI
 396           write(istat,'(i6,$)') jcon+indstart(me)-1
 397           write (istat,'(f10.2,$)') (potE(jcon,k),k=1,nParmSet),&
 398            -entfac(jcon)
 399 #else
 400           write(istat,'(i6,$)') jcon
 401           write (istat,'(f10.2,$)') (potE(jcon,k),k=1,nParmSet),&
 402             -entfac(jcon)
 403 #endif
 404           write (istat,'(f8.3,2f6.3,$)') &
 405             rms_nat,qnat,rmsang/(nres-3)
 406           do j=1,nlevel
 407             write(istat,'(1x,$,20(i3,$))') &
 408               (ncont_nat(1,k,j),k=1,nfrag(j))
 409             if (j.lt.3) then
 410               write(istat,'(1x,$,20(f5.1,f5.2$))') &
 411                 (rmsfrag(k,j),qfrag(k,j),k=1,nfrag(j))
 412             else
 413               write(istat,'(1x,$,20(f5.1$))') &
 414                 (rmsfrag(k,j),k=1,nfrag(j))
 415             endif
 416             write(istat,'(1x,$,20(i1,$))') &
 417               (iclass(k,j),k=1,nfrag(j))
 418           enddo
 419           if (binary) then
 420             write (istat,'("  ",$)')
 421             do j=1,nlevel
 422               write (istat,'(100(i1,$))')(iclass(k,j),&
 423                  k=1,nfrag(j))
 424               if (j.lt.nlevel) write(iout,'(".",$)')
 425             enddo
 426             write (istat,*)
 427           else
 428             write (istat,'(i10)') iscore
 429           endif
 430       endif
 431       RETURN
 432       END subroutine conf_compar
 433 !-----------------------------------------------------------------------------
 434 ! angnorm.f
 435 !-----------------------------------------------------------------------------
 436       subroutine add_angpair(ici,icj,nang_pair,iang_pair)
 437
 438 !      implicit real*8 (a-h,o-z)
 439 !      include 'DIMENSIONS'
 440 !      include 'COMMON.IOUNITS'
 441 !      include 'COMMON.CHAIN'
 442       integer :: ici,icj,nang_pair,iang_pair(2,nres)
 443       integer :: i,ian1,ian2
 444 !      write (iout,*) "add_angpair: ici",ici," icj",icj,
 445 !     &  " nang_pair",nang_pair
 446       ian1=ici+2
 447       if (ian1.lt.4 .or. ian1.gt.nres) return
 448       ian2=icj+2
 449 !      write (iout,*) "ian1",ian1," ian2",ian2
 450       if (ian2.lt.4 .or. ian2.gt.nres) return
 451       do i=1,nang_pair
 452         if (ian1.eq.iang_pair(1,i) .and. ian2.eq.iang_pair(2,i)) return
 453       enddo
 454       nang_pair=nang_pair+1
 455       iang_pair(1,nang_pair)=ian1
 456       iang_pair(2,nang_pair)=ian2
 457       return
 458       end subroutine add_angpair
 459 !-------------------------------------------------------------------------
 460       subroutine angnorm(jfrag,ishif1,ishif2,diffang_max,angn,fract,lprn)
 461
 462       use geometry_data, only:nstart_sup,nend_sup,phi,theta,&
 463                               rad2deg,dwapi
 464 !      implicit real*8 (a-h,o-z)
 465 !      include 'DIMENSIONS'
 466 !      include 'DIMENSIONS.ZSCOPT'
 467 !      include 'DIMENSIONS.COMPAR'
 468 !      include 'COMMON.IOUNITS'
 469 !      include 'COMMON.VAR'
 470 !      include 'COMMON.COMPAR'
 471 !      include 'COMMON.CHAIN'
 472 !      include 'COMMON.GEO'
 473       real(kind=8) :: pinorm,deltang
 474       logical :: lprn
 475       integer :: jfrag,ishif1,ishif2,nn,npart,nn4,nne
 476       real(kind=8) :: diffang_max,angn,fract,ff
 477       integer :: i,j,nbeg,nend,ll,longest
 478       if (lprn) write (iout,'(80(1h*))')
 479       angn=0.0d0
 480       nn = 0
 481       fract = 1.0d0
 482       npart = npiece(jfrag,1)
 483       nn4 = nstart_sup+3
 484       nne = min0(nend_sup,nres)
 485       if (lprn) write (iout,*) "nn4",nn4," nne",nne
 486       do i=1,npart
 487         nbeg = ifrag(1,i,jfrag) + 3 - ishif1
 488         if (nbeg.lt.nn4) nbeg=nn4
 489         nend = ifrag(2,i,jfrag) + 1 - ishif2
 490         if (nend.gt.nne) nend=nne
 491         if (nend.ge.nbeg) then
 492         nn = nn + nend - nbeg + 1
 493         if (lprn) write (iout,*) "i=",i," nbeg",nbeg," nend",nend,&
 494           " nn",nn," ishift1",ishif1," ishift2",ishif2
 495         if (lprn) write (iout,*) "angles"
 496         longest=0
 497         ll = 0
 498         do j=nbeg,nend
 499             write (iout,*) "ishif1",ishif1,j
 500 !          deltang = pinorm(phi(j)-phi_ref(j+ishif1))
 501           deltang=spherang(phi_ref(j+ishif1),theta_ref(j-1+ishif1),&
 502             theta_ref(j+ishif1),phi(j),theta(j-1),theta(j))
 503           if (dabs(deltang).gt.diffang_max) then
 504             if (ll.gt.longest) longest = ll
 505             ll = 0
 506           else
 507             ll=ll+1
 508           endif
 509           if (ll.gt.longest) longest = ll
 510           if (lprn) write (iout,'(i5,3f10.5)')j,rad2deg*phi(j),&
 511            rad2deg*phi_ref(j+ishif1),rad2deg*deltang
 512           angn=angn+dabs(deltang)
 513         enddo
 514         longest=longest+3
 515         ff = dfloat(longest)/dfloat(nend - nbeg + 4)
 516         if (lprn) write (iout,*)"segment",i," longest fragment within",&
 517           diffang_max*rad2deg,":",longest," fraction",ff
 518         if (ff.lt.fract) fract = ff
 519         endif
 520       enddo
 521       if (nn.gt.0) then
 522         angn = angn/nn
 523       else
 524         angn = dwapi
 525       endif
 526       if (lprn) write (iout,*) "nn",nn," norm",rad2deg*angn,&
 527         " fract",fract
 528       return
 529       end subroutine angnorm
 530 !-------------------------------------------------------------------------
 531       subroutine angnorm2(jfrag,ishif1,ishif2,ncont,icont,lprn,&
 532         diffang_max,anorm,fract)
 533
 534       use geometry_data, only:nstart_sup,nend_sup,phi,theta,&
 535                               rad2deg
 536 !      implicit real*8 (a-h,o-z)
 537 !      include 'DIMENSIONS'
 538 !      include 'DIMENSIONS.ZSCOPT'
 539 !      include 'DIMENSIONS.COMPAR'
 540 !      include 'COMMON.IOUNITS'
 541 !      include 'COMMON.VAR'
 542 !      include 'COMMON.COMPAR'
 543 !      include 'COMMON.CHAIN'
 544 !      include 'COMMON.GEO'
 545       integer :: ncont,icont(2,ncont),longest
 546       real(kind=8) :: anorm,diffang_max,fract
 547       integer :: npiece_c,ifrag_c(2,maxpiece),ishift_c(maxpiece)
 548       real(kind=8) :: pinorm
 549       logical :: lprn
 550       integer :: jfrag,ishif1,ishif2
 551       integer :: nn,nn4,nne,npart,i,j,jstart,jend,ic1,ic2,idi,iic
 552       integer :: nbeg,nend,ll
 553       real(kind=8) :: angn,ishifc,deltang,ff
 554
 555       if (lprn) write (iout,'(80(1h*))')
 556 !
 557 ! Determine the segments for which angles will be compared
 558 !
 559       nn4 = nstart_sup+3
 560       nne = min0(nend_sup,nres)
 561       if (lprn) write (iout,*) "nn4",nn4," nne",nne
 562       npart=npiece(jfrag,1)
 563       npiece_c=0
 564       do i=1,npart
 565 !        write (iout,*) "i",i," ifrag",ifrag(1,i,jfrag),ifrag(2,i,jfrag)
 566         if (icont(1,ncont).lt.ifrag(1,i,jfrag) .or. &
 567           icont(1,1).gt.ifrag(2,i,jfrag)) goto 11
 568         jstart=1
 569         do while (jstart.lt.ncont .and. &
 570          icont(1,jstart).lt.ifrag(1,i,jfrag))
 571 !          write (iout,*) "jstart",jstart," icont",icont(1,jstart),
 572 !     &     " ifrag",ifrag(1,i,jfrag)
 573           jstart=jstart+1
 574         enddo
 575 !        write (iout,*) "jstart",jstart," icont",icont(1,jstart),
 576 !     &   " ifrag",ifrag(1,i,jfrag)
 577         if (icont(1,jstart).lt.ifrag(1,i,jfrag)) goto 11
 578         npiece_c=npiece_c+1
 579         ic1=icont(1,jstart)
 580         ifrag_c(1,npiece_c)=icont(1,jstart)
 581         jend=ncont
 582         do while (jend.gt.1 .and. icont(1,jend).gt.ifrag(2,i,jfrag))
 583 !          write (iout,*) "jend",jend," icont",icont(1,jend),
 584 !     &     " ifrag",ifrag(2,i,jfrag)
 585           jend=jend-1
 586         enddo
 587 !        write (iout,*) "jend",jend," icont",icont(1,jend),
 588 !     &   " ifrag",ifrag(2,i,jfrag)
 589         ic2=icont(1,jend)
 590         ifrag_c(2,npiece_c)=icont(1,jend)+1
 591         ishift_c(npiece_c)=ishif1
 592 !        write (iout,*) "1: i",i," jstart:",jstart," jend",jend,
 593 !     &    " ic1",ic1," ic2",ic2,
 594 !     &    " ifrag",ifrag(1,i,jfrag),ifrag(2,i,jfrag)
 595    11   continue
 596         if (ncont.eq.1 .or. icont(2,ncont).gt.icont(2,1)) then
 597           idi=1
 598         else
 599           idi=-1
 600         endif
 601 !        write (iout,*) "idi",idi
 602         if (idi.eq.1) then
 603           if (icont(2,1).gt.ifrag(2,i,jfrag) .or. &
 604             icont(2,ncont).lt.ifrag(1,i,jfrag)) goto 12
 605           jstart=1
 606           do while (jstart.lt.ncont .and. &
 607            icont(2,jstart).lt.ifrag(1,i,jfrag))
 608 !           write (iout,*) "jstart",jstart," icont",icont(2,jstart),
 609 !     &     " ifrag",ifrag(1,i,jfrag)
 610             jstart=jstart+1
 611           enddo
 612 !          write (iout,*) "jstart",jstart," icont",icont(2,jstart),
 613 !     &     " ifrag",ifrag(1,i,jfrag)
 614           if (icont(2,jstart).lt.ifrag(1,i,jfrag)) goto 12
 615           npiece_c=npiece_c+1
 616           ic1=icont(2,jstart)
 617           ifrag_c(2,npiece_c)=icont(2,jstart)+1
 618           jend=ncont
 619           do while (jend.gt.1 .and. icont(2,jend).gt.ifrag(2,i,jfrag))
 620 !            write (iout,*) "jend",jend," icont",icont(2,jend),
 621 !     &     " ifrag",ifrag(2,i,jfrag)
 622             jend=jend-1
 623           enddo
 624 !          write (iout,*) "jend",jend," icont",icont(2,jend),
 625 !     &     " ifrag",ifrag(2,i,jfrag)
 626         else if (idi.eq.-1) then
 627           if (icont(2,ncont).gt.ifrag(2,i,jfrag) .or. &
 628               icont(2,1).lt.ifrag(1,i,jfrag)) goto 12
 629           jstart=ncont
 630           do while (jstart.gt.ncont .and. &
 631            icont(2,jstart).lt.ifrag(1,i,jfrag))
 632 !           write (iout,*) "jstart",jstart," icont",icont(2,jstart),
 633 !     &     " ifrag",ifrag(1,i,jfrag)
 634             jstart=jstart-1
 635           enddo
 636 !          write (iout,*) "jstart",jstart," icont",icont(2,jstart),
 637 !     &     " ifrag",ifrag(1,i,jfrag)
 638           if (icont(2,jstart).lt.ifrag(1,i,jfrag)) goto 12
 639           npiece_c=npiece_c+1
 640           ic1=icont(2,jstart)
 641           ifrag_c(2,npiece_c)=icont(2,jstart)+1
 642           jend=1
 643           do while (jend.lt.ncont .and. &
 644              icont(2,jend).gt.ifrag(2,i,jfrag))
 645 !             write (iout,*) "jend",jend," icont",icont(2,jend),
 646 !     &         " ifrag",ifrag(2,i,jfrag)
 647             jend=jend+1
 648           enddo
 649 !          write (iout,*) "jend",jend," icont",icont(2,jend),
 650 !     &     " ifrag",ifrag(2,i,jfrag)
 651         endif
 652         ic2=icont(2,jend)
 653         if (ic2.lt.ic1) then
 654           iic = ic1
 655           ic1 = ic2
 656           ic2 = iic
 657         endif
 658 !        write (iout,*) "2: i",i," ic1",ic1," ic2",ic2,
 659 !     &    " jstart:",jstart," jend",jend,
 660 !     &    " ifrag",ifrag(1,i,jfrag),ifrag(2,i,jfrag)
 661         ifrag_c(1,npiece_c)=ic1
 662         ifrag_c(2,npiece_c)=ic2+1
 663         ishift_c(npiece_c)=ishif2
 664    12   continue
 665       enddo
 666       if (lprn) then
 667         write (iout,*) "Before merge: npiece_c",npiece_c
 668         do i=1,npiece_c
 669           write (iout,*) ifrag_c(1,i),ifrag_c(2,i),ishift_c(i)
 670         enddo
 671       endif
 672 !
 673 ! Merge overlapping segments (e.g., avoid splitting helices)
 674 !
 675       i=1
 676       do while (i .lt. npiece_c)
 677         if (ishift_c(i).eq.ishift_c(i+1) .and. &
 678            ifrag_c(2,i).gt.ifrag_c(1,i+1)) then
 679            ifrag_c(2,i)=ifrag_c(2,i+1)
 680            do j=i+1,npiece_c
 681              ishift_c(j)=ishift_c(j+1)
 682              ifrag_c(1,j)=ifrag_c(1,j+1)
 683              ifrag_c(2,j)=ifrag_c(2,j+1)
 684            enddo
 685            npiece_c=npiece_c-1
 686         else
 687           i=i+1
 688         endif
 689       enddo
 690       if (lprn) then
 691         write (iout,*) "After merge: npiece_c",npiece_c
 692         do i=1,npiece_c
 693           write (iout,*) ifrag_c(1,i),ifrag_c(2,i),ishift_c(i)
 694         enddo
 695       endif
 696 !
 697 ! Compare angles
 698 !
 699       angn=0.0d0
 700       anorm=0
 701       nn = 0
 702       fract = 1.0d0
 703       npart = npiece_c
 704       do i=1,npart
 705         ishifc=ishift_c(i)
 706         nbeg = ifrag_c(1,i) + 3 - ishifc
 707         if (nbeg.lt.nn4) nbeg=nn4
 708         nend = ifrag_c(2,i)  - ishifc + 1
 709         if (nend.gt.nne) nend=nne
 710         if (nend.ge.nbeg) then
 711         nn = nn + nend - nbeg + 1
 712         if (lprn) write (iout,*) "i=",i," nbeg",nbeg," nend",nend,&
 713           " nn",nn," ishifc",ishifc
 714         if (lprn) write (iout,*) "angles"
 715         longest=0
 716         ll = 0
 717         do j=nbeg,nend-ishifc
 718            if (j.gt.nres) cycle
 719 !          deltang = pinorm(phi(j)-phi_ref(j+ishifc))
 720            write (iout,*) "i=",j," nn",nn," ishifc",ishifc
 721           deltang=spherang(phi_ref(j+ishifc),theta_ref(j-1+ishifc),&
 722             theta_ref(j+ishifc),phi(j),theta(j-1),theta(j))
 723           if (dabs(deltang).gt.diffang_max) then
 724             if (ll.gt.longest) longest = ll
 725             ll = 0
 726           else
 727             ll=ll+1
 728           endif
 729           if (ll.gt.longest) longest = ll
 730           if (lprn) write (iout,'(i5,3f10.5)')j,rad2deg*phi(j),&
 731            rad2deg*phi_ref(j+ishifc),rad2deg*deltang
 732           angn=angn+dabs(deltang)
 733         enddo
 734         longest=longest+3
 735         ff = dfloat(longest)/dfloat(nend - nbeg + 4)
 736         if (lprn) write (iout,*)"segment",i," longest fragment within",&
 737           diffang_max*rad2deg,":",longest," fraction",ff
 738         if (ff.lt.fract) fract = ff
 739         endif
 740       enddo
 741       if (nn.gt.0) anorm = angn/nn
 742       if (lprn) write (iout,*) "nn",nn," norm",anorm," fract:",fract
 743       return
 744       end subroutine angnorm2
 745 !-------------------------------------------------------------------------
 746       real(kind=8) function angnorm1(nang_pair,iang_pair,lprn)
 747
 748       use geometry_data, only:phi,theta,rad2deg
 749 !      implicit real*8 (a-h,o-z)
 750 !      include 'DIMENSIONS'
 751 !      include 'DIMENSIONS.ZSCOPT'
 752 !      include 'DIMENSIONS.COMPAR'
 753 !      include 'COMMON.IOUNITS'
 754 !      include 'COMMON.VAR'
 755 !      include 'COMMON.COMPAR'
 756 !      include 'COMMON.CHAIN'
 757 !      include 'COMMON.GEO'
 758       logical :: lprn
 759       integer :: nang_pair,iang_pair(2,nres)
 760       real(kind=8) :: pinorm
 761       integer :: j,ia1,ia2
 762       real(kind=8) :: angn,deltang
 763       angn=0.0d0
 764       if (lprn) write (iout,'(80(1h*))')
 765       if (lprn) write (iout,*) "nang_pair",nang_pair
 766       if (lprn) write (iout,*) "angles"
 767       do j=1,nang_pair
 768         ia1 = iang_pair(1,j)
 769         ia2 = iang_pair(2,j)
 770 !        deltang = pinorm(phi(ia1)-phi_ref(ia2))
 771          deltang=spherang(phi_ref(ia2),theta_ref(ia2-1),&
 772             theta_ref(ia2),phi(ia2),theta(ia2-1),theta(ia2))
 773         if (lprn) write (iout,'(3i5,3f10.5)')j,ia1,ia2,rad2deg*phi(ia1),&
 774          rad2deg*phi_ref(ia2),rad2deg*deltang
 775         angn=angn+dabs(deltang)
 776       enddo
 777       if (lprn) &
 778       write (iout,*)"nang_pair",nang_pair," angn",rad2deg*angn/nang_pair
 779       angnorm1 = angn/nang_pair
 780       return
 781       end function angnorm1
 782 !------------------------------------------------------------------------------
 783       subroutine angnorm12(diff)
 784
 785       use geometry_data, only:phi,theta,nstart_sup,nend_sup
 786 !      implicit real*8 (a-h,o-z)
 787 !      include 'DIMENSIONS'
 788 !      include 'DIMENSIONS.ZSCOPT'
 789 !      include 'DIMENSIONS.COMPAR'
 790 !      include 'COMMON.IOUNITS'
 791 !      include 'COMMON.VAR'
 792 !      include 'COMMON.COMPAR'
 793 !      include 'COMMON.CHAIN'
 794 !      include 'COMMON.GEO'
 795       real(kind=8) :: pinorm,diff
 796       integer :: nn4,nne,j
 797       diff=0.0d0
 798       nn4 = nstart_sup+3
 799       nne = min0(nend_sup,nres)
 800 !      do j=nn4-1,nne
 801 !        diff = diff+rad2deg*dabs(pinorm(theta(j)-theta_ref(j)))
 802 !      enddo
 803       do j=nn4,nne
 804 !        diff = diff+rad2deg*dabs(pinorm(phi(j)-phi_ref(j)))
 805          diff=diff+spherang(phi_ref(j),theta_ref(j-1),&
 806             theta_ref(j),phi(j),theta(j-1),theta(j))
 807       enddo
 808       return
 809       end subroutine angnorm12
 810 !--------------------------------------------------------------------------------
 811       real(kind=8) function spherang(gam1,theta11,theta12,&
 812          gam2,theta21,theta22)
 813 !      implicit none
 814       use geometry, only:arcos
 815       real(kind=8) :: gam1,theta11,theta12,gam2,theta21,theta22,&
 816         x1,x2,xmed,f1,f2,fmed
 817       real(kind=8) :: tolx=1.0d-4, tolf=1.0d-4
 818       real(kind=8) :: sumcos
 819 !el      real(kind=8) :: pinorm,sumangp !arcos,
 820       integer :: it,maxit=100
 821 ! Calculate the difference of the angles of two superposed 4-redidue fragments
 822 !
 823 !       O      P
 824 !        \    /
 825 !     O'--C--C
 826 !             \
 827 !              P'
 828 !
 829 ! The fragment O'-C-C-P' is rotated by angle fi about the C-C axis
 830 ! to achieve the minimum difference between the O'-C-O and P-C-P angles;
 831 ! the sum of these angles is the difference returned by the function.
 832 !
 833 ! 4/28/04 AL
 834 ! If thetas match, take the difference of gamma and exit.
 835       if (dabs(theta11-theta12).lt.tolx &
 836        .and. dabs(theta21-theta22).lt.tolx) then
 837         spherang=dabs(pinorm(gam2-gam1))
 838         return
 839       endif
 840 ! If the gammas are the same, take the difference of thetas and exit.
 841       x1=0.0d0
 842       x2=0.5d0*pinorm(gam2-gam1)
 843       if (dabs(x2) .lt. tolx) then
 844         spherang=dabs(theta11-theta21)+dabs(theta12-theta22)
 845         return
 846       else if (x2.lt.0.0d0) then
 847         x1=x2
 848         x2=0.0d0
 849       endif
 850 ! Else apply regula falsi method to compute optimum overlap of the terminal Calphas
 851       f1=sumangp(gam1,theta11,theta12,gam2,theta21,theta22,x1)
 852       f2=sumangp(gam1,theta11,theta12,gam2,theta21,theta22,x2)
 853       do it=1,maxit
 854         xmed=x1-f1*(x2-x1)/(f2-f1)
 855         fmed=sumangp(gam1,theta11,theta12,gam2,theta21,theta22,xmed)
 856 !        write (*,*) 'it',it,' xmed ',xmed,' fmed ',fmed
 857         if ( (dabs(xmed-x1).lt.tolx .or. dabs(x2-xmed).lt.tolx) &
 858              .and. dabs(fmed).lt.tolf ) then
 859           x1=xmed
 860           f1=fmed
 861           goto 10
 862         else if ( fmed*f1.lt.0.0d0 ) then
 863           x2=xmed
 864           f2=fmed
 865         else
 866           x1=xmed
 867           f1=fmed
 868         endif
 869       enddo
 870    10 continue
 871       spherang=arcos(dcos(theta11)*dcos(theta12) &
 872        +dsin(theta11)*dsin(theta12)*dcos(x1))+ &
 873        arcos(dcos(theta21)*dcos(theta22)+ &
 874        dsin(theta21)*dsin(theta22)*dcos(gam2-gam1+x1))
 875       return
 876       end function spherang
 877 !--------------------------------------------------------------------------------
 878       real(kind=8) function sumangp(gam1,theta11,theta12,gam2,&
 879        theta21,theta22,fi)
 880 !      implicit none
 881       real(kind=8) :: gam1,theta11,theta12,gam2,theta21,theta22,fi,&
 882        cost11,cost12,cost21,cost22,sint11,sint12,sint21,sint22,cosd1,&
 883        cosd2
 884 ! derivarive of the sum of the difference of the angles of a 4-residue fragment.
 885 !      real(kind=8) :: arcos
 886       cost11=dcos(theta11)
 887       cost12=dcos(theta12)
 888       cost21=dcos(theta21)
 889       cost22=dcos(theta22)
 890       sint11=dsin(theta11)
 891       sint12=dsin(theta12)
 892       sint21=dsin(theta21)
 893       sint22=dsin(theta22)
 894       cosd1=cost11*cost12+sint11*sint12*dcos(fi)
 895       cosd2=cost21*cost22+sint21*sint22*dcos(gam2-gam1+fi)
 896       sumangp=sint11*sint12*dsin(fi)/dsqrt(1.0d0-cosd1*cosd1) &
 897        +sint21*sint22*dsin(gam2-gam1+fi)/dsqrt(1.0d0-cosd2*cosd2)
 898       return
 899       end function sumangp
 900 !-----------------------------------------------------------------------------
 901 ! contact.f
 902 !-----------------------------------------------------------------------------
 903       subroutine contact(lprint,ncont,icont,ist,ien)
 904
 905       use calc_data
 906       use geometry_data, only:c,dc,dc_norm
 907       use energy_data, only:itype,maxcont,molnum
 908 !      implicit none
 909 !      include 'DIMENSIONS'
 910 !      include 'DIMENSIONS.ZSCOPT'
 911 !      include 'COMMON.CONTROL'
 912 !      include 'COMMON.IOUNITS'
 913 !      include 'COMMON.CHAIN'
 914 !      include 'COMMON.INTERACT'
 915 !      include 'COMMON.FFIELD'
 916 !      include 'COMMON.NAMES'
 917 !      include 'COMMON.CALC'
 918 !      include 'COMMON.CONTPAR'
 919 !      include 'COMMON.LOCAL'
 920       integer :: ist,ien,kkk,iti,itj,itypi,itypj,i1,i2,it1,it2,mnum,mnum2
 921       real(kind=8) :: csc !el,dist
 922       real(kind=8),dimension(maxcont) :: cscore,omt1,omt2,omt12,&
 923           ddsc,ddla,ddlb
 924       integer :: ncont,mhum
 925       integer,dimension(2,maxcont) :: icont
 926       real(kind=8) :: u,v,a(3),b(3),dla,dlb
 927       logical :: lprint
 928 !el-------
 929       dla=0.0d0
 930       dlb=0.0d0
 931 !el------
 932       ncont=0
 933       kkk=3
 934       if (lprint) then
 935       do i=1,nres
 936         mnum=molnum(i)
 937         write (iout,110) restyp(itype(i,mnum),mnum),i,c(1,i),c(2,i),&
 938           c(3,i),dc(1,nres+i),dc(2,nres+i),dc(3,nres+i),&
 939           dc_norm(1,nres+i),dc_norm(2,nres+i),dc_norm(3,nres+i)
 940       enddo
 941       endif
 942   110 format (a,'(',i3,')',9f8.3)
 943       do i=ist,ien-kkk
 944         mnum=molnum(i)
 945         iti=iabs(itype(i,mnum))
 946         if (iti.le.0 .or. iti.gt.ntyp_molec(mnum)) cycle
 947         do j=i+kkk,ien
 948           mnum2=molnum(j)
 949           itj=iabs(itype(j,mnum2))
 950           if (itj.le.0 .or. itj.gt.ntyp_molec(mnum2)) cycle
 951           itypi=iti
 952           itypj=itj
 953           xj = c(1,nres+j)-c(1,nres+i)
 954           yj = c(2,nres+j)-c(2,nres+i)
 955           zj = c(3,nres+j)-c(3,nres+i)
 956           dxi = dc_norm(1,nres+i)
 957           dyi = dc_norm(2,nres+i)
 958           dzi = dc_norm(3,nres+i)
 959           dxj = dc_norm(1,nres+j)
 960           dyj = dc_norm(2,nres+j)
 961           dzj = dc_norm(3,nres+j)
 962           do k=1,3
 963             a(k)=dc(k,nres+i)
 964             b(k)=dc(k,nres+j)
 965           enddo
 966 !          write (iout,*) (a(k),k=1,3),(b(k),k=1,3)
 967           if (icomparfunc.eq.1) then
 968             call contfunc(csc,iti,itj)
 969           else if (icomparfunc.eq.2) then
 970             call scdist(csc,iti,itj)
 971           else if (icomparfunc.eq.3 .or. icomparfunc.eq.5) then
 972             csc = dist(nres+i,nres+j)
 973           else if (icomparfunc.eq.4) then
 974             call odlodc(c(1,i),c(1,j),a,b,u,v,dla,dlb,csc)
 975           else
 976             write (*,*) "Error - Unknown sidechain contact function"
 977             write (iout,*) "Error - Unknown sidechain contact function"
 978           endif
 979           if (csc.lt.sc_cutoff(iti,itj)) then
 980 !            write(iout,*) "i",i," j",j," dla",dla,dsc(iti),
 981 !     &      " dlb",dlb,dsc(itj)," csc",csc,sc_cutoff(iti,itj),
 982 !     &      dxi,dyi,dzi,dxi**2+dyi**2+dzi**2,
 983 !     &      dxj,dyj,dzj,dxj**2+dyj**2+dzj**2,om1,om2,om12,
 984 !     &      xj,yj,zj
 985 !            write(iout,*)'egb',itypi,itypj,chi1,chi2,chip1,chip2,
 986 !     &       sig0ij,rij,rrij,om1,om2,om12,chiom1,chiom2,chiom12,
 987 !     &       chipom1,chipom2,chipom12,sig,eps2rt,rij_shift,e2,evdw,
 988 !     &       csc
 989             ncont=ncont+1
 990             cscore(ncont)=csc
 991             icont(1,ncont)=i
 992             icont(2,ncont)=j
 993             omt1(ncont)=om1
 994             omt2(ncont)=om2
 995             omt12(ncont)=om12
 996             ddsc(ncont)=1.0d0/rij
 997             ddla(ncont)=dla
 998             ddlb(ncont)=dlb
 999           endif
1000         enddo
1001       enddo
1002       if (lprint) then
1003         write (iout,'(a)') 'Contact map:'
1004         do i=1,ncont
1005
1006           i1=icont(1,i)
1007           i2=icont(2,i)
1008 !          mnum=molnum(i1)
1009           it1=itype(i1,1)
1010           it2=itype(i2,1)
1011           if (mnum.eq.0) mnum=1
1012           if ((it1.eq.0).or.(it2.eq.0)) then
1013           it1=1
1014           it2=1
1015           endif
1016           print *,"CONTACT",i1,i2,mnum,it1,it2
1017           write (iout,'(i3,2x,a,i4,2x,a,i4,5f8.3,3f10.5)') &
1018            i,restyp(it1,mnum),i1,restyp(it2,mnum),i2,cscore(i),&
1019            sc_cutoff(iabs(it1),iabs(it2)),ddsc(i),ddla(i),ddlb(i),&
1020            omt1(i),omt2(i),omt12(i)
1021         enddo
1022       endif
1023       return
1024       end subroutine contact
1025 #else
1026 !----------------------------------------------------------------------------
1027       subroutine contact(lprint,ncont,icont)
1028
1029       use energy_data, only: nnt,nct,itype,ipot,maxcont,sigma,sigmaii,molnum
1030 !      include 'DIMENSIONS'
1031 !      include 'COMMON.IOUNITS'
1032 !      include 'COMMON.CHAIN'
1033 !      include 'COMMON.INTERACT'
1034 !      include 'COMMON.FFIELD'
1035 !      include 'COMMON.NAMES'
1036       real(kind=8) :: facont=1.569D0  ! facont = (2/(1-sqrt(1-1/4)))**(1/6)
1037       integer :: ncont,icont(2,maxcont)
1038       logical :: lprint
1039       integer :: kkk,i,j,i1,i2,it1,it2,iti,itj
1040       real(kind=8) :: rcomp
1041       integer :: mnum,mnum2
1042       ncont=0
1043       kkk=3
1044 !     print *,'nnt=',nnt,' nct=',nct
1045       do i=nnt+kkk,nct
1046         mnum=molnum(i)
1047         iti=iabs(itype(i,1))
1048         do j=nnt,i-kkk
1049           mnum2=molnum(j)
1050           itj=iabs(itype(j,1))
1051           if (ipot.ne.4) then
1052 !           rcomp=sigmaii(iti,itj)+1.0D0
1053             rcomp=facont*sigmaii(iti,itj)
1054           else
1055 !           rcomp=sigma(iti,itj)+1.0D0
1056             rcomp=facont*sigma(iti,itj)
1057           endif
1058 !         rcomp=6.5D0
1059 !         print *,'rcomp=',rcomp,' dist=',dist(nres+i,nres+j)
1060           if (dist(nres+i,nres+j).lt.rcomp) then
1061             ncont=ncont+1
1062             icont(1,ncont)=i
1063             icont(2,ncont)=j
1064           endif
1065         enddo
1066       enddo
1067       if (lprint) then
1068         write (iout,'(a)') 'Contact map:'
1069         do i=1,ncont
1070            mnum=molnum(i)
1071           i1=icont(1,i)
1072           i2=icont(2,i)
1073           it1=itype(i1,1)
1074           it2=itype(i2,1)
1075           write(iout,*) "test",i1,i2,it1,it2
1076           write (iout,'(i3,2x,a,i4,2x,a,i4)') &
1077            i,restyp(it1,mnum),i1,restyp(it2,mnum),i2
1078         enddo
1079       endif
1080       return
1081       end subroutine contact
1082 #endif
1083 !----------------------------------------------------------------------------
1084       real(kind=8) function contact_fract(ncont,ncont_ref,&
1085                                            icont,icont_ref)
1086
1087       use energy_data, only:maxcont
1088 !      implicit none
1089 !      include 'DIMENSIONS'
1090 !      include 'COMMON.IOUNITS'
1091       integer :: i,j,nmatch
1092       integer :: ncont,ncont_ref
1093       integer,dimension(2,maxcont) :: icont,icont_ref
1094       nmatch=0
1095 !     print *,'ncont=',ncont,' ncont_ref=',ncont_ref
1096 !     write (iout,'(20i4)') (icont_ref(1,i),i=1,ncont_ref)
1097 !     write (iout,'(20i4)') (icont_ref(2,i),i=1,ncont_ref)
1098 !     write (iout,'(20i4)') (icont(1,i),i=1,ncont)
1099 !     write (iout,'(20i4)') (icont(2,i),i=1,ncont)
1100       do i=1,ncont
1101         do j=1,ncont_ref
1102           if (icont(1,i).eq.icont_ref(1,j) .and. &
1103               icont(2,i).eq.icont_ref(2,j)) nmatch=nmatch+1
1104         enddo
1105       enddo
1106 !     print *,' nmatch=',nmatch
1107 !     contact_fract=dfloat(nmatch)/dfloat(max0(ncont,ncont_ref))
1108       contact_fract=dfloat(nmatch)/dfloat(ncont_ref)
1109       return
1110       end function contact_fract
1111 #ifndef CLUSTER
1112 !------------------------------------------------------------------------------
1113       subroutine pept_cont(lprint,ncont,icont)
1114
1115       use geometry_data, only:c
1116       use energy_data, only:maxcont,nnt,nct,itype,molnum
1117 !      implicit none
1118 !      include 'DIMENSIONS'
1119 !      include 'DIMENSIONS.ZSCOPT'
1120 !      include 'COMMON.IOUNITS'
1121 !      include 'COMMON.CHAIN'
1122 !      include 'COMMON.INTERACT'
1123 !      include 'COMMON.FFIELD'
1124 !      include 'COMMON.NAMES'
1125       integer :: ncont,icont(2,maxcont)
1126       integer :: i,j,k,kkk,i1,i2,it1,it2,mnum
1127       logical :: lprint
1128 !el      real(kind=8) :: dist
1129       real(kind=8) :: rcomp=5.5d0
1130       ncont=0
1131       kkk=0
1132       print *,'Entering pept_cont: nnt=',nnt,' nct=',nct
1133       do i=nnt,nct-3
1134         do k=1,3
1135           c(k,2*nres+1)=0.5d0*(c(k,i)+c(k,i+1))
1136         enddo
1137         do j=i+2,nct-1
1138           do k=1,3
1139             c(k,2*nres+2)=0.5d0*(c(k,j)+c(k,j+1))
1140           enddo
1141           if (dist(2*nres+1,2*nres+2).lt.rcomp) then
1142             ncont=ncont+1
1143             icont(1,ncont)=i
1144             icont(2,ncont)=j
1145           endif
1146         enddo
1147       enddo
1148       if (lprint) then
1149         write (iout,'(a)') 'PP contact map:'
1150         do i=1,ncont
1151           mnum=molnum(i)
1152           i1=icont(1,i)
1153           i2=icont(2,i)
1154           it1=itype(i1,mnum)
1155           it2=itype(i2,mnum)
1156           write (iout,'(i3,2x,a,i4,2x,a,i4)') &
1157            i,restyp(it1,mnum),i1,restyp(it2,mnum),i2
1158         enddo
1159       endif
1160       return
1161       end subroutine pept_cont
1162 !-----------------------------------------------------------------------------
1163 ! cont_frag.f
1164 !-----------------------------------------------------------------------------
1165       subroutine contacts_between_fragments(lprint,is,ncont,icont,&
1166          ncont_interfrag,icont_interfrag)
1167
1168       use energy_data, only:itype,maxcont,molnum
1169 !      include 'DIMENSIONS'
1170 !      include 'DIMENSIONS.ZSCOPT'
1171 !      include 'DIMENSIONS.COMPAR'
1172 !      include 'COMMON.INTERACT'
1173 !      include 'COMMON.COMPAR'
1174 !      include 'COMMON.IOUNITS'
1175 !      include 'COMMON.CHAIN'
1176 !      include 'COMMON.NAMES'
1177       integer :: icont(2,maxcont),ncont_interfrag(mmaxfrag),&
1178         icont_interfrag(2,maxcont,mmaxfrag)
1179       logical :: OK1,OK2,lprint
1180       integer :: is,ncont,i,j,ind,nc,k,ic1,ic2,l,i1,i2,it1,it2,mnum
1181 ! Determine the contacts that occur within a fragment and between fragments.
1182       do i=1,nfrag(1)
1183         do j=1,i
1184           ind = icant(i,j)
1185           nc=0
1186 !          write (iout,*) "i",i,(ifrag(1,k,i),ifrag(2,k,i)
1187 !     &      ,k=1,npiece(i,1))
1188 !          write (iout,*) "j",j,(ifrag(1,k,j),ifrag(2,k,j)
1189 !     &      ,k=1,npiece(j,1))
1190 !          write (iout,*) "ncont",ncont
1191           do k=1,ncont
1192             ic1=icont(1,k)
1193             ic2=icont(2,k)
1194             OK1=.false.
1195             l=0
1196             do while (.not.OK1 .and. l.lt.npiece(j,1))
1197               l=l+1
1198               OK1=ic1.ge.ifrag(1,l,j)-is .and. &
1199                ic1.le.ifrag(2,l,j)+is
1200             enddo
1201             OK2=.false.
1202             l=0
1203             do while (.not.OK2 .and. l.lt.npiece(i,1))
1204               l=l+1
1205               OK2=ic2.ge.ifrag(1,l,i)-is .and. &
1206                ic2.le.ifrag(2,l,i)+is
1207             enddo
1208 !            write(iout,*) "k",k," ic1",ic1," ic2",ic2," OK1",OK1,
1209 !     &        " OK2",OK2
1210             if (OK1.and.OK2) then
1211               nc=nc+1
1212               icont_interfrag(1,nc,ind)=ic1
1213               icont_interfrag(2,nc,ind)=ic2
1214 !              write (iout,*) "nc",nc," ic1",ic1," ic2",ic2
1215             endif
1216           enddo
1217           ncont_interfrag(ind)=nc
1218 !          do k=1,ncont_interfrag(ind)
1219 !              i1=icont_interfrag(1,k,ind)
1220 !              i2=icont_interfrag(2,k,ind)
1221 !              it1=itype(i1)
1222 !              it2=itype(i2)
1223 !              write (iout,'(i3,2x,a,i4,2x,a,i4)')
1224 !     &          i,restyp(it1),i1,restyp(it2),i2
1225 !          enddo
1226         enddo
1227       enddo
1228       if (lprint) then
1229         write (iout,*) "Contacts within fragments:"
1230         do i=1,nfrag(1)
1231           write (iout,*) "Fragment",i," (",(ifrag(1,k,i),&
1232            ifrag(2,k,i),k=1,npiece(i,1)),")"
1233           ind=icant(i,i)
1234           do k=1,ncont_interfrag(ind)
1235             i1=icont_interfrag(1,k,ind)
1236             i2=icont_interfrag(2,k,ind)
1237             mnum=molnum(i1)
1238             it1=itype(i1,mnum)
1239             it2=itype(i2,molnum(i2))
1240             write (iout,'(i3,2x,a,i4,2x,a,i4)') &
1241               i,restyp(it1,mnum),i1,restyp(it2,molnum(i2)),i2
1242           enddo
1243         enddo
1244         write (iout,*)
1245         write (iout,*) "Contacts between fragments:"
1246         do i=1,nfrag(1)
1247         do j=1,i-1
1248           ind = icant(i,j)
1249           write (iout,*) "Fragments",i," (",(ifrag(1,k,i),&
1250            ifrag(2,k,i),k=1,npiece(i,1)),") and",j," (",&
1251            (ifrag(1,k,j),ifrag(2,k,j),k=1,npiece(j,1)),")"
1252           write (iout,*) "Number of contacts",&
1253            ncont_interfrag(ind)
1254           ind=icant(i,j)
1255           do k=1,ncont_interfrag(ind)
1256             i1=icont_interfrag(1,k,ind)
1257             i2=icont_interfrag(2,k,ind)
1258             mnum=molnum(i1)
1259             it1=itype(i1,mnum)
1260             it2=itype(i2,molnum(i2))
1261             write (iout,'(i3,2x,a,i4,2x,a,i4)') &
1262               i,restyp(it1,mnum),i1,restyp(it2,molnum(i2)),i2
1263           enddo
1264         enddo
1265         enddo
1266       endif
1267       return
1268       end subroutine contacts_between_fragments
1269 !-----------------------------------------------------------------------------
1270 ! contfunc.f
1271 !-----------------------------------------------------------------------------
1272       subroutine contfunc(cscore,itypi,itypj)
1273 !
1274 ! This subroutine calculates the contact function based on
1275 ! the Gay-Berne potential of interaction.
1276 !
1277       use calc_data
1278 !      implicit real*8 (a-h,o-z)
1279 !      include 'DIMENSIONS'
1280 !      include 'COMMON.CONTPAR'
1281 !      include 'COMMON.CALC'
1282       integer :: expon=6
1283       integer :: itypi,itypj
1284       real(kind=8) :: cscore,sig0ij,rrij,sig,rij_shift,evdw,e2
1285 !
1286       sig0ij=sig_comp(itypi,itypj)
1287       chi1=chi_comp(itypi,itypj)
1288       chi2=chi_comp(itypj,itypi)
1289       chi12=chi1*chi2
1290       chip1=chip_comp(itypi,itypj)
1291       chip2=chip_comp(itypj,itypi)
1292       chip12=chip1*chip2
1293       rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
1294       rij=dsqrt(rrij)
1295 ! Calculate angle-dependent terms of the contact function
1296       erij(1)=xj*rij
1297       erij(2)=yj*rij
1298       erij(3)=zj*rij
1299       om1=dxi*erij(1)+dyi*erij(2)+dzi*erij(3)
1300       om2=dxj*erij(1)+dyj*erij(2)+dzj*erij(3)
1301       om12=dxi*dxj+dyi*dyj+dzi*dzj
1302       chiom12=chi12*om12
1303 !      print *,'egb',itypi,itypj,chi1,chi2,chip1,chip2,
1304 !     &  sig0ij,
1305 !     &  rij,rrij,om1,om2,om12
1306 ! Calculate eps1(om12)
1307       faceps1=1.0D0-om12*chiom12
1308       faceps1_inv=1.0D0/faceps1
1309       eps1=dsqrt(faceps1_inv)
1310 ! Following variable is eps1*deps1/dom12
1311       eps1_om12=faceps1_inv*chiom12
1312 ! Calculate sigma(om1,om2,om12)
1313       om1om2=om1*om2
1314       chiom1=chi1*om1
1315       chiom2=chi2*om2
1316       facsig=om1*chiom1+om2*chiom2-2.0D0*om1om2*chiom12
1317       sigsq=1.0D0-facsig*faceps1_inv
1318 ! Calculate eps2 and its derivatives in om1, om2, and om12.
1319       chipom1=chip1*om1
1320       chipom2=chip2*om2
1321       chipom12=chip12*om12
1322       facp=1.0D0-om12*chipom12
1323       facp_inv=1.0D0/facp
1324       facp1=om1*chipom1+om2*chipom2-2.0D0*om1om2*chipom12
1325 ! Following variable is the square root of eps2
1326       eps2rt=1.0D0-facp1*facp_inv
1327       sigsq=1.0D0/sigsq
1328       sig=sig0ij*dsqrt(sigsq)
1329       rij_shift=1.0D0/rij-sig+sig0ij
1330       if (rij_shift.le.0.0D0) then
1331         evdw=1.0D1
1332         cscore = -dlog(evdw+1.0d-6)
1333         return
1334       endif
1335       rij_shift=1.0D0/rij_shift
1336       e2=(rij_shift*sig0ij)**expon
1337       evdw=dabs(eps1*eps2rt**2*e2)
1338       if (evdw.gt.1.0d1) evdw = 1.0d1
1339       cscore = -dlog(evdw+1.0d-6)
1340       return
1341       end subroutine contfunc
1342 !------------------------------------------------------------------------------
1343       subroutine scdist(cscore,itypi,itypj)
1344 !
1345 ! This subroutine calculates the contact distance
1346 !
1347       use calc_data
1348 !      implicit real*8 (a-h,o-z)
1349 !      include 'DIMENSIONS'
1350 !      include 'COMMON.CONTPAR'
1351 !      include 'COMMON.CALC'
1352       integer :: itypi,itypj
1353       real(kind=8) :: cscore,rrij
1354
1355       chi1=chi_comp(itypi,itypj)
1356       chi2=chi_comp(itypj,itypi)
1357       chi12=chi1*chi2
1358       rrij=xj*xj+yj*yj+zj*zj
1359       rij=dsqrt(rrij)
1360 ! Calculate angle-dependent terms of the contact function
1361       erij(1)=xj/rij
1362       erij(2)=yj/rij
1363       erij(3)=zj/rij
1364       om1=dxi*erij(1)+dyi*erij(2)+dzi*erij(3)
1365       om2=dxj*erij(1)+dyj*erij(2)+dzj*erij(3)
1366       om12=dxi*dxj+dyi*dyj+dzi*dzj
1367       chiom12=chi12*om12
1368       om1om2=om1*om2
1369       chiom1=chi1*om1
1370       chiom2=chi2*om2
1371       cscore=dsqrt(rrij+chi1**2+chi2**2+2*rij*(chiom2-chiom1)-2*chiom12)
1372       return
1373       end subroutine scdist
1374 !------------------------------------------------------------------------------
1375 ! elecont.f
1376 !------------------------------------------------------------------------------
1377       subroutine elecont(lprint,ncont,icont,ist,ien)
1378
1379       use geometry_data, only:c
1380       use energy_data, only:maxcont,rpp,epp,itype,itel,vblinv,vblinv2,molnum
1381 !      implicit none
1382 !      include 'DIMENSIONS'
1383 !      include 'DIMENSIONS.ZSCOPT'
1384 !      include 'DIMENSIONS.COMPAR'
1385 !      include 'COMMON.IOUNITS'
1386 !      include 'COMMON.CHAIN'
1387 !      include 'COMMON.INTERACT'
1388 !      include 'COMMON.FFIELD'
1389 !      include 'COMMON.NAMES'
1390 !      include 'COMMON.LOCAL'
1391       logical :: lprint
1392       integer :: i,j,k,ist,ien,iteli,itelj,ind,i1,i2,it1,it2,ic1,ic2
1393       real(kind=8) :: rri,xi,yi,zi,dxi,dyi,dzi,xmedi,ymedi,zmedi,&
1394         xj,yj,zj,dxj,dyj,dzj,aaa,bbb,ael6i,ael3i,rrmij,rmij,r3ij,r6ij,&
1395         vrmij,cosa,cosb,cosg,fac,ev1,ev2,fac3,fac4,evdwij,el1,el2,&
1396         eesij,ees,evdw,ene
1397       real(kind=8),dimension(2,2) :: elpp6c=reshape((/-0.2379d0,&
1398                        -0.2056d0,-0.2056d0,-0.0610d0/),shape(elpp6c))
1399       real(kind=8),dimension(2,2) :: elpp3c=reshape((/ 0.0503d0,&
1400                         0.0000d0, 0.0000d0, 0.0692d0/),shape(elpp3c))
1401       real(kind=8),dimension(2,2) :: ael6c,ael3c,appc,bppc
1402       real(kind=8) :: elcutoff=-0.3d0
1403       real(kind=8) :: elecutoff_14=-0.5d0
1404       integer :: ncont,icont(2,maxcont),mnum
1405       real(kind=8) :: econt(maxcont)
1406 !
1407 ! Load the constants of peptide bond - peptide bond interactions.
1408 ! Type 1 - ordinary peptide bond, type 2 - alkylated peptide bond (e.g.
1409 ! proline) - determined by averaging ECEPP energy.
1410 !
1411 ! as of 7/06/91.
1412 !
1413 !      data epp    / 0.3045d0, 0.3649d0, 0.3649d0, 0.5743d0/
1414 !      data rpp    / 4.5088d0, 4.5395d0, 4.5395d0, 4.4846d0/
1415 !el      data (elpp6c)   /-0.2379d0,-0.2056d0,-0.2056d0,-0.0610d0/
1416 !el      data (elpp3c)   / 0.0503d0, 0.0000d0, 0.0000d0, 0.0692d0/
1417 !el      data (elcutoff) /-0.3d0/
1418 !el      data (elecutoff_14) /-0.5d0/
1419       ees=0.0d0
1420       evdw=0.0d0
1421       if (lprint) write (iout,'(a)') &
1422         "Constants of electrostatic interaction energy expression."
1423       do i=1,2
1424         do j=1,2
1425         rri=rpp(i,j)**6
1426         appc(i,j)=epp(i,j)*rri*rri
1427         bppc(i,j)=-2.0*epp(i,j)*rri
1428         ael6c(i,j)=elpp6c(i,j)*4.2**6
1429         ael3c(i,j)=elpp3c(i,j)*4.2**3
1430         if (lprint) &
1431         write (iout,'(2i2,4e15.4)') i,j,appc(i,j),bppc(i,j),ael6c(i,j),&
1432                                      ael3c(i,j)
1433         enddo
1434       enddo
1435       ncont=0
1436       do 1 i=ist,ien-2
1437         xi=c(1,i)
1438         yi=c(2,i)
1439         zi=c(3,i)
1440         dxi=c(1,i+1)-c(1,i)
1441         dyi=c(2,i+1)-c(2,i)
1442         dzi=c(3,i+1)-c(3,i)
1443         xmedi=xi+0.5*dxi
1444         ymedi=yi+0.5*dyi
1445         zmedi=zi+0.5*dzi
1446         do 4 j=i+2,ien-1
1447 !          ind=ind+1
1448           iteli=itel(i)
1449           itelj=itel(j)
1450           if (j.eq.i+2 .and. itelj.eq.2) iteli=2
1451           if (iteli.eq.2 .and. itelj.eq.2 &
1452             .or.iteli.eq.0 .or.itelj.eq.0) goto 4
1453           aaa=appc(iteli,itelj)
1454           bbb=bppc(iteli,itelj)
1455           ael6i=ael6c(iteli,itelj)
1456           ael3i=ael3c(iteli,itelj)
1457           dxj=c(1,j+1)-c(1,j)
1458           dyj=c(2,j+1)-c(2,j)
1459           dzj=c(3,j+1)-c(3,j)
1460           xj=c(1,j)+0.5*dxj-xmedi
1461           yj=c(2,j)+0.5*dyj-ymedi
1462           zj=c(3,j)+0.5*dzj-zmedi
1463           rrmij=1.0/(xj*xj+yj*yj+zj*zj)
1464           rmij=sqrt(rrmij)
1465           r3ij=rrmij*rmij
1466           r6ij=r3ij*r3ij
1467           vrmij=vblinv*rmij
1468           cosa=(dxi*dxj+dyi*dyj+dzi*dzj)*vblinv2
1469           cosb=(xj*dxi+yj*dyi+zj*dzi)*vrmij
1470           cosg=(xj*dxj+yj*dyj+zj*dzj)*vrmij
1471           fac=cosa-3.0*cosb*cosg
1472           ev1=aaa*r6ij*r6ij
1473           ev2=bbb*r6ij
1474           fac3=ael6i*r6ij
1475           fac4=ael3i*r3ij
1476           evdwij=ev1+ev2
1477           el1=fac3*(4.0+fac*fac-3.0*(cosb*cosb+cosg*cosg))
1478           el2=fac4*fac
1479           eesij=el1+el2
1480           if (j.gt.i+2 .and. eesij.le.elcutoff .or. &
1481               j.eq.i+2 .and. eesij.le.elecutoff_14) then
1482              ncont=ncont+1
1483              icont(1,ncont)=i
1484              icont(2,ncont)=j
1485              econt(ncont)=eesij
1486           endif
1487           ees=ees+eesij
1488           evdw=evdw+evdwij
1489     4   continue
1490     1 continue
1491       if (lprint) then
1492         write (iout,*) 'Total average electrostatic energy: ',ees
1493         write (iout,*) 'VDW energy between peptide-group centers: ',evdw
1494         write (iout,*)
1495         write (iout,*) 'Electrostatic contacts before pruning: '
1496         do i=1,ncont
1497           i1=icont(1,i)
1498           i2=icont(2,i)
1499           it1=itype(i1,molnum(i1))
1500           it2=itype(i2,molnum(i1))
1501           write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)') &
1502            i,restyp(it1,molnum(i1)),i1,restyp(it2,molnum(i1)),i2,econt(i)
1503         enddo
1504       endif
1505 ! For given residues keep only the contacts with the greatest energy.
1506       i=0
1507       do while (i.lt.ncont)
1508         i=i+1
1509         ene=econt(i)
1510         ic1=icont(1,i)
1511         ic2=icont(2,i)
1512         j=i
1513         do while (j.lt.ncont)
1514           j=j+1
1515           if (ic1.eq.icont(1,j).and.iabs(icont(2,j)-ic2).le.2 .or. &
1516               ic2.eq.icont(2,j).and.iabs(icont(1,j)-ic1).le.2) then
1517 !            write (iout,*) "i",i," j",j," ic1",ic1," ic2",ic2,
1518 !     &       " jc1",icont(1,j)," jc2",icont(2,j)," ncont",ncont
1519             if (econt(j).lt.ene .and. icont(2,j).ne.icont(1,j)+2) then
1520               if (ic1.eq.icont(1,j)) then
1521                 do k=1,ncont
1522                   if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.icont(2,j)&
1523                      .and. iabs(icont(1,k)-ic1).le.2 .and. &
1524                      econt(k).lt.econt(j) ) goto 21
1525                 enddo
1526               else if (ic2.eq.icont(2,j) ) then
1527                 do k=1,ncont
1528                   if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.icont(1,j)&
1529                      .and. iabs(icont(2,k)-ic2).le.2 .and. &
1530                      econt(k).lt.econt(j) ) goto 21
1531                 enddo
1532               endif
1533 ! Remove ith contact
1534               do k=i+1,ncont
1535                 icont(1,k-1)=icont(1,k)
1536                 icont(2,k-1)=icont(2,k)
1537                 econt(k-1)=econt(k)
1538               enddo
1539               i=i-1
1540               ncont=ncont-1
1541 !              write (iout,*) "ncont",ncont
1542 !              do k=1,ncont
1543 !                write (iout,*) icont(1,k),icont(2,k)
1544 !              enddo
1545               goto 20
1546             else if (econt(j).gt.ene .and. ic2.ne.ic1+2) &
1547             then
1548               if (ic1.eq.icont(1,j)) then
1549                 do k=1,ncont
1550                   if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.ic2 &
1551                      .and. iabs(icont(1,k)-icont(1,j)).le.2 .and. &
1552                      econt(k).lt.econt(i) ) goto 21
1553                 enddo
1554               else if (ic2.eq.icont(2,j) ) then
1555                 do k=1,ncont
1556                   if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.ic1 &
1557                      .and. iabs(icont(2,k)-icont(2,j)).le.2 .and. &
1558                      econt(k).lt.econt(i) ) goto 21
1559                 enddo
1560               endif
1561 ! Remove jth contact
1562               do k=j+1,ncont
1563                 icont(1,k-1)=icont(1,k)
1564                 icont(2,k-1)=icont(2,k)
1565                 econt(k-1)=econt(k)
1566               enddo
1567               ncont=ncont-1
1568 !              write (iout,*) "ncont",ncont
1569 !              do k=1,ncont
1570 !                write (iout,*) icont(1,k),icont(2,k)
1571 !              enddo
1572               j=j-1
1573             endif
1574           endif
1575    21     continue
1576         enddo
1577    20   continue
1578       enddo
1579       if (lprint) then
1580         write (iout,*)
1581         write (iout,*) 'Electrostatic contacts after pruning: '
1582         do i=1,ncont
1583 !          mnum=molnum(i)
1584           i1=icont(1,i)
1585           i2=icont(2,i)
1586             mnum=molnum(i1)
1587             it1=itype(i1,mnum)
1588             it2=itype(i2,molnum(i2))
1589           write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)') &
1590            i,restyp(it1,mnum),i1,restyp(it2,molnum(i2)),i2,econt(i)
1591         enddo
1592       endif
1593       return
1594       end subroutine elecont
1595 !------------------------------------------------------------------------------
1596 ! match_contact.f
1597 !------------------------------------------------------------------------------
1598       subroutine match_contact(ishif1,ishif2,nc_match,nc_match1_max,&
1599          ncont_ref,icont_ref,ncont,icont,jfrag,n_shif1,n_shif2,&
1600          nc_frac,nc_req_set,istr,llocal,lprn)
1601
1602       use energy_data, only:maxcont
1603 !      implicit real*8 (a-h,o-z)
1604 !      include 'DIMENSIONS'
1605 !      include 'COMMON.IOUNITS'
1606       integer :: ncont_ref,ncont,ishift,ishif2,nc_match
1607       integer,dimension(2,maxcont) :: icont_ref,icont !(2,maxcont)
1608       real(kind=8) :: nc_frac
1609       logical :: llocal,lprn
1610       integer :: ishif1,nc_match1_max,jfrag,n_shif1,n_shif2,&
1611                  nc_req_set,istr,nc_match_max
1612       integer :: i,nc_req,nc_match1,is,js
1613       nc_match_max=0
1614       do i=1,ncont_ref
1615         nc_match_max=nc_match_max+ &
1616          min0(icont_ref(2,i)-icont_ref(1,i)-1,3)
1617       enddo
1618       if (istr.eq.3) then
1619         nc_req=0
1620       else if (nc_req_set.eq.0) then
1621         nc_req=nc_match_max*nc_frac
1622       else
1623         nc_req = dmin1(nc_match_max*nc_frac+0.5d0,&
1624           dfloat(nc_req_set)+1.0d-7)
1625       endif
1626 !      write (iout,*) "match_contact: nc_req:",nc_req
1627 !      write (iout,*) "nc_match_max",nc_match_max
1628 !      write (iout,*) "jfrag",jfrag," n_shif1",n_shif1,
1629 !     &   " n_shif2",n_shif2
1630 ! Match current contact map against reference contact map; exit, if at least
1631 ! half of the contacts match
1632       call ncont_match(nc_match,nc_match1,0,0,ncont_ref,icont_ref,&
1633           ncont,icont,jfrag,llocal,lprn)
1634       nc_match1_max=nc_match1
1635       if (lprn .and. nc_match.gt.0) write (iout,*) &
1636         "Shift:",0,0," nc_match1",nc_match1,&
1637         " nc_match=",nc_match," req'd",nc_req
1638       if (nc_req.gt.0 .and. nc_match.ge.nc_req .or. &
1639           nc_req.eq.0 .and. nc_match.eq.1) then
1640          ishif1=0
1641          ishif2=0
1642          return
1643       endif
1644 ! If sufficient matches are not found, try to shift contact maps up to three
1645 ! positions.
1646       if (n_shif1.gt.0) then
1647       do is=1,n_shif1
1648 ! The following four tries help to find shifted beta-sheet patterns
1649 ! Shift "left" strand backward
1650         call ncont_match(nc_match,nc_match1,-is,0,ncont_ref,&
1651           icont_ref,ncont,icont,jfrag,llocal,lprn)
1652         if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
1653         if (lprn .and. nc_match.gt.0) write (iout,*) &
1654           "Shift:",-is,0," nc_match1",nc_match1,&
1655           " nc_match=",nc_match," req'd",nc_req
1656         if (nc_req.gt.0 .and. nc_match.ge.nc_req .or. &
1657            nc_req.eq.0 .and. nc_match.eq.1) then
1658           ishif1=-is
1659           ishif2=0
1660           return
1661         endif
1662 ! Shift "left" strand forward
1663         call ncont_match(nc_match,nc_match1,is,0,ncont_ref,&
1664             icont_ref,ncont,icont,jfrag,llocal,lprn)
1665         if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
1666         if (lprn .and. nc_match.gt.0) write (iout,*) &
1667          "Shift:",is,0," nc_match1",nc_match1,&
1668          " nc_match=",nc_match," req'd",nc_req
1669         if (nc_req.gt.0 .and. nc_match.ge.nc_req .or. &
1670            nc_req.eq.0 .and. nc_match.eq.1) then
1671           ishif1=is
1672           ishif2=0
1673           return
1674         endif
1675       enddo
1676       if (nc_req.eq.0) return
1677 ! Shift "right" strand backward
1678       do is=1,n_shif1
1679         call ncont_match(nc_match,nc_match1,0,-is,ncont_ref,&
1680            icont_ref,ncont,icont,jfrag,llocal,lprn)
1681         if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
1682         if (lprn .and. nc_match.gt.0) write (iout,*) &
1683           "Shift:",0,-is," nc_match1",nc_match1,&
1684           " nc_match=",nc_match," req'd",nc_req
1685         if (nc_match.ge.nc_req) then
1686           ishif1=0
1687           ishif2=-is
1688           return
1689         endif
1690 ! Shift "right" strand upward
1691         call ncont_match(nc_match,nc_match1,0,is,ncont_ref,&
1692           icont_ref,ncont,icont,jfrag,llocal,lprn)
1693         if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
1694         if (lprn .and. nc_match.gt.0) write (iout,*) &
1695           "Shift:",0,is," nc_match1",nc_match1,&
1696           " nc_match=",nc_match," req'd",nc_req
1697         if (nc_match.ge.nc_req) then
1698           ishif1=0
1699           ishif2=is
1700           return
1701         endif
1702       enddo ! is
1703 ! Now try to shift both residues in contacts.
1704       do is=1,n_shif1
1705         do js=1,is
1706           if (js.ne.is) then
1707             call ncont_match(nc_match,nc_match1,-is,-js,ncont_ref,&
1708               icont_ref,ncont,icont,jfrag,llocal,lprn)
1709             if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
1710             if (lprn .and. nc_match.gt.0) write (iout,*) &
1711                "Shift:",-is,-js," nc_match1",nc_match1,&
1712                " nc_match=",nc_match," req'd",nc_req
1713             if (nc_match.ge.nc_req) then
1714               ishif1=-is
1715               ishif2=-js
1716               return
1717             endif
1718             call ncont_match(nc_match,nc_match1,is,js,ncont_ref,&
1719               icont_ref,ncont,icont,jfrag,llocal,lprn)
1720             if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
1721             if (lprn .and. nc_match.gt.0) write (iout,*) &
1722               "Shift:",is,js," nc_match1",nc_match1,&
1723               " nc_match=",nc_match," req'd",nc_req
1724             if (nc_match.ge.nc_req) then
1725               ishif1=is
1726               ishif2=js
1727               return
1728             endif
1729 !
1730             call ncont_match(nc_match,nc_match1,-js,-is,ncont_ref,&
1731               icont_ref,ncont,icont,jfrag,llocal,lprn)
1732             if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
1733             if (lprn .and. nc_match.gt.0) write (iout,*) &
1734               "Shift:",-js,-is," nc_match1",nc_match1,&
1735               " nc_match=",nc_match," req'd",nc_req
1736             if (nc_match.ge.nc_req) then
1737               ishif1=-js
1738               ishif2=-is
1739               return
1740             endif
1741 !
1742             call ncont_match(nc_match,nc_match1,js,is,ncont_ref,&
1743               icont_ref,ncont,icont,jfrag,llocal,lprn)
1744             if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
1745             if (lprn .and. nc_match.gt.0) write (iout,*) &
1746               "Shift:",js,is," nc_match1",nc_match1,&
1747               " nc_match=",nc_match," req'd",nc_req
1748             if (nc_match.ge.nc_req) then
1749               ishif1=js
1750               ishif2=is
1751               return
1752             endif
1753           endif
1754 !
1755           if (is+js.le.n_shif1) then
1756           call ncont_match(nc_match,nc_match1,-is,js,ncont_ref,&
1757             icont_ref,ncont,icont,jfrag,llocal,lprn)
1758           if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
1759           if (lprn .and. nc_match.gt.0) write (iout,*) &
1760            "Shift:",-is,js," nc_match1",nc_match1,&
1761            " nc_match=",nc_match," req'd",nc_req
1762           if (nc_match.ge.nc_req) then
1763             ishif1=-is
1764             ishif2=js
1765             return
1766           endif
1767 !
1768           call ncont_match(nc_match,nc_match1,js,-is,ncont_ref,&
1769             icont_ref,ncont,icont,jfrag,llocal,lprn)
1770           if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
1771           if (lprn .and. nc_match.gt.0) write (iout,*) &
1772            "Shift:",js,-is," nc_match1",nc_match1,&
1773            " nc_match=",nc_match," req'd",nc_req
1774           if (nc_match.ge.nc_req) then
1775             ishif1=js
1776             ishif2=-is
1777             return
1778           endif
1779           endif
1780 !
1781         enddo !js
1782       enddo !is
1783       endif
1784
1785       if (n_shif2.gt.0) then
1786       do is=1,n_shif2
1787         call ncont_match(nc_match,nc_match1,-is,-is,ncont_ref,&
1788           icont_ref,ncont,icont,jfrag,llocal,lprn)
1789         if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
1790         if (lprn .and. nc_match.gt.0) write (iout,*) &
1791            "Shift:",-is,-is," nc_match1",nc_match1,&
1792            " nc_match=",nc_match," req'd",nc_req
1793         if (nc_match.ge.nc_req) then
1794           ishif1=-is
1795           ishif2=-is
1796           return
1797         endif
1798         call ncont_match(nc_match,nc_match1,is,is,ncont_ref,&
1799           icont_ref,ncont,icont,jfrag,llocal,lprn)
1800         if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
1801         if (lprn .and. nc_match.gt.0) write (iout,*) &
1802           "Shift:",is,is," nc_match1",nc_match1,&
1803           " nc_match=",nc_match," req'd",nc_req
1804         if (nc_match.ge.nc_req) then
1805           ishif1=is
1806           ishif2=is
1807           return
1808         endif
1809       enddo
1810       endif
1811 ! If this point is reached, the contact maps are different.
1812       nc_match=0
1813       ishif1=0
1814       ishif2=0
1815       return
1816       end subroutine match_contact
1817 !-------------------------------------------------------------------------
1818       subroutine ncont_match(nc_match,nc_match1,ishif1,ishif2,&
1819          ncont_ref,icont_ref,ncont,icont,jfrag,llocal,lprn)
1820
1821       use energy_data, only:nnt,nct,maxcont
1822 !      implicit real*8 (a-h,o-z)
1823 !      include 'DIMENSIONS'
1824 !      include 'DIMENSIONS.ZSCOPT'
1825 !      include 'DIMENSIONS.COMPAR'
1826 !      include 'COMMON.IOUNITS'
1827 !      include 'COMMON.INTERACT'
1828 !      include 'COMMON.GEO'
1829 !      include 'COMMON.COMPAR'
1830       logical :: llocal,lprn
1831       integer ncont_ref,ncont,ishift,ishif2,nang_pair
1832       integer,dimension(2,maxcont) :: icont_ref,icont,icont_match !(2,maxcont)
1833       integer,dimension(2,nres) :: iang_pair !(2,maxres)
1834       integer :: nc_match,nc_match1,ishif1,jfrag
1835       integer :: i,j,ic1,ic2
1836       real(kind=8) :: diffang,fract,rad2deg
1837
1838 ! Compare the contact map against the reference contact map; they're stored
1839 ! in ICONT and ICONT_REF, respectively. The current contact map can be shifted.
1840       if (lprn) write (iout,'(80(1h*))')
1841       nc_match=0
1842       nc_match1=0
1843 ! Check the local structure by comparing dihedral angles.
1844 !      write (iout,*) "ncont_match: ncont_ref",ncont_ref," llocal",llocal
1845       if (llocal .and. ncont_ref.eq.0) then
1846 ! If there are no contacts just compare the dihedral angles and exit.
1847         call angnorm(jfrag,ishif1,ishif2,ang_cut1(jfrag),diffang,fract,&
1848           lprn)
1849         if (lprn) write (iout,*) "diffang:",diffang*rad2deg,&
1850          " ang_cut:",ang_cut(jfrag)*rad2deg," fract",fract
1851         if (diffang.le.ang_cut(jfrag) .and. fract.ge.frac_min(jfrag)) &
1852         then
1853           nc_match=1
1854         else
1855           nc_match=0
1856         endif
1857         return
1858       endif
1859       nang_pair=0
1860       do i=1,ncont
1861         ic1=icont(1,i)+ishif1
1862         ic2=icont(2,i)+ishif2
1863 !        write (iout,*) "i",i," ic1",ic1," ic2",ic2
1864         if (ic1.lt.nnt .or. ic2.gt.nct) goto 10
1865         do j=1,ncont_ref
1866           if (ic1.eq.icont_ref(1,j).and.ic2.eq.icont_ref(2,j)) then
1867             nc_match=nc_match+min0(icont_ref(2,j)-icont_ref(1,j)-1,3)
1868             nc_match1=nc_match1+1
1869             icont_match(1,nc_match1)=ic1
1870             icont_match(2,nc_match1)=ic2
1871 !            call add_angpair(icont(1,i),icont_ref(1,j),
1872 !     &         nang_pair,iang_pair)
1873 !            call add_angpair(icont(2,i),icont_ref(2,j),
1874 !     &         nang_pair,iang_pair)
1875             if (lprn) write (iout,*) "Contacts:",icont(1,i),icont(2,i),&
1876              " match",icont_ref(1,j),icont_ref(2,j),&
1877              " shifts",ishif1,ishif2
1878             goto 10
1879           endif
1880         enddo
1881    10   continue
1882       enddo
1883       if (lprn) then
1884         write (iout,*) "nc_match",nc_match," nc_match1",nc_match1
1885         write (iout,*) "icont_match"
1886         do i=1,nc_match1
1887           write (iout,*) icont_match(1,i),icont_match(2,i)
1888         enddo
1889       endif
1890       if (llocal .and. nc_match.gt.0) then
1891         call angnorm2(jfrag,ishif1,ishif2,nc_match1,icont_match,lprn,&
1892           ang_cut1(jfrag),diffang,fract)
1893         if (lprn) write (iout,*) "diffang:",diffang*rad2deg,&
1894          " ang_cut:",ang_cut(jfrag)*rad2deg,&
1895          " ang_cut1",ang_cut1(jfrag)*rad2deg
1896         if (diffang.gt.ang_cut(jfrag) &
1897           .or. fract.lt.frac_min(jfrag)) nc_match=0
1898       endif
1899 !      if (nc_match.gt.0) then
1900 !        diffang = angnorm1(nang_pair,iang_pair,lprn)
1901 !        if (diffang.gt.ang_cut(jfrag)) nc_match=0
1902 !      endif
1903       if (lprn) write (iout,*) "ishif1",ishif1," ishif2",ishif2,&
1904          " diffang",rad2deg*diffang," nc_match",nc_match
1905       return
1906       end subroutine ncont_match
1907 !------------------------------------------------------------------------------
1908       subroutine match_secondary(jfrag,isecstr,nsec_match,lprn)
1909 ! This subroutine compares the secondary structure (isecstr) of fragment jfrag
1910 ! conformation considered to that of the reference conformation.
1911 ! Returns the number of equivalent residues (nsec_match).
1912 !      implicit real*8 (a-h,o-z)
1913 !      include 'DIMENSIONS'
1914 !      include 'DIMENSIONS.ZSCOPT'
1915 !      include 'DIMENSIONS.COMPAR'
1916 !      include 'COMMON.IOUNITS'
1917 !      include 'COMMON.CHAIN'
1918 !      include 'COMMON.PEPTCONT'
1919 !      include 'COMMON.COMPAR'
1920       logical :: lprn
1921       integer :: isecstr(nres)
1922       integer :: jfrag,nsec_match,npart,i,j
1923       npart = npiece(jfrag,1)
1924       nsec_match=0
1925       if (lprn) then
1926         write (iout,*) "match_secondary jfrag",jfrag," ifrag",&
1927               (ifrag(1,i,jfrag),ifrag(2,i,jfrag),i=1,npart)
1928         write (iout,'(80i1)') (isec_ref(j),j=1,nres)
1929         write (iout,'(80i1)') (isecstr(j),j=1,nres)
1930       endif
1931       do i=1,npart
1932         do j=ifrag(1,i,jfrag),ifrag(2,i,jfrag)
1933 ! The residue has equivalent conformational state to that of the reference
1934 ! structure, if:
1935 !  a) the conformational states are equal or
1936 !  b) the reference state is a coil and that of the conformation considered
1937 !     is a strand or
1938 !  c) the conformational state of the conformation considered is a strand
1939 !     and that of the reference conformation is a coil.
1940 ! 10/28/02 - case (b) deleted.
1941           if (isecstr(j).eq.isec_ref(j) .or. &
1942 !     &        isecstr(j).eq.0 .and. isec_ref(j).eq.1 .or.
1943               isec_ref(j).eq.0 .and. isecstr(j).eq.1) &
1944             nsec_match=nsec_match+1
1945         enddo
1946       enddo
1947       return
1948       end subroutine match_secondary
1949 !------------------------------------------------------------------------------
1950 ! odlodc.f
1951 !------------------------------------------------------------------------------
1952       subroutine odlodc(r1,r2,a,b,uu,vv,aa,bb,dd)
1953
1954       use energy_data, only:ncont_ref,icont_ref!,nsccont_frag_ref,&
1955 !                            isccont_frag_ref
1956 !      implicit real*8 (a-h,o-z)
1957       real(kind=8),dimension(3) :: r1,r2,a,b,x,y
1958       real(kind=8) :: uu,vv,aa,bb,dd
1959       real(kind=8) :: ab,ar,br,det,dd1,dd2,dd3,dd4,dd5
1960 !el      odl(u,v) = (r1(1)-r2(1))**2+(r1(2)-r2(2))**2+(r1(3)-r2(3))**2 &
1961 !el       + 2*ar*u - 2*br*v - 2*ab*u*v + aa*u**2 + bb*v**2
1962 !      print *,"r1",(r1(i),i=1,3)
1963 !      print *,"r2",(r2(i),i=1,3)
1964 !      print *,"a",(a(i),i=1,3)
1965 !      print *,"b",(b(i),i=1,3)
1966       aa = a(1)**2+a(2)**2+a(3)**2
1967       bb = b(1)**2+b(2)**2+b(3)**2
1968       ab = a(1)*b(1)+a(2)*b(2)+a(3)*b(3)
1969       ar = a(1)*(r1(1)-r2(1))+a(2)*(r1(2)-r2(2))+a(3)*(r1(3)-r2(3))
1970       br = b(1)*(r1(1)-r2(1))+b(2)*(r1(2)-r2(2))+b(3)*(r1(3)-r2(3))
1971       det = aa*bb-ab**2
1972 !      print *,'aa',aa,' bb',bb,' ab',ab,' ar',ar,' br',br,' det',det
1973       uu = (-ar*bb+br*ab)/det
1974       vv = (br*aa-ar*ab)/det
1975 !      print *,u,v
1976       uu=dmin1(uu,1.0d0)
1977       uu=dmax1(uu,0.0d0)
1978       vv=dmin1(vv,1.0d0)
1979       vv=dmax1(vv,0.0d0)
1980 !el      dd1 = odl(uu,vv)
1981       dd1 = odl(uu,vv,r1,r2,ar,br,ab,aa,bb)
1982 !el      dd2 = odl(0.0d0,0.0d0)
1983       dd2 = odl(0.0d0,0.0d0,r1,r2,ar,br,ab,aa,bb)
1984 !el      dd3 = odl(0.0d0,1.0d0)
1985       dd3 = odl(0.0d0,1.0d0,r1,r2,ar,br,ab,aa,bb)
1986 !el      dd4 = odl(1.0d0,0.0d0)
1987       dd4 = odl(1.0d0,0.0d0,r1,r2,ar,br,ab,aa,bb)
1988 !el      dd5 = odl(1.0d0,1.0d0)
1989       dd5 = odl(1.0d0,1.0d0,r1,r2,ar,br,ab,aa,bb)
1990       dd = dsqrt(dmin1(dd1,dd2,dd3,dd4,dd5))
1991       if (dd.eq.dd2) then
1992         uu=0.0d0
1993         vv=0.0d0
1994       else if (dd.eq.dd3) then
1995         uu=0.0d0
1996         vv=1.0d0
1997       else if (dd.eq.dd4) then
1998         uu=1.0d0
1999         vv=0.0d0
2000       else if (dd.eq.dd5) then
2001         uu=1.0d0
2002         vv=1.0d0
2003       endif
2004 ! Control check
2005 !      do i=1,3
2006 !        x(i)=r1(i)+u*a(i)
2007 !        y(i)=r2(i)+v*b(i)
2008 !      enddo
2009 !      dd1 = (x(1)-y(1))**2+(x(2)-y(2))**2+(x(3)-y(3))**2
2010 !      dd1 = dsqrt(dd1)
2011       aa = dsqrt(aa)
2012       bb = dsqrt(bb)
2013 !      write (8,*) uu,vv,dd,dd1
2014 !      print *,dd,dd1
2015       return
2016       end subroutine odlodc
2017 !------------------------------------------------------------------------------
2018       real(kind=8) function odl(u,v,r1,r2,ar,br,ab,aa,bb)
2019
2020       real(kind=8),dimension(3) :: r1,r2
2021       real(kind=8) :: aa,bb,u,v,ar,br,ab
2022
2023       odl = (r1(1)-r2(1))**2+(r1(2)-r2(2))**2+(r1(3)-r2(3))**2 &
2024        + 2*ar*u - 2*br*v - 2*ab*u*v + aa*u**2 + bb*v**2
2025
2026       end function odl
2027 !------------------------------------------------------------------------------
2028 ! proc_cont.f
2029 !------------------------------------------------------------------------------
2030       subroutine proc_cont
2031
2032       use geometry_data, only:rad2deg
2033       use energy_data, only:ncont_ref,icont_ref!,nsccont_frag_ref,&
2034 !                            isccont_frag_ref
2035 !      implicit real*8 (a-h,o-z)
2036 !      include 'DIMENSIONS'
2037 !      include 'DIMENSIONS.ZSCOPT'
2038 !      include 'DIMENSIONS.COMPAR'
2039 !      include 'COMMON.IOUNITS'
2040 !      include 'COMMON.TIME1'
2041 !      include 'COMMON.SBRIDGE'
2042 !      include 'COMMON.CONTROL'
2043 !      include 'COMMON.COMPAR'
2044 !      include 'COMMON.CHAIN'
2045 !      include 'COMMON.HEADER'
2046 !      include 'COMMON.CONTACTS1'
2047 !      include 'COMMON.PEPTCONT'
2048 !      include 'COMMON.GEO'
2049       integer :: i,j,k,ind,len_cut,ndigit,length_frag
2050
2051       write (iout,*) "proc_cont: nlevel",nlevel
2052       if (nlevel.lt.0) then
2053         write (iout,*) "call define_fragments"
2054         call define_fragments
2055       else
2056         write (iout,*) "call secondary2"
2057         call secondary2(.true.,.false.,ncont_pept_ref,icont_pept_ref,&
2058            isec_ref)
2059       endif
2060       write (iout,'(80(1h=))')
2061       write (iout,*) "Electrostatic contacts"
2062       call contacts_between_fragments(.true.,0,ncont_pept_ref,&
2063        icont_pept_ref,ncont_frag_ref(1),icont_frag_ref(1,1,1))
2064       write (iout,'(80(1h=))')
2065       write (iout,*) "Side chain contacts"
2066       call contacts_between_fragments(.true.,0,ncont_ref,&
2067        icont_ref,nsccont_frag_ref(1),isccont_frag_ref(1,1,1))
2068       if (nlevel.lt.0) then
2069         do i=1,nfrag(1)
2070           ind=icant(i,i)
2071           len_cut=1000
2072           if (istruct(i).le.1) then
2073             len_cut=max0(len_frag(i,1)*4/5,3)
2074           else if (istruct(i).eq.2 .or. istruct(i).eq.4) then
2075             len_cut=max0(len_frag(i,1)*2/5,3)
2076           endif
2077           write (iout,*) "i",i," istruct",istruct(i)," ncont_frag",&
2078             ncont_frag_ref(ind)," len_cut",len_cut,&
2079             " icont_single",icont_single," iloc_single",iloc_single
2080           iloc(i)=iloc_single
2081           if (iloc(i).gt.0) write (iout,*) &
2082            "Local structure used to compare structure of fragment",i,&
2083            " to native."
2084           if (istruct(i).ne.3 .and. istruct(i).ne.0 &
2085               .and. icont_single.gt.0 .and. &
2086               ncont_frag_ref(ind).ge.len_cut) then
2087             write (iout,*) "Electrostatic contacts used to compare",&
2088              " structure of fragment",i," to native."
2089             ielecont(i,1)=1
2090             isccont(i,1)=0
2091           else if (icont_single.gt.0 .and. nsccont_frag_ref(ind) &
2092             .ge.len_cut) then
2093             write (iout,*) "Side chain contacts used to compare",&
2094              " structure of fragment",i," to native."
2095             isccont(i,1)=1
2096             ielecont(i,1)=0
2097           else
2098             write (iout,*) "Contacts not used to compare",&
2099              " structure of fragment",i," to native."
2100             ielecont(i,1)=0
2101             isccont(i,1)=0
2102             nc_req_setf(i,1)=0
2103           endif
2104           if (irms_single.gt.0 .or. isccont(i,1).eq.0 &
2105                .and. ielecont(i,1).eq.0) then
2106             write (iout,*) "RMSD used to compare",&
2107              " structure of fragment",i," to native."
2108             irms(i,1)=1
2109           else
2110             write (iout,*) "RMSD not used to compare",&
2111              " structure of fragment",i," to native."
2112             irms(i,1)=0
2113           endif
2114         enddo
2115       endif
2116       if (nlevel.lt.-1) then
2117         call define_pairs
2118         nlevel = -nlevel
2119         if (nlevel.gt.3) nlevel=3
2120         if (nlevel.eq.3) then
2121           nfrag(3)=1
2122           npiece(1,3)=nfrag(1)
2123           do i=1,nfrag(1)
2124             ipiece(i,1,3)=i
2125           enddo
2126           ielecont(1,3)=0
2127           isccont(1,3)=0
2128           irms(1,3)=1
2129           n_shift(1,1,3)=0
2130           n_shift(2,1,3)=0
2131         endif
2132       else if (nlevel.eq.-1) then
2133         nlevel=1
2134       endif
2135       isnfrag(1)=0
2136       do i=1,nlevel
2137         isnfrag(i+1)=isnfrag(i)+nfrag(i)
2138       enddo
2139       ndigit=3*nfrag(1)
2140       do i=2,nlevel
2141         ndigit=ndigit+2*nfrag(i)
2142       enddo
2143       write (iout,*) "ndigit",ndigit
2144       if (.not.binary .and. ndigit.gt.30) then
2145         write (iout,*) "Highest class too large; switching to",&
2146           " binary representation."
2147         binary=.true.
2148       endif
2149       write (iout,*) "isnfrag",(isnfrag(i),i=1,nlevel+1)
2150       write(iout,*) "rmscut_base_up",rmscut_base_up,&
2151        " rmscut_base_low",rmscut_base_low," rmsup_lim",rmsup_lim
2152       do i=1,nlevel
2153         do j=1,nfrag(i)
2154           length_frag = 0
2155           if (i.eq.1) then
2156             do k=1,npiece(j,i)
2157               length_frag=length_frag+ifrag(2,k,j)-ifrag(1,k,j)+1
2158             enddo
2159           else
2160             do k=1,npiece(j,i)
2161               length_frag=length_frag+len_frag(ipiece(k,j,i),1)
2162             enddo
2163           endif
2164           len_frag(j,i)=length_frag
2165           rmscutfrag(1,j,i)=rmscut_base_up*length_frag
2166           rmscutfrag(2,j,i)=rmscut_base_low*length_frag
2167           if (rmscutfrag(1,j,i).lt.rmsup_lim) &
2168             rmscutfrag(1,j,i)=rmsup_lim
2169           if (rmscutfrag(1,j,i).gt.rmsupup_lim) &
2170             rmscutfrag(1,j,i)=rmsupup_lim
2171         enddo
2172       enddo
2173       write (iout,*) "Level",1," number of fragments:",nfrag(1)
2174       do j=1,nfrag(1)
2175         write (iout,*) npiece(j,1),(ifrag(1,k,j),ifrag(2,k,j),&
2176           k=1,npiece(j,1)),len_frag(j,1),rmscutfrag(1,j,1),&
2177           rmscutfrag(2,j,1),n_shift(1,j,1),n_shift(2,j,1),&
2178           ang_cut(j)*rad2deg,ang_cut1(j)*rad2deg,frac_min(j),&
2179           nc_fragm(j,1),nc_req_setf(j,1),istruct(j)
2180       enddo
2181       do i=2,nlevel
2182         write (iout,*) "Level",i," number of fragments:",nfrag(i)
2183         do j=1,nfrag(i)
2184           write (iout,*) npiece(j,i),(ipiece(k,j,i),&
2185             k=1,npiece(j,i)),len_frag(j,i),rmscutfrag(1,j,i),&
2186             rmscutfrag(2,j,i),n_shift(1,j,i),n_shift(2,j,i),&
2187             nc_fragm(j,i),nc_req_setf(j,i)
2188         enddo
2189       enddo
2190       return
2191       end subroutine proc_cont
2192 !------------------------------------------------------------------------------
2193 ! define_pairs.f
2194 !------------------------------------------------------------------------------
2195       subroutine define_pairs
2196
2197 !      use energy_data, only:nsccont_frag_ref
2198 !      implicit real*8 (a-h,o-z)
2199 !      include 'DIMENSIONS'
2200 !      include 'DIMENSIONS.ZSCOPT'
2201 !      include 'DIMENSIONS.COMPAR'
2202 !      include 'COMMON.IOUNITS'
2203 !      include 'COMMON.TIME1'
2204 !      include 'COMMON.SBRIDGE'
2205 !      include 'COMMON.CONTROL'
2206 !      include 'COMMON.COMPAR'
2207 !      include 'COMMON.FRAG'
2208 !      include 'COMMON.CHAIN'
2209 !      include 'COMMON.HEADER'
2210 !      include 'COMMON.GEO'
2211 !      include 'COMMON.CONTACTS1'
2212 !      include 'COMMON.PEPTCONT'
2213       integer :: j,k,i,length_frag,ind,ll1,ll2,len_cut
2214
2215       do j=1,nfrag(1)
2216         length_frag = 0
2217         do k=1,npiece(j,1)
2218           length_frag=length_frag+ifrag(2,k,j)-ifrag(1,k,j)+1
2219         enddo
2220         len_frag(j,1)=length_frag
2221         write (iout,*) "Fragment",j," length",len_frag(j,1)
2222       enddo
2223       nfrag(2)=0
2224       do i=1,nfrag(1)
2225         do j=i+1,nfrag(1)
2226           ind = icant(i,j)
2227           if (istruct(i).le.1 .or. istruct(j).le.1) then
2228             if (istruct(i).le.1) then
2229               ll1=len_frag(i,1)
2230             else
2231               ll1=len_frag(i,1)/2
2232             endif
2233             if (istruct(j).le.1) then
2234               ll2=len_frag(j,1)
2235             else
2236               ll2=len_frag(j,1)/2
2237             endif
2238             len_cut=max0(min0(ll1*2/3,ll2*4/5),3)
2239           else
2240             if (istruct(i).eq.2 .or. istruct(i).eq.4) then
2241               ll1=len_frag(i,1)/2
2242             else
2243               ll1=len_frag(i,1)
2244             endif
2245             if (istruct(j).eq.2 .or. istruct(j).eq.4) then
2246               ll2=len_frag(j,1)/2
2247             else
2248               ll2=len_frag(j,1)
2249             endif
2250             len_cut=max0(min0(ll1*4/5,ll2)*4/5,3)
2251           endif
2252           write (iout,*) "Fragments",i,j," structure",istruct(i),&
2253              istruct(j)," # contacts",&
2254              ncont_frag_ref(ind),nsccont_frag_ref(ind),&
2255              " lengths",len_frag(i,1),len_frag(j,1),&
2256              " ll1",ll1," ll2",ll2," len_cut",len_cut
2257           if ((istruct(i).eq.1 .or. istruct(j).eq.1) .and. &
2258             nsccont_frag_ref(ind).ge.len_cut ) then
2259             if (istruct(i).eq.1 .and. istruct(j).eq.1) then
2260               write (iout,*) "Adding pair of helices",i,j,&
2261               " based on SC contacts"
2262             else
2263               write (iout,*) "Adding helix+strand/sheet pair",i,j,&
2264               " based on SC contacts"
2265             endif
2266             nfrag(2)=nfrag(2)+1
2267             if (icont_pair.gt.0) then
2268               write (iout,*)  "# SC contacts will be used",&
2269               " in comparison."
2270               isccont(nfrag(2),2)=1
2271             endif
2272             if (irms_pair.gt.0) then
2273               write (iout,*)  "Fragment RMSD will be used",&
2274               " in comparison."
2275               irms(nfrag(2),2)=1
2276             endif
2277             npiece(nfrag(2),2)=2
2278             ipiece(1,nfrag(2),2)=i
2279             ipiece(2,nfrag(2),2)=j
2280             ielecont(nfrag(2),2)=0
2281             n_shift(1,nfrag(2),2)=nshift_pair
2282             n_shift(2,nfrag(2),2)=nshift_pair
2283             nc_fragm(nfrag(2),2)=ncfrac_pair
2284             nc_req_setf(nfrag(2),2)=ncreq_pair
2285           else if ((istruct(i).ge.2 .and. istruct(i).le.4) &
2286              .and. (istruct(j).ge.2 .and. istruct(i).le.4) &
2287              .and. ncont_frag_ref(ind).ge.len_cut ) then
2288             nfrag(2)=nfrag(2)+1
2289             write (iout,*) "Adding pair strands/sheets",i,j,&
2290               " based on pp contacts"
2291             if (icont_pair.gt.0) then
2292               write (iout,*) "# pp contacts will be used",&
2293               " in comparison."
2294               ielecont(nfrag(2),2)=1
2295             endif
2296             if (irms_pair.gt.0) then
2297               write (iout,*)  "Fragment RMSD will be used",&
2298               " in comparison."
2299               irms(nfrag(2),2)=1
2300             endif
2301             npiece(nfrag(2),2)=2
2302             ipiece(1,nfrag(2),2)=i
2303             ipiece(2,nfrag(2),2)=j
2304             ielecont(nfrag(2),2)=1
2305             isccont(nfrag(2),2)=0
2306             n_shift(1,nfrag(2),2)=nshift_pair
2307             n_shift(2,nfrag(2),2)=nshift_pair
2308             nc_fragm(nfrag(2),2)=ncfrac_bet
2309             nc_req_setf(nfrag(2),2)=ncreq_bet
2310           endif
2311         enddo
2312       enddo
2313       write (iout,*) "Pairs found"
2314       do i=1,nfrag(2)
2315         write (iout,*) ipiece(1,i,2),ipiece(2,i,2)
2316       enddo
2317       return
2318       end subroutine define_pairs
2319 !------------------------------------------------------------------------------
2320 ! icant.f
2321 !------------------------------------------------------------------------------
2322       INTEGER FUNCTION ICANT(I,J)
2323       integer :: i,j
2324       IF (I.GE.J) THEN
2325         ICANT=(I*(I-1))/2+J
2326       ELSE
2327         ICANT=(J*(J-1))/2+I
2328       ENDIF
2329       RETURN
2330       END FUNCTION ICANT
2331 !------------------------------------------------------------------------------
2332 ! mysort.f
2333 !------------------------------------------------------------------------------
2334       subroutine imysort(n, m, mm, x, y, z, z1, z2, z3, z4, z5, z6)
2335 !      implicit none
2336       integer :: n,m,mm
2337       integer :: x(m,mm,n),y(n),z(n),z1(2,n),z6(n),xmin,xtemp
2338       real(kind=8) :: z2(n),z3(n),z4(n),z5(n)
2339       real(kind=8) :: xxtemp
2340       integer :: i,j,k,imax
2341       do i=1,n
2342         xmin=x(1,1,i)
2343         imax=i
2344         do j=i+1,n
2345           if (x(1,1,j).lt.xmin) then
2346             imax=j
2347             xmin=x(1,1,j)
2348           endif
2349         enddo
2350         xxtemp=z2(imax)
2351         z2(imax)=z2(i)
2352         z2(i)=xxtemp
2353         xxtemp=z3(imax)
2354         z3(imax)=z3(i)
2355         z3(i)=xxtemp
2356         xxtemp=z4(imax)
2357         z4(imax)=z4(i)
2358         z4(i)=xxtemp
2359         xxtemp=z5(imax)
2360         z5(imax)=z5(i)
2361         z5(i)=xxtemp
2362         xtemp=y(imax)
2363         y(imax)=y(i)
2364         y(i)=xtemp
2365         xtemp=z(imax)
2366         z(imax)=z(i)
2367         z(i)=xtemp
2368         xtemp=z6(imax)
2369         z6(imax)=z6(i)
2370         z6(i)=xtemp
2371         do j=1,2
2372           xtemp=z1(j,imax)
2373           z1(j,imax)=z1(j,i)
2374           z1(j,i)=xtemp
2375         enddo
2376         do j=1,m
2377           do k=1,mm
2378             xtemp=x(j,k,imax)
2379             x(j,k,imax)=x(j,k,i)
2380             x(j,k,i)=xtemp
2381           enddo
2382         enddo
2383       enddo
2384       return
2385       end subroutine imysort
2386 !------------------------------------------------------------------------------
2387 ! qwolynes.f
2388 !-------------------------------------------------------------------------------
2389       real(kind=8) function qwolynes(ilevel,jfrag)
2390
2391       use geometry_data, only:cref,nperm
2392       use control_data, only:symetr
2393       use energy_data, only:nnt,nct,itype,molnum
2394 !      implicit none
2395 !      include 'DIMENSIONS'
2396 !      include 'DIMENSIONS.ZSCOPT'
2397 !      include 'DIMENSIONS.COMPAR'
2398 !      include 'COMMON.IOUNITS'
2399 !      include 'COMMON.COMPAR'
2400 !      include 'COMMON.CHAIN'
2401 !      include 'COMMON.INTERACT'
2402 !      include 'COMMON.CONTROL'
2403       integer :: ilevel,jfrag,kkk
2404       integer :: i,j,jl,k,l,il,kl,nl,np,ip,kp
2405       integer :: nsep=3
2406       real(kind=8),dimension(:),allocatable :: tempus !(maxperm)
2407       real(kind=8) :: maxiQ !dist,
2408       real(kind=8) :: qq,qqij,qqijCM,dij,d0ij,dijCM,d0ijCM
2409       logical :: lprn=.false.
2410       real(kind=8) :: x !el sigm
2411 !el      sigm(x)=0.25d0*x
2412       nperm=1
2413       maxiQ=0
2414       do i=1,symetr
2415       nperm=i*nperm
2416       enddo
2417 !      write (iout,*) "QWolyes: " jfrag",jfrag,
2418 !     &  " ilevel",ilevel
2419       allocate(tempus(nperm))
2420       do kkk=1,nperm
2421       qq = 0.0d0
2422       if (ilevel.eq.0) then
2423         if (lprn) write (iout,*) "Q computed for whole molecule"
2424         nl=0
2425         do il=nnt+nsep,nct
2426           do jl=nnt,il-nsep
2427             dij=0.0d0
2428             dijCM=0.0d0
2429             d0ij=0.0d0
2430             d0ijCM=0.0d0
2431             qqij=0.0d0
2432             qqijCM=0.0d0
2433             nl=nl+1
2434             d0ij=dsqrt((cref(1,jl,kkk)-cref(1,il,kkk))**2+ &
2435                        (cref(2,jl,kkk)-cref(2,il,kkk))**2+ &
2436                        (cref(3,jl,kkk)-cref(3,il,kkk))**2)
2437             dij=dist(il,jl)
2438             qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2)
2439             if (itype(il,molnum(il)).ne.10 .or. itype(jl,molnum(jl)).ne.10) then
2440               nl=nl+1
2441               d0ijCM=dsqrt( &
2442                      (cref(1,jl+nres,kkk)-cref(1,il+nres,kkk))**2+ &
2443                      (cref(2,jl+nres,kkk)-cref(2,il+nres,kkk))**2+ &
2444                      (cref(3,jl+nres,kkk)-cref(3,il+nres,kkk))**2)
2445               dijCM=dist(il+nres,jl+nres)
2446               qqijCM = dexp(-0.5d0*((dijCM-d0ijCM)/(sigm(d0ijCM)))**2)
2447             endif
2448             qq = qq+qqij+qqijCM
2449             if (lprn) then
2450               write (iout,*) "il",il," jl",jl,&
2451                " itype",itype(il,molnum(il)),itype(jl,molnum(jl))
2452               write (iout,*)"d0ij",d0ij," dij",dij," d0ijCM",d0ijCM,&
2453                " dijCM",dijCM," qqij",qqij," qqijCM",qqijCM
2454             endif
2455           enddo
2456         enddo
2457         qq = qq/nl
2458         if (lprn) write (iout,*) "nl",nl," qq",qq
2459       else if (ilevel.eq.1) then
2460         if (lprn) write (iout,*) "Level",ilevel," fragment",jfrag
2461         nl=0
2462 !        write (iout,*) "nlist_frag",nlist_frag(jfrag)
2463         do i=2,nlist_frag(jfrag)
2464           do j=1,i-1
2465             il=list_frag(i,jfrag)
2466             jl=list_frag(j,jfrag)
2467             if (iabs(il-jl).gt.nsep) then
2468               dij=0.0d0
2469               dijCM=0.0d0
2470               d0ij=0.0d0
2471               d0ijCM=0.0d0
2472               qqij=0.0d0
2473               qqijCM=0.0d0
2474               nl=nl+1
2475               d0ij=dsqrt((cref(1,jl,kkk)-cref(1,il,kkk))**2+ &
2476                          (cref(2,jl,kkk)-cref(2,il,kkk))**2+ &
2477                          (cref(3,jl,kkk)-cref(3,il,kkk))**2)
2478               dij=dist(il,jl)
2479               qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2)
2480               if (itype(il,molnum(il)).ne.10 .or. itype(jl,molnum(jl)).ne.10) then
2481                 nl=nl+1
2482                 d0ijCM=dsqrt( &
2483                        (cref(1,jl+nres,kkk)-cref(1,il+nres,kkk))**2+ &
2484                        (cref(2,jl+nres,kkk)-cref(2,il+nres,kkk))**2+ &
2485                        (cref(3,jl+nres,kkk)-cref(3,il+nres,kkk))**2)
2486                 dijCM=dist(il+nres,jl+nres)
2487                qqijCM = dexp(-0.5d0*((dijCM-d0ijCM)/(sigm(d0ijCM)))**2)
2488               endif
2489               qq = qq+qqij+qqijCM
2490               if (lprn) then
2491                 write (iout,*) "i",i," j",j," il",il," jl",jl,&
2492                  " itype",itype(il,molnum(il)),itype(jl,molnum(jl))
2493                 write (iout,*)"d0ij",d0ij," dij",dij," d0ijCM",d0ijCM,&
2494                  " dijCM",dijCM," qqij",qqij," qqijCM",qqijCM
2495               endif
2496             endif
2497           enddo
2498         enddo
2499         qq = qq/nl
2500         if (lprn) write (iout,*) "nl",nl," qq",qq
2501       else if (ilevel.eq.2) then
2502         np=npiece(jfrag,ilevel)
2503         nl=0
2504         do i=2,np
2505           ip=ipiece(i,jfrag,ilevel)
2506           do j=1,nlist_frag(ip)
2507             il=list_frag(j,ip)
2508             do k=1,i-1
2509               kp=ipiece(k,jfrag,ilevel)
2510               do l=1,nlist_frag(kp)
2511                 kl=list_frag(l,kp)
2512                 if (iabs(kl-il).gt.nsep) then
2513                   nl=nl+1
2514                   dij=0.0d0
2515                   dijCM=0.0d0
2516                   d0ij=0.0d0
2517                   d0ijCM=0.0d0
2518                   qqij=0.0d0
2519                   qqijCM=0.0d0
2520                   d0ij=dsqrt((cref(1,kl,kkk)-cref(1,il,kkk))**2+ &
2521                              (cref(2,kl,kkk)-cref(2,il,kkk))**2+ &
2522                              (cref(3,kl,kkk)-cref(3,il,kkk))**2)
2523                   dij=dist(il,kl)
2524                   qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2)
2525                   if (itype(il,molnum(il)).ne.10 .or. itype(kl,molnum(kl)).ne.10) then
2526                     nl=nl+1
2527                     d0ijCM=dsqrt( &
2528                        (cref(1,kl+nres,kkk)-cref(1,il+nres,kkk))**2+ &
2529                        (cref(2,kl+nres,kkk)-cref(2,il+nres,kkk))**2+ &
2530                        (cref(3,kl+nres,kkk)-cref(3,il+nres,kkk))**2)
2531                     dijCM=dist(il+nres,kl+nres)
2532                     qqijCM = dexp(-0.5d0*((dijCM-d0ijCM)/ &
2533                       (sigm(d0ijCM)))**2)
2534                   endif
2535                   qq = qq+qqij+qqijCM
2536                   if (lprn) then
2537                     write (iout,*) "i",i," j",j," k",k," l",l," il",il,&
2538                       " kl",kl," itype",itype(il,molnum(il)), &
2539                          itype(kl,molnum(kl))
2540                     write (iout,*) " d0ij",d0ij," dij",dij," d0ijCM",&
2541                     d0ijCM," dijCM",dijCM," qqij",qqij," qqijCM",qqijCM
2542                   endif
2543                 endif
2544               enddo  ! l
2545             enddo    ! k
2546           enddo      ! j
2547         enddo        ! i
2548         qq = qq/nl
2549         if (lprn) write (iout,*) "nl",nl," qq",qq
2550       else
2551         write (iout,*)"Error: Q can be computed only for level 1 and 2."
2552       endif
2553       tempus(kkk)=qq
2554       enddo
2555       do kkk=1,nperm
2556        if (maxiQ.le.tempus(kkk)) maxiQ=tempus(kkk)
2557       enddo
2558       qwolynes=1.0d0-maxiQ
2559       deallocate(tempus)
2560       return
2561       end function qwolynes
2562 !-------------------------------------------------------------------------------
2563       real(kind=8) function sigm(x)
2564       real(kind=8) :: x
2565       sigm=0.25d0*x
2566       return
2567       end function sigm
2568 !-------------------------------------------------------------------------------
2569       subroutine fragment_list
2570 !      implicit none
2571 !      include 'DIMENSIONS'
2572 !      include 'DIMENSIONS.ZSCOPT'
2573 !      include 'DIMENSIONS.COMPAR'
2574 !      include 'COMMON.IOUNITS'
2575 !      include 'COMMON.COMPAR'
2576       logical :: lprn=.true.
2577       integer :: i,ilevel,j,k,jfrag
2578       do jfrag=1,nfrag(1)
2579         nlist_frag(jfrag)=0
2580         do i=1,npiece(jfrag,1)
2581           if (lprn) write (iout,*) "jfrag=",jfrag,&
2582             "i=",i," fragment",ifrag(1,i,jfrag),&
2583             ifrag(2,i,jfrag)
2584           do j=ifrag(1,i,jfrag),ifrag(2,i,jfrag)
2585             do k=1,nlist_frag(jfrag)
2586               if (list_frag(k,jfrag).eq.j) goto 10
2587             enddo
2588             nlist_frag(jfrag)=nlist_frag(jfrag)+1
2589             list_frag(nlist_frag(jfrag),jfrag)=j
2590           enddo
2591   10      continue
2592         enddo
2593       enddo
2594       write (iout,*) "Fragment list"
2595       do j=1,nfrag(1)
2596         write (iout,*)"Fragment",j," list",(list_frag(k,j),&
2597          k=1,nlist_frag(j))
2598       enddo
2599       return
2600       end subroutine fragment_list
2601 !-------------------------------------------------------------------------------
2602       real(kind=8) function rmscalc(ishif,i,j,jcon,lprn)
2603
2604       use w_comm_local
2605       use control_data, only:symetr
2606       use geometry_data, only:nperm
2607 !      implicit real*8 (a-h,o-z)
2608 !      include 'DIMENSIONS'
2609 !      include 'DIMENSIONS.ZSCOPT'
2610 !      include 'DIMENSIONS.COMPAR'
2611 !      include 'COMMON.IOUNITS'
2612 !      include 'COMMON.COMPAR'
2613 !      include 'COMMON.CHAIN'
2614 !      include 'COMMON.INTERACT'
2615 !      include 'COMMON.VAR'
2616 !      include 'COMMON.CONTROL'
2617       real(kind=8) :: przes(3),obrot(3,3)
2618 !el      real(kind=8) :: creff(3,nres*2),cc(3,nres*2)
2619 !el      logical :: iadded(nres)
2620 !el      integer :: inumber(2,nres)
2621 !el      common /ccc/ creff,cc,iadded,inumber
2622       logical :: lprn
2623       logical :: non_conv
2624       integer :: ishif,i,j,jcon,idup,kkk,l,k,kk
2625       real(kind=8) :: rminrms,rms
2626       if (lprn) then
2627         write (iout,*) "i",i," j",j," jcont",jcon," ishif",ishif
2628         write (iout,*) "npiece",npiece(j,i)
2629         call flush(iout)
2630       endif
2631 !      write (iout,*) "symetr",symetr
2632 !      call flush(iout)
2633       nperm=1
2634       do idup=1,symetr
2635       nperm=nperm*idup
2636       enddo
2637 !      write (iout,*) "nperm",nperm
2638 !      call flush(iout)
2639       do kkk=1,nperm
2640       idup=0
2641       do l=1,nres
2642         iadded(l)=.false.
2643       enddo
2644 !      write (iout,*) "kkk",kkk
2645 !      call flush(iout)
2646       do k=1,npiece(j,i)
2647         if (i.eq.1) then
2648           if (lprn) then
2649             write (iout,*) "Level 1: j=",j,"k=",k," adding fragment",&
2650                ifrag(1,k,j),ifrag(2,k,j)
2651             call flush(iout)
2652           endif
2653           call cprep(ifrag(1,k,j),ifrag(2,k,j),ishif,idup,kkk)
2654 !          write (iout,*) "Exit cprep"
2655 !          call flush(iout)
2656 !          write (iout,*) "ii=",ii
2657         else
2658           kk = ipiece(k,j,i)
2659 !          write (iout,*) "kk",kk," npiece",npiece(kk,1)
2660           do l=1,npiece(kk,1)
2661             if (lprn) then
2662               write (iout,*) "Level",i,": j=",j,"k=",k," kk=",kk,&
2663                 " l=",l," adding fragment",&
2664                 ifrag(1,l,kk),ifrag(2,l,kk)
2665               call flush(iout)
2666             endif
2667             call cprep(ifrag(1,l,kk),ifrag(2,l,kk),ishif,idup,kkk)
2668 !            write (iout,*) "After cprep"
2669 !            call flush(iout)
2670           enddo
2671         endif
2672       enddo
2673       enddo
2674       if (lprn) then
2675         write (iout,*) "tuszukaj"
2676         do kkk=1,nperm
2677           do k=1,idup
2678             write(iout,'(5i4,2(3f10.5,5x))') i,j,k,inumber(1,k),&
2679               inumber(2,k),(creff(l,k),l=1,3),(cc(l,k),l=1,3)
2680           enddo
2681         enddo
2682         call flush(iout)
2683       endif
2684       rminrms=1.0d10
2685       do kkk=1,nperm
2686       call fitsq(rms,cc(1,1),creff(1,1),idup,przes,obrot,non_conv)
2687       if (non_conv) then
2688         print *,'Error: FITSQ non-convergent, jcon',jcon,i
2689         rms = 1.0d10
2690       else if (rms.lt.-1.0d-6) then
2691         print *,'Error: rms^2 = ',rms,jcon,i
2692         rms = 1.0d10
2693       else if (rms.ge.1.0d-6 .and. rms.lt.0) then
2694         rms = 0.0d0
2695       endif
2696 !      write (iout,*) "rmsmin", rminrms, "rms", rms
2697       if (rms.le.rminrms) rminrms=rms
2698       enddo
2699       rmscalc = dsqrt(rminrms)
2700 !      write (iout, *) "analysys", rmscalc,anatemp
2701       return
2702       end function rmscalc
2703 !-------------------------------------------------------------------------
2704       subroutine cprep(if1,if2,ishif,idup,kwa)
2705
2706       use w_comm_local
2707       use control_data, only:symetr
2708       use geometry_data, only:nperm,cref,c
2709 !      implicit real*8 (a-h,o-z)
2710 !      include 'DIMENSIONS'
2711 !      include 'DIMENSIONS.ZSCOPT'
2712 !      include 'DIMENSIONS.COMPAR'
2713 !      include 'COMMON.CONTROL'
2714 !      include 'COMMON.IOUNITS'
2715 !      include 'COMMON.COMPAR'
2716 !      include 'COMMON.CHAIN'
2717 !      include 'COMMON.INTERACT'
2718 !      include 'COMMON.VAR'
2719       real(kind=8) :: przes(3),obrot(3,3)
2720 !el      real(kind=8) :: creff(3,nres*2),cc(3,nres*2)
2721 !el      logical :: iadded(nres)
2722 !el      integer :: inumber(2,nres)
2723       integer :: iistrart,kwa,blar
2724 !el      common /ccc/ creff,cc,iadded,inumber
2725       integer :: if1,if2,ishif,idup,kkk,l,m
2726 !      write (iout,*) "Calling cprep symetr",symetr," kwa",kwa
2727       nperm=1
2728       do blar=1,symetr
2729       nperm=nperm*blar
2730       enddo
2731 !      write (iout,*) "nperm",nperm
2732       kkk=kwa
2733 !      ii=0
2734       do l=if1,if2
2735 !        write (iout,*) "l",l," iadded",iadded(l)
2736 !        call flush(iout)
2737         if (l+ishif.gt.1 .and. l+ishif.le.nres .and. .not.iadded(l)) &
2738         then
2739           idup=idup+1
2740           iadded(l)=.true.
2741           inumber(1,idup)=l
2742           inumber(2,idup)=l+ishif
2743           do m=1,3
2744             creff(m,idup)=cref(m,l,kkk)
2745             cc(m,idup)=c(m,l+ishif)
2746           enddo
2747         endif
2748       enddo
2749       return
2750       end subroutine cprep
2751 !-------------------------------------------------------------------------
2752       real(kind=8) function rmsnat(jcon)
2753
2754       use control_data, only:symetr
2755       use geometry_data, only:nperm,cref,c
2756       use energy_data, only:itype,molnum
2757 !      implicit real*8 (a-h,o-z)
2758 !      include 'DIMENSIONS'
2759 !      include 'DIMENSIONS.ZSCOPT'
2760 !      include 'DIMENSIONS.COMPAR'
2761 !      include 'COMMON.IOUNITS'
2762 !      include 'COMMON.COMPAR'
2763 !      include 'COMMON.CHAIN'
2764 !      include 'COMMON.INTERACT'
2765 !      include 'COMMON.VAR'
2766 !      include 'COMMON.CONTROL'
2767       real(kind=8) :: przes(3),obrot(3,3),cc(3,2*nres),ccref(3,2*nres)
2768       logical :: non_conv
2769       integer :: ishif,i,j,resprzesun,jcon,kkk,nnsup
2770       real(kind=8) :: rminrms,rmsminsing,rms
2771       rminrms=10.0d10
2772       rmsminsing=10d10
2773       nperm=1
2774       do i=1,symetr
2775        nperm=nperm*i
2776       enddo
2777       do kkk=1,nperm
2778        nnsup=0
2779        do i=1,nres
2780         if (itype(i,molnum(i)).ne.ntyp1_molec(molnum(i))) then
2781           nnsup=nnsup+1
2782           do j=1,3
2783             cc(j,nnsup)=c(j,i)
2784             ccref(j,nnsup)=cref(j,i,kkk)
2785           enddo
2786         endif
2787        enddo
2788        call fitsq(rms,cc(1,1),ccref(1,1),nnsup,przes,obrot,non_conv)
2789        if (non_conv) then
2790         print *,'Error: FITSQ non-convergent, jcon',jcon,i
2791         rms=1.0d10
2792        else if (rms.lt.-1.0d-6) then
2793         print *,'Error: rms^2 = ',rms,jcon,i
2794         rms = 1.0d10
2795        else if (rms.ge.1.0d-6 .and. rms.lt.0) then
2796         rms=0.0d0
2797        endif
2798        if (rms.le.rminrms) rminrms=rms
2799 !       write (iout,*) "kkk",kkk," rmsnat",rms , rminrms
2800       enddo
2801       rmsnat = dsqrt(rminrms)
2802 !      write (iout,*)  "analysys",rmsnat, anatemp
2803 !      liczenie rmsdla pojedynczego lancucha
2804       return
2805       end function rmsnat
2806 !-------------------------------------------------------------------------------
2807       subroutine define_fragments
2808
2809       use geometry_data, only:rad2deg
2810       use energy_data, only:itype,molnum
2811       use compare_data, only:nhfrag,nbfrag,bfrag,hfrag
2812 !      implicit real*8 (a-h,o-z)
2813 !      include 'DIMENSIONS'
2814 !      include 'DIMENSIONS.ZSCOPT'
2815 !      include 'DIMENSIONS.COMPAR'
2816 !      include 'COMMON.IOUNITS'
2817 !      include 'COMMON.TIME1'
2818 !      include 'COMMON.FRAG'
2819 !      include 'COMMON.SBRIDGE'
2820 !      include 'COMMON.CONTROL'
2821 !      include 'COMMON.COMPAR'
2822 !      include 'COMMON.CHAIN'
2823 !      include 'COMMON.HEADER'
2824 !      include 'COMMON.GEO'
2825 !      include 'COMMON.CONTACTS'
2826 !      include 'COMMON.PEPTCONT'
2827 !      include 'COMMON.INTERACT'
2828 !      include 'COMMON.NAMES'
2829       integer :: nstrand,istrand(2,nres/2)
2830       integer :: nhairp,ihairp(2,nres/5),mnum
2831       character(len=16) :: strstr(4)=reshape((/'helix','hairpin',&
2832                           'strand','strand pair'/),shape(strstr))
2833       integer :: j,i,ii,i1,i2,i3,i4,it1,it2,it3,it4
2834
2835       write (iout,*) 'NC_FRAC_HEL',ncfrac_hel,' NC_REQ_HEL',ncreq_hel,&
2836                      'NC_FRAC_BET',ncfrac_bet,' NC_REQ_BET',ncreq_bet,&
2837                  'NC_FRAC_PAIR',ncfrac_pair,' NC_REQ_PAIR',ncreq_pair,&
2838         ' RMS_PAIR',irms_pair,' SPLIT_BET',isplit_bet
2839       write (iout,*) 'NSHIFT_HEL',nshift_hel,' NSHIFT_BET',nshift_bet,&
2840         ' NSHIFT_STRAND',nshift_strand,' NSHIFT_PAIR',nshift_pair
2841       write (iout,*) 'ANGCUT_HEL',angcut_hel*rad2deg,&
2842         ' MAXANG_HEL',angcut1_hel*rad2deg
2843       write (iout,*) 'ANGCUT_BET',angcut_bet*rad2deg,&
2844                      ' MAXANG_BET',angcut1_bet*rad2deg
2845       write (iout,*) 'ANGCUT_STRAND',angcut_strand*rad2deg,&
2846                      ' MAXANG_STRAND',angcut1_strand*rad2deg
2847       write (iout,*) 'FRAC_MIN',frac_min_set
2848 ! Find secondary structure elements (helices and beta-sheets)
2849       call secondary2(.true.,.false.,ncont_pept_ref,icont_pept_ref,&
2850          isec_ref)
2851 ! Define primary fragments. First include the helices.
2852       nhairp=0
2853       nstrand=0
2854 ! Merge helices
2855 ! AL 12/23/03 - to avoid splitting helices into very small fragments
2856       if (merge_helices) then
2857       write (iout,*) "Before merging helices: nhfrag",nhfrag
2858       do i=1,nhfrag
2859         write (2,*) hfrag(1,i),hfrag(2,i)
2860       enddo
2861       i=1
2862       do while (i.lt.nhfrag)
2863         if (hfrag(1,i+1)-hfrag(2,i).le.1) then
2864           nhfrag=nhfrag-1
2865           hfrag(2,i)=hfrag(2,i+1)
2866           do j=i+1,nhfrag
2867             hfrag(1,j)=hfrag(1,j+1)
2868             hfrag(2,j)=hfrag(2,j+1)
2869           enddo
2870         endif
2871         i=i+1
2872       enddo
2873       write (iout,*) "After merging helices: nhfrag",nhfrag
2874       do i=1,nhfrag
2875         write (2,*) hfrag(1,i),hfrag(2,i)
2876       enddo
2877       endif
2878       nfrag(1)=nhfrag
2879       do i=1,nhfrag
2880         npiece(i,1)=1
2881         ifrag(1,1,i)=hfrag(1,i)
2882         ifrag(2,1,i)=hfrag(2,i)
2883         n_shift(1,i,1)=0
2884         n_shift(2,i,1)=nshift_hel
2885         ang_cut(i)=angcut_hel
2886         ang_cut1(i)=angcut1_hel
2887         frac_min(i)=frac_min_set
2888         nc_fragm(i,1)=ncfrac_hel
2889         nc_req_setf(i,1)=ncreq_hel
2890         istruct(i)=1
2891       enddo
2892       write (iout,*) "isplit_bet",isplit_bet
2893       if (isplit_bet.gt.1) then
2894 ! Split beta-sheets into strands and store strands as primary fragments.
2895         call split_beta(nbfrag,bfrag,nstrand,istrand,nhairp,ihairp)
2896         do i=1,nstrand
2897           ii=i+nfrag(1)
2898           npiece(ii,1)=1
2899           ifrag(1,1,ii)=istrand(1,i)
2900           ifrag(2,1,ii)=istrand(2,i)
2901           n_shift(1,ii,1)=nshift_strand
2902           n_shift(2,ii,1)=nshift_strand
2903           ang_cut(ii)=angcut_strand
2904           ang_cut1(ii)=angcut1_strand
2905           frac_min(ii)=frac_min_set
2906           nc_fragm(ii,1)=0
2907           nc_req_setf(ii,1)=0
2908           istruct(ii)=3
2909         enddo
2910         nfrag(1)=nfrag(1)+nstrand
2911       else if (isplit_bet.eq.1) then
2912 ! Split only far beta-sheets; does not split hairpins.
2913         call find_and_remove_hairpins(nbfrag,bfrag,nhairp,ihairp)
2914         call split_beta(nbfrag,bfrag,nstrand,istrand,nhairp,ihairp)
2915         do i=1,nhairp
2916           ii=i+nfrag(1)
2917           npiece(ii,1)=1
2918           ifrag(1,1,ii)=ihairp(1,i)
2919           ifrag(2,1,ii)=ihairp(2,i)
2920           n_shift(1,ii,1)=nshift_bet
2921           n_shift(2,ii,1)=nshift_bet
2922           ang_cut(ii)=angcut_bet
2923           ang_cut1(ii)=angcut1_bet
2924           frac_min(ii)=frac_min_set
2925           nc_fragm(ii,1)=ncfrac_bet
2926           nc_req_setf(ii,1)=ncreq_bet
2927           istruct(ii)=2
2928         enddo
2929         nfrag(1)=nfrag(1)+nhairp
2930         do i=1,nstrand
2931           ii=i+nfrag(1)
2932           npiece(ii,1)=1
2933           ifrag(1,1,ii)=istrand(1,i)
2934           ifrag(2,1,ii)=istrand(2,i)
2935           n_shift(1,ii,1)=nshift_strand
2936           n_shift(2,ii,1)=nshift_strand
2937           ang_cut(ii)=angcut_strand
2938           ang_cut1(ii)=angcut1_strand
2939           frac_min(ii)=frac_min_set
2940           nc_fragm(ii,1)=0
2941           nc_req_setf(ii,1)=0
2942           istruct(ii)=3
2943         enddo
2944         nfrag(1)=nfrag(1)+nstrand
2945       else
2946 ! Do not split beta-sheets; each pair of strands is a primary element.
2947         call find_and_remove_hairpins(nbfrag,bfrag,nhairp,ihairp)
2948         do i=1,nhairp
2949           ii=i+nfrag(1)
2950           npiece(ii,1)=1
2951           ifrag(1,1,ii)=ihairp(1,i)
2952           ifrag(2,1,ii)=ihairp(2,i)
2953           n_shift(1,ii,1)=nshift_bet
2954           n_shift(2,ii,1)=nshift_bet
2955           ang_cut(ii)=angcut_bet
2956           ang_cut1(ii)=angcut1_bet
2957           frac_min(ii)=frac_min_set
2958           nc_fragm(ii,1)=ncfrac_bet
2959           nc_req_setf(ii,1)=ncreq_bet
2960           istruct(ii)=2
2961         enddo
2962         nfrag(1)=nfrag(1)+nhairp
2963         do i=1,nbfrag
2964           ii=i+nfrag(1)
2965           npiece(ii,1)=2
2966           ifrag(1,1,ii)=bfrag(1,i)
2967           ifrag(2,1,ii)=bfrag(2,i)
2968           if (bfrag(3,i).lt.bfrag(4,i)) then
2969             ifrag(1,2,ii)=bfrag(3,i)
2970             ifrag(2,2,ii)=bfrag(4,i)
2971           else
2972             ifrag(1,2,ii)=bfrag(4,i)
2973             ifrag(2,2,ii)=bfrag(3,i)
2974           endif
2975           n_shift(1,ii,1)=nshift_bet
2976           n_shift(2,ii,1)=nshift_bet
2977           ang_cut(ii)=angcut_bet
2978           ang_cut1(ii)=angcut1_bet
2979           frac_min(ii)=frac_min_set
2980           nc_fragm(ii,1)=ncfrac_bet
2981           nc_req_setf(ii,1)=ncreq_bet
2982           istruct(ii)=4
2983         enddo
2984         nfrag(1)=nfrag(1)+nbfrag
2985       endif
2986       write (iout,*) "The following primary fragments were found:"
2987       write (iout,*) "Helices:",nhfrag
2988       do i=1,nhfrag
2989         mnum=molnum(i)
2990         i1=ifrag(1,1,i)
2991         i2=ifrag(2,1,i)
2992         it1=itype(i1,mnum)
2993         it2=itype(i2,mnum)
2994         write (iout,'(i3,2x,a,i4,2x,a,i4)') &
2995              i,restyp(it1,mnum),i1,restyp(it2,mnum),i2
2996       enddo
2997       write (iout,*) "Hairpins:",nhairp
2998       do i=nhfrag+1,nhfrag+nhairp
2999         i1=ifrag(1,1,i)
3000         i2=ifrag(2,1,i)
3001         it1=itype(i1,molnum(i1))
3002         it2=itype(i2,molnum(i2))
3003         write (iout,'(i3,2x,a,i4,2x,a,i4,2x)') &
3004              i,restyp(it1,molnum(i1)),i1,restyp(it2,molnum(i2)),i2
3005       enddo
3006       write (iout,*) "Far strand pairs:",nbfrag
3007       do i=nhfrag+nhairp+1,nhfrag+nhairp+nbfrag
3008         i1=ifrag(1,1,i)
3009         i2=ifrag(2,1,i)
3010         it1=itype(i1,molnum(i1))
3011         it2=itype(i2,molnum(i1))
3012         i3=ifrag(1,2,i)
3013         i4=ifrag(2,2,i)
3014         it3=itype(i3,molnum(i3))
3015         it4=itype(i4,molnum(i4))
3016         write (iout,'(i3,2x,a,i4,2x,a,i4," and ",a,i4,2x,a,i4)') &
3017              i,restyp(it1,molnum(i1)),i1,restyp(it2,molnum(i2)),i2,&
3018                restyp(it3,molnum(i3)),i3,restyp(it4,molnum(i4)),i4
3019       enddo
3020       write (iout,*) "Strands:",nstrand
3021       do i=nhfrag+nhairp+nbfrag+1,nfrag(1)
3022         mnum=molnum(i)
3023         i1=ifrag(1,1,i)
3024         i2=ifrag(2,1,i)
3025         it1=itype(i1,mnum)
3026         it2=itype(i2,mnum)
3027         write (iout,'(i3,2x,a,i4,2x,a,i4)') &
3028              i,restyp(it1,mnum),i1,restyp(it2,mnum),i2
3029       enddo
3030       call imysort(nfrag(1),2,maxpiece,ifrag(1,1,1),npiece(1,1),&
3031         istruct(1),n_shift(1,1,1),ang_cut(1),ang_cut1(1),frac_min(1),&
3032         nc_fragm(1,1),nc_req_setf(1,1))
3033       write (iout,*) "Fragments after sorting:"
3034       do i=1,nfrag(1)
3035         mnum=molnum(i)
3036         i1=ifrag(1,1,i)
3037         i2=ifrag(2,1,i)
3038         it1=itype(i1,mnum)
3039         it2=itype(i2,mnum)
3040         write (iout,'(i3,2x,a,i4,2x,a,i4,$)') &
3041              i,restyp(it1,molnum(it1)),i1,restyp(it2,molnum(it2)),i2
3042         if (npiece(i,1).eq.1) then
3043           write (iout,'(2x,a)') strstr(istruct(i))
3044         else
3045           i1=ifrag(1,2,i)
3046           i2=ifrag(2,2,i)
3047           it1=itype(i1,mnum)
3048           it2=itype(i2,mnum)
3049           write (iout,'(2x,a,i4,2x,a,i4,2x,a)') &
3050              restyp(it1,molnum(it1)),i1,restyp(it2,molnum(it2)),i2,strstr(istruct(i))
3051         endif
3052       enddo
3053       return
3054       end subroutine define_fragments
3055 !------------------------------------------------------------------------------
3056       subroutine find_and_remove_hairpins(nbfrag,bfrag,nhairp,ihairp)
3057 !      implicit real*8 (a-h,o-z)
3058 !      include 'DIMENSIONS'
3059 !      include 'DIMENSIONS.ZSCOPT'
3060 !      include 'DIMENSIONS.COMPAR'
3061 !      include 'COMMON.IOUNITS'
3062       integer :: nbfrag,bfrag(4,nres/3)
3063       integer :: nhairp,ihairp(2,nres/5)
3064       integer :: i,j,k
3065       write (iout,*) "Entered find_and_remove_hairpins"
3066       write (iout,*) "nbfrag",nbfrag
3067       do i=1,nbfrag
3068         write (iout,*) i,(bfrag(k,i),k=1,4)
3069       enddo
3070       nhairp=0
3071       i=1
3072       do while (i.le.nbfrag)
3073         write (iout,*) "check hairpin:",i,(bfrag(j,i),j=1,4)
3074         if (bfrag(3,i).gt.bfrag(4,i) .and. bfrag(4,i)-bfrag(2,i).lt.5) &
3075         then
3076           write (iout,*) "Found hairpin:",i,bfrag(1,i),bfrag(3,i)
3077           nhairp=nhairp+1
3078           ihairp(1,nhairp)=bfrag(1,i)
3079           ihairp(2,nhairp)=bfrag(3,i)
3080           nbfrag=nbfrag-1
3081           do j=i,nbfrag
3082             do k=1,4
3083               bfrag(k,j)=bfrag(k,j+1)
3084             enddo
3085           enddo
3086         else
3087           i=i+1
3088         endif
3089       enddo
3090       write (iout,*) "After finding hairpins:"
3091       write (iout,*) "nhairp",nhairp
3092       do i=1,nhairp
3093         write (iout,*) i,ihairp(1,i),ihairp(2,i)
3094       enddo
3095       write (iout,*) "nbfrag",nbfrag
3096       do i=1,nbfrag
3097         write (iout,*) i,(bfrag(k,i),k=1,4)
3098       enddo
3099       return
3100       end subroutine find_and_remove_hairpins
3101 !------------------------------------------------------------------------------
3102       subroutine split_beta(nbfrag,bfrag,nstrand,istrand,nhairp,ihairp)
3103 !      implicit real*8 (a-h,o-z)
3104 !      include 'DIMENSIONS'
3105 !      include 'DIMENSIONS.ZSCOPT'
3106 !      include 'DIMENSIONS.COMPAR'
3107 !      include 'COMMON.IOUNITS'
3108       integer :: nbfrag,bfrag(4,nres/3)
3109       integer :: nstrand,istrand(2,nres/2)
3110       integer :: nhairp,ihairp(2,nres/5)
3111       logical :: found
3112       integer :: i,k
3113       write (iout,*) "Entered split_beta"
3114       write (iout,*) "nbfrag",nbfrag
3115       do i=1,nbfrag
3116         write (iout,*) i,(bfrag(k,i),k=1,4)
3117       enddo
3118       nstrand=0
3119       do i=1,nbfrag
3120         write (iout,*) "calling add_strand:",i,bfrag(1,i),bfrag(2,i)
3121         call add_strand(nstrand,istrand,nhairp,ihairp,&
3122            bfrag(1,i),bfrag(2,i),found)
3123         if (bfrag(3,i).lt.bfrag(4,i)) then
3124           write (iout,*) "calling add_strand:",i,bfrag(3,i),bfrag(4,i)
3125           call add_strand(nstrand,istrand,nhairp,ihairp,&
3126            bfrag(3,i),bfrag(4,i),found)
3127         else
3128           write (iout,*) "calling add_strand:",i,bfrag(4,i),bfrag(3,i)
3129           call add_strand(nstrand,istrand,nhairp,ihairp,&
3130             bfrag(4,i),bfrag(3,i),found)
3131         endif
3132       enddo
3133       nbfrag=0
3134       write (iout,*) "Strands found:",nstrand
3135       do i=1,nstrand
3136         write (iout,*) i,istrand(1,i),istrand(2,i)
3137       enddo
3138       return
3139       end subroutine split_beta
3140 !------------------------------------------------------------------------------
3141       subroutine add_strand(nstrand,istrand,nhairp,ihairp,is1,is2,found)
3142 !      implicit real*8 (a-h,o-z)
3143 !      include 'DIMENSIONS'
3144 !      include 'DIMENSIONS.ZSCOPT'
3145 !      include 'DIMENSIONS.COMPAR'
3146 !      include 'COMMON.IOUNITS'
3147       integer :: nstrand,istrand(2,nres/2)
3148       integer :: nhairp,ihairp(2,nres/5)
3149       logical :: found
3150       integer :: is1,is2,j,idelt
3151       found=.false.
3152       do j=1,nhairp
3153         idelt=(ihairp(2,j)-ihairp(1,j))/6
3154         if (is1.lt.ihairp(2,j)-idelt.and.is2.gt.ihairp(1,j)+idelt) then
3155           write (iout,*) "strand",is1,is2," is part of hairpin",&
3156             ihairp(1,j),ihairp(2,j)
3157           return
3158         endif
3159       enddo
3160       do j=1,nstrand
3161         idelt=(istrand(2,j)-istrand(1,j))/3
3162         if (is1.lt.istrand(2,j)-idelt.and.is2.gt.istrand(1,j)+idelt) &
3163         then
3164 ! The strand already exists in the array; update its ends if necessary.
3165           write (iout,*) "strand",is1,is2," found at position",j,&
3166            ":",istrand(1,j),istrand(2,j)
3167           istrand(1,j)=min0(istrand(1,j),is1)
3168           istrand(2,j)=max0(istrand(2,j),is2)
3169           return
3170         endif
3171       enddo
3172 ! The strand has not been found; add it to the array.
3173       write (iout,*) "strand",is1,is2," added to the array."
3174       found=.true.
3175       nstrand=nstrand+1
3176       istrand(1,nstrand)=is1
3177       istrand(2,nstrand)=is2
3178       return
3179       end subroutine add_strand
3180 !------------------------------------------------------------------------------
3181       subroutine secondary2(lprint,lprint_sec,ncont,icont,isecstr)
3182
3183       use geometry_data, only:anatemp,rad2deg,phi,nstart_sup,nend_sup
3184       use energy_data, only:itype,maxcont,molnum
3185       use compare_data, only:bfrag,hfrag,nbfrag,nhfrag
3186       use compare, only:freeres
3187 !      implicit real*8 (a-h,o-z)
3188 !      include 'DIMENSIONS'
3189 !      include 'DIMENSIONS.ZSCOPT'
3190 !      include 'COMMON.IOUNITS'
3191 !      include 'COMMON.FRAG'
3192 !      include 'COMMON.VAR'
3193 !      include 'COMMON.GEO'
3194 !      include 'COMMON.CHAIN'
3195 !      include 'COMMON.NAMES'
3196 !      include 'COMMON.INTERACT'
3197       integer :: ncont,icont(2,maxcont),isec(nres,4),nsec(nres+1),&
3198         isecstr(nres+1)
3199       logical :: lprint,lprint_sec,not_done !el,freeres
3200       integer :: i,j,ii1,jj1,i1,j1,ij,k,ien,ist
3201       integer :: nstrand,nbeta,nhelix,iii1,jjj1,mnum
3202       real(kind=8) :: p1,p2
3203 !rel      external freeres
3204       character(len=1) :: csec(0:2)=reshape((/'-','E','H'/),shape(csec))
3205       if (lprint) then
3206         write (iout,*) "entered secondary2",ncont
3207         write (iout,*) "nstart_sup",nstart_sup," nend_sup",nend_sup
3208         do i=1,ncont
3209           write (iout,*) icont(1,i),icont(2,i)
3210         enddo
3211       endif
3212       do i=1,nres
3213         isecstr(i)=0
3214       enddo
3215       nbfrag=0
3216       nhfrag=0
3217       do i=1,nres
3218         isec(i,1)=0
3219         isec(i,2)=0
3220         nsec(i)=0
3221       enddo
3222
3223 ! finding parallel beta
3224 !d      write (iout,*) '------- looking for parallel beta -----------'
3225       nbeta=0
3226       nstrand=0
3227       do i=1,ncont
3228         i1=icont(1,i)
3229         j1=icont(2,i)
3230         if (i1.ge.nstart_sup .and. i1.le.nend_sup &
3231            .and. j1.gt.nstart_sup .and. j1.le.nend_sup) then
3232 !d        write (iout,*) "parallel",i1,j1
3233         if(j1-i1.gt.5 .and. freeres(i1,j1,nsec,isec)) then
3234           ii1=i1
3235           jj1=j1
3236 !d          write (iout,*) i1,j1
3237           not_done=.true.
3238           do while (not_done)
3239            i1=i1+1
3240            j1=j1+1
3241             do j=1,ncont
3242               if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j) .and. &
3243                    freeres(i1,j1,nsec,isec)) goto 5
3244             enddo
3245             not_done=.false.
3246   5         continue
3247 !d            write (iout,*) i1,j1,not_done
3248           enddo
3249           j1=j1-1
3250           i1=i1-1
3251           if (i1-ii1.gt.1) then
3252             ii1=max0(ii1-1,1)
3253             jj1=max0(jj1-1,1)
3254             nbeta=nbeta+1
3255             if(lprint)write(iout,'(a,i3,4i4)')'parallel beta',&
3256                      nbeta,ii1,i1,jj1,j1
3257
3258             nbfrag=nbfrag+1
3259             bfrag(1,nbfrag)=ii1+1
3260             bfrag(2,nbfrag)=i1+1
3261             bfrag(3,nbfrag)=jj1+1
3262             bfrag(4,nbfrag)=min0(j1+1,nres)
3263
3264             do ij=ii1,i1
3265              nsec(ij)=nsec(ij)+1
3266              isec(ij,nsec(ij))=nbeta
3267             enddo
3268             do ij=jj1,j1
3269              nsec(ij)=nsec(ij)+1
3270              isec(ij,nsec(ij))=nbeta
3271             enddo
3272
3273            if(lprint_sec) then
3274             nstrand=nstrand+1
3275             if (nbeta.le.9) then
3276               write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
3277                 "DefPropRes 'strand",nstrand,&
3278                 "' 'num = ",ii1-1,"..",i1-1,"'"
3279             else
3280               write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
3281                 "DefPropRes 'strand",nstrand,&
3282                 "' 'num = ",ii1-1,"..",i1-1,"'"
3283             endif
3284             nstrand=nstrand+1
3285             if (nbeta.le.9) then
3286               write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
3287                 "DefPropRes 'strand",nstrand,&
3288                 "' 'num = ",jj1-1,"..",j1-1,"'"
3289             else
3290               write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
3291                 "DefPropRes 'strand",nstrand,&
3292                 "' 'num = ",jj1-1,"..",j1-1,"'"
3293             endif
3294               write(12,'(a8,4i4)') &
3295                 "SetNeigh",ii1-1,i1-1,jj1-1,j1-1
3296            endif
3297           endif
3298         endif
3299         endif ! i1.ge.nstart_sup .and. i1.le.nend_sup .and. i2.gt.nstart_sup .and. i2.le.nend_sup
3300       enddo
3301
3302 ! finding antiparallel beta
3303 !d      write (iout,*) '--------- looking for antiparallel beta ---------'
3304
3305       do i=1,ncont
3306         i1=icont(1,i)
3307         j1=icont(2,i)
3308         if (freeres(i1,j1,nsec,isec)) then
3309           ii1=i1
3310           jj1=j1
3311 !d          write (iout,*) i1,j1
3312
3313           not_done=.true.
3314           do while (not_done)
3315            i1=i1+1
3316            j1=j1-1
3317             do j=1,ncont
3318               if (i1.eq.icont(1,j).and.j1.eq.icont(2,j) .and. &
3319                    freeres(i1,j1,nsec,isec)) goto 6
3320             enddo
3321             not_done=.false.
3322   6         continue
3323 !d            write (iout,*) i1,j1,not_done
3324           enddo
3325           i1=i1-1
3326           j1=j1+1
3327           if (i1-ii1.gt.1) then
3328
3329             nbfrag=nbfrag+1
3330             bfrag(1,nbfrag)=ii1
3331             bfrag(2,nbfrag)=min0(i1+1,nres)
3332             bfrag(3,nbfrag)=min0(jj1+1,nres)
3333             bfrag(4,nbfrag)=j1
3334
3335             nbeta=nbeta+1
3336             iii1=max0(ii1-1,1)
3337             do ij=iii1,i1
3338              nsec(ij)=nsec(ij)+1
3339              if (nsec(ij).le.2) then
3340               isec(ij,nsec(ij))=nbeta
3341              endif
3342             enddo
3343             jjj1=max0(j1-1,1)
3344             do ij=jjj1,jj1
3345              nsec(ij)=nsec(ij)+1
3346              if (nsec(ij).le.2) then
3347               isec(ij,nsec(ij))=nbeta
3348              endif
3349             enddo
3350
3351
3352            if (lprint_sec) then
3353             write (iout,'(a,i3,4i4)')'antiparallel beta',&
3354                          nbeta,ii1-1,i1,jj1,j1-1
3355             nstrand=nstrand+1
3356             if (nstrand.le.9) then
3357               write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
3358                 "DefPropRes 'strand",nstrand,&
3359                 "' 'num = ",ii1-2,"..",i1-1,"'"
3360             else
3361               write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
3362                 "DefPropRes 'strand",nstrand,&
3363                 "' 'num = ",ii1-2,"..",i1-1,"'"
3364             endif
3365             nstrand=nstrand+1
3366             if (nstrand.le.9) then
3367               write(12,'(a18,i1,a9,i3,a2,i3,a1)') &
3368                 "DefPropRes 'strand",nstrand,&
3369                 "' 'num = ",j1-2,"..",jj1-1,"'"
3370             else
3371               write(12,'(a18,i2,a9,i3,a2,i3,a1)') &
3372                 "DefPropRes 'strand",nstrand,&
3373                 "' 'num = ",j1-2,"..",jj1-1,"'"
3374             endif
3375               write(12,'(a8,4i4)') &
3376                 "SetNeigh",ii1-2,i1-1,jj1-1,j1-2
3377            endif
3378           endif
3379         endif
3380       enddo
3381
3382 !d      write (iout,*) "After beta:",nbfrag
3383 !d      do i=1,nbfrag
3384 !d        write (iout,*) (bfrag(j,i),j=1,4)
3385 !d      enddo
3386
3387       if (nstrand.gt.0.and.lprint_sec) then
3388         write(12,'(a27,$)') "DefPropRes 'sheet' 'strand1"
3389         do i=2,nstrand
3390          if (i.le.9) then
3391           write(12,'(a9,i1,$)') " | strand",i
3392          else
3393           write(12,'(a9,i2,$)') " | strand",i
3394          endif
3395         enddo
3396         write(12,'(a1)') "'"
3397       endif
3398
3399
3400 ! finding alpha or 310 helix
3401
3402       nhelix=0
3403       do i=1,ncont
3404         i1=icont(1,i)
3405         j1=icont(2,i)
3406         p1=phi(i1+2)*rad2deg
3407         p2=0.0
3408         if (j1+2.le.nres) p2=phi(j1+2)*rad2deg
3409
3410
3411         if (j1.eq.i1+3 .and. &
3412              ((p1.ge.10.and.p1.le.80).or.i1.le.2).and. &
3413              ((p2.ge.10.and.p2.le.80).or.j1.le.2.or.j1.ge.nres-3) )then
3414 !d          if (j1.eq.i1+3) write (iout,*) "found 1-4 ",i1,j1,p1,p2
3415 !o          if (j1.eq.i1+4) write (iout,*) "found 1-5 ",i1,j1,p1,p2
3416           ii1=i1
3417           jj1=j1
3418           if (nsec(ii1).eq.0) then
3419             not_done=.true.
3420           else
3421             not_done=.false.
3422           endif
3423           do while (not_done)
3424             i1=i1+1
3425             j1=j1+1
3426             do j=1,ncont
3427               if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j)) goto 10
3428             enddo
3429             not_done=.false.
3430   10        continue
3431             p1=phi(i1+2)*rad2deg
3432             p2=phi(j1+2)*rad2deg
3433             if (p1.lt.10.or.p1.gt.80.or.p2.lt.10.or.p2.gt.80) &
3434                                     not_done=.false.
3435
3436 !d           write (iout,*) i1,j1,not_done,p1,p2
3437           enddo
3438           j1=j1+1
3439           if (j1-ii1.gt.4) then
3440             nhelix=nhelix+1
3441 !d            write (iout,*)'helix',nhelix,ii1,j1
3442
3443             nhfrag=nhfrag+1
3444             hfrag(1,nhfrag)=ii1
3445             hfrag(2,nhfrag)=j1
3446
3447             do ij=ii1,j1
3448              nsec(ij)=-1
3449             enddo
3450            if (lprint_sec) then
3451             write (iout,'(a,i3,2i4)') "Helix",nhelix,ii1-1,j1-1
3452             if (nhelix.le.9) then
3453               write(12,'(a17,i1,a9,i3,a2,i3,a1)') &
3454                 "DefPropRes 'helix",nhelix,&
3455                 "' 'num = ",ii1-1,"..",j1-2,"'"
3456             else
3457               write(12,'(a17,i2,a9,i3,a2,i3,a1)') &
3458                 "DefPropRes 'helix",nhelix,&
3459                 "' 'num = ",ii1-1,"..",j1-2,"'"
3460             endif
3461            endif
3462           endif
3463         endif
3464       enddo
3465
3466       if (nhelix.gt.0.and.lprint_sec) then
3467         write(12,'(a26,$)') "DefPropRes 'helix' 'helix1"
3468         do i=2,nhelix
3469          if (nhelix.le.9) then
3470           write(12,'(a8,i1,$)') " | helix",i
3471          else
3472           write(12,'(a8,i2,$)') " | helix",i
3473          endif
3474         enddo
3475         write(12,'(a1)') "'"
3476       endif
3477
3478       if (lprint_sec) then
3479        write(12,'(a37)') "DefPropRes 'coil' '! (helix | sheet)'"
3480        write(12,'(a20)') "XMacStand ribbon.mac"
3481       endif
3482
3483       if (lprint) then
3484
3485         write(iout,*) 'UNRES seq:',anatemp
3486         do j=1,nbfrag
3487          write(iout,*) 'beta ',(bfrag(i,j),i=1,4)
3488         enddo
3489
3490         do j=1,nhfrag
3491          write(iout,*) 'helix ',(hfrag(i,j),i=1,2),anatemp
3492         enddo
3493
3494       endif
3495
3496       do j=1,nbfrag
3497         do k=min0(bfrag(1,j),bfrag(2,j)),max0(bfrag(1,j),bfrag(2,j))
3498           isecstr(k)=1
3499         enddo
3500         do k=min0(bfrag(3,j),bfrag(4,j)),max0(bfrag(3,j),bfrag(4,j))
3501           isecstr(k)=1
3502         enddo
3503       enddo
3504       do j=1,nhfrag
3505         do k=hfrag(1,j),hfrag(2,j)
3506           isecstr(k)=2
3507         enddo
3508       enddo
3509       if (lprint) then
3510         write (iout,*)
3511         write (iout,*) "Secondary structure"
3512         do i=1,nres,80
3513           mnum=molnum(i)
3514           ist=i
3515           ien=min0(i+79,nres)
3516           write (iout,*)
3517           write (iout,'(8(7x,i3))') (k,k=ist+9,ien,10)
3518           write (iout,'(80a1)') (onelet(itype(k,mnum)),k=ist,ien)
3519           write (iout,'(80a1)') (csec(isecstr(k)),k=ist,ien)
3520         enddo
3521         write (iout,*)
3522       endif
3523       return
3524       end subroutine secondary2
3525 !-------------------------------------------------
3526 !      logical function freeres(i,j,nsec,isec)
3527 !      include 'DIMENSIONS'
3528 !      integer :: isec(nres,4),nsec(nres)
3529 !      integer :: i,j,k,l
3530 !      freeres=.false.
3531 !
3532 !      if (nsec(i).gt.1.or.nsec(j).gt.1) return
3533 !      do k=1,nsec(i)
3534 !        do l=1,nsec(j)
3535 !          if (isec(i,k).eq.isec(j,l)) return
3536 !        enddo
3537 !      enddo
3538 !      freeres=.true.
3539 !      return
3540 !      end function freeres
3541 !-------------------------------------------------
3542        subroutine alloc_compar_arrays(nfrg,nlev)
3543
3544        use energy_data, only:maxcont
3545        use w_comm_local
3546        integer :: nfrg,nlev
3547
3548 !write(iout,*) "in alloc conpar arrays: nlevel=", nlevel," nfrag(1)=",nfrag(1)
3549 !------------------------
3550 ! commom.contacts
3551 !      common /contacts/
3552       allocate(nsccont_frag_ref(mmaxfrag)) !(mmaxfrag) !wham
3553       allocate(isccont_frag_ref(2,maxcont,mmaxfrag)) !(2,maxcont,mmaxfrag) !wham
3554 !------------------------
3555 ! COMMON.COMPAR
3556 !      common /compar/
3557       allocate(rmsfrag(nfrg,nlev+1),nc_fragm(nfrg,nlev+1)) !(maxfrag,maxlevel)
3558       allocate(qfrag(nfrg,2)) !(maxfrag,2)
3559       allocate(rmscutfrag(2,nfrg,nlev+1)) !(2,maxfrag,maxlevel)
3560       allocate(ang_cut(nfrg),ang_cut1(nfrg),frac_min(nfrg)) !(maxfrag)
3561       allocate(nc_req_setf(nfrg,nlev+1),npiece(nfrg,nlev+1),&
3562         ielecont(nfrg,nlev+1),isccont(nfrg,nlev+1),irms(nfrg,nlev+1),&
3563         ishifft(nfrg,nlev+1),len_frag(nfrg,nlev+1)) !(maxfrag,maxlevel)
3564       allocate(ncont_nat(2,nfrg,nlev+1))
3565       allocate(n_shift(2,nfrg,nlev+1)) !(2,maxfrag,maxlevel)
3566 !      allocate(nfrag(nlev)) !(maxlevel)
3567       allocate(isnfrag(nlev+2)) !(maxlevel+1)
3568       allocate(ifrag(2,maxpiece,nfrg)) !(2,maxpiece,maxfrag)
3569       allocate(ipiece(maxpiece,nfrg,2:nlev+1)) !(maxpiece,maxfrag,2:maxlevel)
3570       allocate(istruct(nfrg),iloc(nfrg),nlist_frag(nfrg)) !(maxfrag)
3571       allocate(iclass(nlev*nfrg,nlev+1)) !(maxlevel*maxfrag,maxlevel)
3572       allocate(list_frag(nres,nfrg)) !(maxres,maxfrag)
3573 !------------------------
3574 ! COMMON.PEPTCONT
3575 !      common /peptcont/
3576 !      integer,dimension(:,:),allocatable :: icont_pept_ref !(2,maxcont)
3577       allocate(ncont_frag_ref(mmaxfrag)) !(mmaxfrag)
3578       allocate(icont_frag_ref(2,maxcont,mmaxfrag)) !(2,maxcont,mmaxfrag)
3579 !      integer,dimension(:),allocatable :: isec_ref !(maxres)
3580 !------------------------
3581 !      module w_comm_local
3582 !      common /ccc/
3583       allocate(creff(3,2*nres),cc(3,2*nres)) !(3,nres*2)
3584       allocate(iadded(nres)) !(nres)
3585       allocate(inumber(2,nres)) !(2,nres)
3586
3587
3588 !-------------------------------------------------------------------------------
3589       end subroutine alloc_compar_arrays
3590 #endif
3591 !-------------------------------------------------------------------------------
3592       end module conform_compar
3593