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