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