source/cluster/wham/src-HCD-5D/wrtclust.f

   1       SUBROUTINE WRTCLUST(NCON,ICUT,PRINTANG,PRINTPDB,printmol2,ib)
   2       implicit real*8 (a-h,o-z)
   3       include 'DIMENSIONS'
   4       include 'sizesclu.dat'
   5       parameter (num_in_line=5)
   6       LOGICAL PRINTANG(max_cut)
   7       integer PRINTPDB(max_cut),printmol2(max_cut)
   8       include 'COMMON.CONTROL'
   9       include 'COMMON.HEADER'
  10       include 'COMMON.CHAIN'
  11       include 'COMMON.VAR'
  12       include 'COMMON.CLUSTER'
  13       include 'COMMON.IOUNITS'
  14       include 'COMMON.GEO'
  15       include 'COMMON.FREE'
  16       include 'COMMON.TEMPFAC'
  17       include 'COMMON.FFIELD'
  18       include 'COMMON.SBRIDGE'
  19       include 'COMMON.TORCNSTR'
  20       include 'COMMON.SAXS'
  21       CHARACTER*64 prefixp,NUMM,MUMM,EXTEN,extmol
  22       character*120 cfname
  23       character*8 ctemper
  24       DATA EXTEN /'.pdb'/,extmol /'.mol2'/,NUMM /'000'/,MUMM /'000'/
  25       external ilen
  26       logical viol_nmr
  27       integer ib,list_peak_viol(maxdim_cont)
  28       double precision Esaxs_all(maxgr),Pcalc_all(maxsaxs,maxgr)
  29
  30       do i=1,64
  31         cfname(i:i)=" "
  32       enddo
  33 c      print *,"calling WRTCLUST",ncon
  34 c      write (iout,*) "ICUT",icut," PRINTPDB ",PRINTPDB(icut)
  35       rewind 80
  36       call flush(iout)
  37       temper=1.0d0/(beta_h(ib)*1.987d-3)
  38       if (temper.lt.100.0d0) then
  39         write(ctemper,'(f3.0)') temper
  40         ctemper(3:3)=" "
  41       else if (temper.lt.1000.0) then
  42         write (ctemper,'(f4.0)') temper
  43         ctemper(4:4)=" "
  44       else
  45         write (ctemper,'(f5.0)') temper
  46         ctemper(5:5)=" "
  47       endif
  48
  49       do i=1,ncon*(ncon-1)/2
  50         read (80) diss(i)
  51       enddo
  52       close(80,status='delete')
  53 C
  54 C  PRINT OUT THE RESULTS OF CLUSTER ANALYSIS
  55 C
  56       ii1= index(intinname,'/')
  57       ii2=ii1
  58       ii1=ii1+1
  59       do while (ii2.gt.0)
  60         ii1=ii1+ii2
  61         ii2=index(intinname(ii1:),'/')
  62       enddo
  63       ii = ii1+index(intinname(ii1:),'.')-1
  64       if (ii.eq.0) then
  65         ii=ilen(intinname)
  66       else
  67         ii=ii-1
  68       endif
  69       prefixp=intinname(ii1:ii)
  70 cd    print *,icut,printang(icut),printpdb(icut),printmol2(icut)
  71 cd    print *,'ecut=',ecut
  72       WRITE (iout,100) NGR
  73       DO 19 IGR=1,NGR
  74       WRITE (iout,200) IGR,totfree_gr(igr)/beta_h(ib),LICZ(IGR)
  75       NRECORD=LICZ(IGR)/num_in_line
  76       IND1=1
  77       DO 63 IRECORD=1,NRECORD
  78       IND2=IND1+num_in_line-1
  79       WRITE (iout,300) (list_conf(NCONF(IGR,ICO)),
  80      &    totfree(NCONF(IGR,ICO))/beta_h(ib),ICO=IND1,IND2)
  81       IND1=IND2+1
  82    63 CONTINUE
  83       WRITE (iout,300) (list_conf(NCONF(IGR,ICO)),
  84      &   totfree(NCONF(IGR,ICO))/beta_h(ib),ICO=IND1,LICZ(IGR))
  85       IND1=1
  86       ICON=list_conf(NCONF(IGR,1))
  87 c      WRITE (iout,'(16F5.0)') (rad2deg*phiall(IND,icon),IND=4,nphi+3)
  88 C 12/8/93 Estimation of "diameters" of the subsequent families.
  89       ave_dim=0.0
  90       amax_dim=0.0
  91 c      write (iout,*) "ecut",ecut
  92       emin=totfree(nconf(igr,1))
  93       do i=2,licz(igr)
  94         ii=nconf(igr,i)
  95         if (totfree(ii)-emin .gt. ecut) goto 10
  96         do j=1,i-1
  97           jj=nconf(igr,j)
  98           if (jj.eq.1) exit
  99           if (ii.lt.jj) then
 100             ind=ioffset(ncon,ii,jj)
 101           else
 102             ind=ioffset(ncon,jj,ii)
 103           endif
 104 c          write (iout,*) " ncon",ncon,"i",i," j",j," ii",ii," jj",jj,
 105 c     &     " ind",ind," diss",diss(ind)
 106 c          call flush(iout)
 107           curr_dist=dabs(diss(ind)+0.0d0)
 108 c          write(iout,'(i10,4i4,f12.4)') ind,ii,jj,list_conf(ii),
 109 c     &      list_conf(jj),curr_dist
 110           if (curr_dist .gt. amax_dim) amax_dim=curr_dist
 111           ave_dim=ave_dim+curr_dist**2
 112         enddo
 113       enddo
 114    10 if (licz(igr) .gt. 1)
 115      & ave_dim=sqrt(ave_dim/(licz(igr)*(licz(igr)-1)/2))
 116       write (iout,'(/A,F8.1,A,F8.1)')
 117      & 'Max. distance in the family:',amax_dim,
 118      & '; average distance in the family:',ave_dim
 119       rmsave(igr)=0.0d0
 120       gdt_ts_ave(igr)=0.0d0
 121       gdt_ha_ave(igr)=0.0d0
 122       tmscore_ave(igr)=0.0d0
 123       qpart=0.0d0
 124       e1=totfree(nconf(igr,1))
 125       do i=1,licz(igr)
 126         icon=nconf(igr,i)
 127         boltz=dexp(-(totfree(icon)-e1))
 128         rmsave(igr)=rmsave(igr)+boltz*rmstb(icon)
 129         gdt_ts_ave(igr)=gdt_ts_ave(igr)+boltz*gdt_ts_tb(icon)
 130         gdt_ha_ave(igr)=gdt_ha_ave(igr)+boltz*gdt_ha_tb(icon)
 131         tmscore_ave(igr)=tmscore_ave(igr)+boltz*tmscore_tb(icon)
 132         qpart=qpart+boltz
 133 c        write (iout,'(2i5,10f10.5)') i,icon,boltz,rmstb(icon),
 134 c     &    gdt_ts_tb(icon),gdt_ha_tb(icon),tmscore_tb(icon)
 135       enddo
 136 c      write (iout,*) "qpart",qpart
 137       rmsave(igr)=rmsave(igr)/qpart
 138       gdt_ts_ave(igr)=gdt_ts_ave(igr)/qpart
 139       gdt_ha_ave(igr)=gdt_ha_ave(igr)/qpart
 140       tmscore_ave(igr)=tmscore_ave(igr)/qpart
 141       write (iout,'(a,f5.2,a,3(a,f7.4))') "Cluster averages: RMSD",
 142      & rmsave(igr)," A, ",
 143      & "TMscore",tmscore_ave(igr),
 144      & ", GDT_TS",gdt_ts_ave(igr),", GDT_HA",
 145      & gdt_ha_ave(igr)
 146    19 CONTINUE
 147       WRITE (iout,400)
 148       WRITE (iout,500) (list_conf(I),IASS(I),I=1,NCON)
 149 c      print *,icut,printang(icut)
 150       IF (PRINTANG(ICUT) .and. (lprint_cart .or. lprint_int)) then
 151         emin=totfree_gr(1)
 152 c        print *,'emin',emin,' ngr',ngr
 153         if (lprint_cart) then
 154           cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
 155      &      //"K"//".x"
 156         else
 157           cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
 158      &      //"K"//".int"
 159         endif
 160         do igr=1,ngr
 161           icon=nconf(igr,1)
 162           if (totfree_gr(igr)-emin.le.ecut) then
 163             if (lprint_cart) then
 164               call cartout(igr,icon,totfree(icon)/beta_h(ib),
 165      &          totfree_gr(igr)/beta_h(ib),
 166      &          rmstb(icon),cfname)
 167             else
 168 c              print '(a)','calling briefout'
 169               do i=1,2*nres
 170                 do j=1,3
 171                   c(j,i)=allcart(j,i,icon)
 172                 enddo
 173               enddo
 174               call int_from_cart1(.false.)
 175               call briefout(igr,iscore(icon),totfree(icon)/beta_h(ib),
 176      &          totfree_gr(igr),nss_all(icon),ihpb_all(1,icon),
 177      &          jhpb_all(1,icon),cfname)
 178 c              print '(a)','exit briefout'
 179             endif
 180           endif
 181         enddo
 182         close(igeom)
 183       ENDIF
 184       IF (PRINTPDB(ICUT).gt.0) THEN
 185 c Write out a number of conformations from each family in PDB format and
 186 c create InsightII command file for their displaying in different colors
 187         cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
 188      &    //"K_"//'ave'//exten
 189 c        write (iout,*) "cfname",cfname
 190         OPEN(ipdb,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
 191         write (ipdb,'(a,f8.2)')
 192      &    "REMARK AVERAGE CONFORMATIONS AT TEMPERATURE",temper
 193         close (ipdb)
 194         I=1
 195         ICON=NCONF(1,1)
 196         EMIN=totfree_gr(I)
 197         emin1=totfree(icon)
 198         DO WHILE(I.LE.NGR .AND. totfree_gr(i)-EMIN.LE.ECUT)
 199 c          write (iout,*) "i",i," ngr",ngr,totfree_gr(I),EMIN,ecut
 200           write (NUMM,'(bz,i4.4)') i
 201           ncon_lim=min0(licz(i),printpdb(icut))
 202           cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
 203      &      //"K_"//numm(:ilen(numm))//exten
 204           OPEN(ipdb,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
 205           write (ipdb,'("REMARK CLUSTER",i5," FREE ENERGY",1pe14.5,
 206      &     " AVE RMSD",0pf5.2)')
 207      &     i,totfree_gr(i)/beta_h(ib),rmsave(i)
 208 c Write conformations of the family i to PDB files
 209           ncon_out=1
 210           do while (ncon_out.lt.printpdb(icut) .and.
 211      &     ncon_out.lt.licz(i).and.
 212      &     totfree(nconf(i,ncon_out+1))-EMIN1.LE.ECUT)
 213             ncon_out=ncon_out+1
 214 c            write (iout,*) i,ncon_out,nconf(i,ncon_out),
 215 c     &        totfree(nconf(i,ncon_out)),emin1,ecut
 216           enddo
 217 c          write (iout,*) "ncon_out",ncon_out
 218           call flush(iout)
 219           do j=1,nres
 220             tempfac(1,j)=5.0d0
 221             tempfac(2,j)=5.0d0
 222           enddo
 223           do j=1,ncon_out
 224             icon=nconf(i,j)
 225             do ii=1,2*nres
 226               do k=1,3
 227                 c(k,ii)=allcart(k,ii,icon)
 228               enddo
 229             enddo
 230             nss=nss_all(icon)
 231             write (iout,*) "ICON",icon," nss",nss
 232             do k=1,nss
 233               ihpb(k)=ihpb_all(k,icon)
 234               jhpb(k)=jhpb_all(k,icon)
 235               write (iout,*) ihpb(k),jhpb(k)
 236             enddo
 237             call center
 238             call pdbout(totfree(icon)/beta_h(ib),rmstb(icon),titel)
 239             write (ipdb,'("TER")')
 240           enddo
 241           close(ipdb)
 242 c Average structures and structures closest to average
 243           cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
 244      &    //"K_"//'ave'//exten
 245           OPEN(ipdb,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED',
 246      &     position="APPEND")
 247           call ave_coord(i)
 248           write (ipdb,'(a,i5)') "REMARK CLUSTER",i
 249           call center
 250           nss=0
 251           call pdbout(totfree_gr(i)/beta_h(ib),rmsave(i),titel)
 252           write (ipdb,'("TER")')
 253           if (print_fittest.and.(nsaxs.gt.0 .or. nhpb.gt.0
 254      &     .or.npeak.gt.0)) then
 255             call fittest_coord(i)
 256           else
 257             call closest_coord(i)
 258           endif
 259           if (refstr) then
 260 c            write (iout,*) "Calling rmsnat"
 261             rms_closest(i) = rmsnat(i)
 262 c            write (iout,*) "Cluster",i
 263             call TMscore_sub(rmsd,gdt_ts_closest(i),gdt_ha_closest(i),
 264      &      tmscore_closest(i),cfname,.true.)
 265 c            write (iout,*) "WRTCLUST: nsaxs",nsaxs," i",i
 266           endif
 267           if (nsaxs.gt.0 .and. saxs_mode.eq.0) then
 268             call e_saxs(Esaxs_constr)
 269             Cnorm=0.0d0
 270             do j=1,nsaxs-1
 271               Cnorm=Cnorm+(distsaxs(j+1)-distsaxs(j))*
 272      &             (Pcalc(j+1)+Pcalc(j))/2
 273             enddo
 274             do j=1,nsaxs
 275               Pcalc_all(j,i)=Pcalc(j)/Cnorm
 276             enddo
 277 c            write (iout,*) "Pcalc"
 278 c            write (iout,'(f6.2,f10.5)') (distsaxs(j),Pcalc(j),j=1,nsaxs)
 279             Esaxs_all(i)=Esaxs_constr
 280             write (iout,*) "Esaxs",Esaxs_constr
 281           endif
 282           nviolxlink=0
 283           if (link_start.gt.0) then
 284           do j=link_start,link_end
 285             if (irestr_type(j).eq.10 .or. irestr_type(j).eq. 11) then
 286               dxlink=dist(ihpb(j),jhpb(j))
 287               if (dxlink.le.25.0d0) then
 288               write (iout,'(a,i2,2i5,f8.2)') "XLINK-",
 289      &          irestr_type(j),ihpb(j),jhpb(j),
 290      &          dxlink
 291               else
 292               nviolxlink=nviolxlink+1
 293               write (iout,'(a,i2,2i5,f8.2,2h *)') "XLINK-",
 294      &          irestr_type(j),ihpb(j),jhpb(j),
 295      &          dxlink
 296               endif
 297             endif
 298           enddo
 299           if (nviolxlink.gt.0)
 300      &      write (iout,*) nviolxlink," crosslink violations."
 301 c          write (iout,*) "Family",i," rmsd",rmsd,"gdt_ts",
 302 c     &      gdt_ts_closest(i)," gdt_ha",gdt_ha_closest(i),
 303 c     &      "tmscore",tmscore_closest(i)
 304           endif
 305 c Determine # violated NMR restraints
 306           if (link_end_peak.gt.0) then
 307           nviolpeak=0
 308           write (NUMM,'(bz,i4.4)') i
 309           cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
 310      &    //"K_"//NUMM(:ilen(NUMM))//'.nmr'
 311           open(jrms,file=cfname)
 312           do j=link_start_peak,link_end_peak
 313             viol_nmr=.true.
 314             do ip=ipeak(1,j),ipeak(2,j)
 315               ii=ihpb_peak(ip)
 316               jj=jhpb_peak(ip)
 317               dd=dist(ii,jj)
 318 c              iip=ip-ipeak(1,j)+1
 319 C iii and jjj point to the residues for which the distance is assigned.
 320               if (ii.gt.nres) then
 321                 iii=ii-nres
 322                 jjj=jj-nres
 323                 iiib=1
 324               else
 325                 iii=ii
 326                 jjj=jj
 327                 iiib=0
 328               endif
 329               if (dd.lt.dhpb1_peak(ip)) then
 330                 viol_nmr=.false.
 331 c                write (iout,*) j,iii,jjj,iiib
 332                 write (jrms,'(4i6)') j,iii,jjj,iiib
 333               endif
 334             enddo
 335             if (viol_nmr) then
 336               nviolpeak=nviolpeak+1
 337               list_peak_viol(nviolpeak)=j
 338             endif
 339           enddo
 340           if (nviolpeak.gt.0) then
 341            write (iout,'(a,i5,2h (f8.4,2h%))')
 342      &      "Number of violated NMR restraints:",
 343      &      nviolpeak,100*(nviolpeak+0.)/npeak
 344            write (iout,'(a)')"List of violated restraints:"
 345            write (iout,'(16i5)') (list_peak_viol(j),j=1,nviolpeak)
 346           endif
 347           close(jrms)
 348           endif
 349           if (.not.raw_psipred .and. idihconstr_end.gt.0) then
 350             cfname=prefixp(:ilen(prefixp))//"_T"
 351      &      //ctemper(:ilen(ctemper))
 352      &      //"K_"//NUMM(:ilen(NUMM))//'.angle'
 353           open(jrms,file=cfname)
 354             call int_from_cart1(.false.)
 355             nangviol=0
 356             do j=idihconstr_start,idihconstr_end
 357               itori=idih_constr(j)
 358               phii=phi(itori)
 359               difi=pinorm(phii-phi0(j))
 360               if (difi.gt.drange(j) .or. difi.lt.-drange(j))
 361      &          nangviol=nangviol+1
 362               write (jrms,'(i5,3f10.3)') itori,phii*rad2deg,
 363      &          phi0(j)*rad2deg,rad2deg*drange(j)
 364             enddo
 365             write (iout,'(a,i5)')"Number of angle-restraint violations:"
 366      &           ,nangviol
 367             close(jrms)
 368           endif
 369           call center
 370           call pdbout(totfree_gr(i)/beta_h(ib),rms_closest(i),titel)
 371           write (ipdb,'("TER")')
 372           close (ipdb)
 373           I=I+1
 374           ICON=NCONF(I,1)
 375           emin1=totfree(icon)
 376         ENDDO
 377         ngr_print=i-1
 378         if (nsaxs.gt.0 .and. saxs_mode.eq.0) then
 379           cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
 380      &    //"K_"//'ave'//'.dist'
 381           OPEN(99,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
 382           write (99,'(5h#     ,10f10.5)')
 383      &      (Esaxs_all(i)*wsaxs,i=1,ngr_print)
 384           do j=1,nsaxs
 385             write (99,'(f6.2,10f10.5)') distsaxs(j),
 386      &        (Pcalc_all(j,i),i=1,ngr_print)
 387           enddo
 388           close(99)
 389         endif
 390       ENDIF
 391       IF (printmol2(icut).gt.0) THEN
 392 c Write out a number of conformations from each family in PDB format and
 393 c create InsightII command file for their displaying in different colors
 394         I=1
 395         ICON=NCONF(1,1)
 396         EMIN=ENERGY(ICON)
 397         emin1=emin
 398         DO WHILE(I.LE.NGR .AND. totfree_gr(i)-EMIN.LE.ECUT)
 399           write (NUMM,'(bz,i4.4)') i
 400           cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
 401      &      //"K_"//numm(:ilen(numm))//extmol
 402           OPEN(imol2,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
 403           ncon_out=1
 404           do while (ncon_out.lt.printmol2(icut) .and.
 405      &     ncon_out.lt.licz(i).and.
 406      &     totfree(nconf(i,ncon_out+1))-EMIN1.LE.ECUT)
 407             ncon_out=ncon_out+1
 408           enddo
 409           do j=1,ncon_out
 410             icon=nconf(i,j)
 411             do ii=1,2*nres
 412               do k=1,3
 413                 c(k,ii)=allcart(k,ii,icon)
 414               enddo
 415             enddo
 416             CALL MOL2OUT(totfree(icon)/beta_h(ib),'STRUCTURE'//numm)
 417           enddo
 418           CLOSE(imol2)
 419           I=I+1
 420           ICON=NCONF(I,1)
 421           emin1=totfree(icon)
 422         ENDDO
 423       ENDIF
 424       call WRITE_STATS(ICUT,NCON,IB)
 425   100 FORMAT (//'THERE ARE ',I4,' FAMILIES OF CONFORMATIONS')
 426   200 FORMAT (/'FAMILY ',I4,' WITH TOTAL FREE ENERGY',1pE15.5,
 427      & ' CONTAINS ',I4,' CONFORMATION(S): ')
 428 c 300 FORMAT ( 8(I4,F6.1))
 429   300 FORMAT (5(I4,1pe12.3))
 430   400 FORMAT (//'ASSIGNMENT OF CONSECUTIVE CONFORMATIONS TO FAMILIES:')
 431   500 FORMAT (8(2I4,2X))
 432   600 FORMAT ('REMARK FAMILY',I4,' CONFORMATION',I4,' ENERGY ',E15.6)
 433       RETURN
 434       END
 435 c------------------------------------------------------------------------------
 436       subroutine ave_coord(igr)
 437       implicit none
 438       include 'DIMENSIONS'
 439       include 'sizesclu.dat'
 440       include 'COMMON.CONTROL'
 441       include 'COMMON.CLUSTER'
 442       include 'COMMON.CHAIN'
 443       include 'COMMON.INTERACT'
 444       include 'COMMON.VAR'
 445       include 'COMMON.TEMPFAC'
 446       include 'COMMON.IOUNITS'
 447       logical non_conv
 448       double precision przes(3),obrot(3,3)
 449       double precision xx(3,maxres2),csq(3,maxres2)
 450       double precision eref
 451       double precision rmscalc
 452 c      double precision rmscheck
 453       integer i,ii,j,k,icon,jcon,igr,ipermmin
 454       double precision rms,boltz,qpart,cwork(3,maxres2),cref1(3,maxres2)
 455 c      write (iout,*) "AVE_COORD: igr",igr
 456       jcon=nconf(igr,1)
 457       eref=totfree(jcon)
 458       boltz = dexp(-totfree(jcon)+eref)
 459       qpart=boltz
 460       do i=1,2*nres
 461         do j=1,3
 462           c(j,i)=allcart(j,i,jcon)*boltz
 463           cref1(j,i)=allcart(j,i,jcon)
 464           csq(j,i)=allcart(j,i,jcon)**2*boltz
 465         enddo
 466       enddo
 467       DO K=2,LICZ(IGR)
 468         jcon=nconf(igr,k)
 469 c        write (iout,*) "k",k," jcon",jcon
 470         do i=1,2*nres
 471           do j=1,3
 472             cwork(j,i)=allcart(j,i,jcon)
 473           enddo
 474         enddo
 475         rms=rmscalc(cwork(1,1),cref1(1,1),przes,obrot,ipermmin)
 476 c        write (iout,*) "rms",rms," ipermmin",ipermmin
 477 c        do i=1,3
 478 c          write (iout,'(i3,f10.5,5x,3f10.5)')i,przes(i),
 479 c     &      (obrot(i,j),j=1,3)
 480 c        enddo
 481 c        if (rms.lt.0.0) then
 482 c          print *,'error, rms^2 = ',rms,icon,jcon
 483 c          stop
 484 c        endif
 485 c        if (non_conv) print *,non_conv,icon,jcon
 486         boltz=dexp(-totfree(jcon)+eref)
 487         qpart = qpart + boltz
 488         do i=1,2*nres
 489           do j=1,3
 490             xx(j,i)=allcart(j,i,jcon)
 491           enddo
 492         enddo
 493         call matvec(cwork,obrot,xx,2*nres)
 494         do i=1,2*nres
 495 c          write (iout,'(i5,2(3f10.5,5x))') i,(cwork(j,i),j=1,3),
 496 c     &    (allcart(j,i,jcon),j=1,3)
 497           do j=1,3
 498             cwork(j,i)=cwork(j,i)+przes(j)
 499             c(j,i)=c(j,i)+cwork(j,i)*boltz
 500             csq(j,i)=csq(j,i)+cwork(j,i)**2*boltz
 501           enddo
 502         enddo
 503 c rms check
 504 c        rmscheck=0.0d0
 505 c        do i=nnt,nct
 506 c          do j=1,3
 507 c            rmscheck=rmscheck+(cwork(j,i)-cref1(j,i))**2
 508 c          enddo
 509 c        enddo
 510 c        write (iout,*) "rmscheck",dsqrt(rmscheck/(nct-nnt+1)),rms
 511       ENDDO ! K
 512       do i=1,2*nres
 513         do j=1,3
 514           c(j,i)=c(j,i)/qpart
 515           csq(j,i)=csq(j,i)/qpart-c(j,i)**2
 516         enddo
 517 c        write (iout,'(i5,3f10.5)') i,(csq(j,i),j=1,3)
 518       enddo
 519       do i=nnt,nct
 520         tempfac(1,i)=0.0d0
 521         tempfac(2,i)=0.0d0
 522         do j=1,3
 523           tempfac(1,i)=tempfac(1,i)+csq(j,i)
 524           tempfac(2,i)=tempfac(2,i)+csq(j,i+nres)
 525         enddo
 526         tempfac(1,i)=dsqrt(tempfac(1,i))
 527         tempfac(2,i)=dsqrt(tempfac(2,i))
 528       enddo
 529       return
 530       end
 531 c------------------------------------------------------------------------------
 532       subroutine fittest_coord(igr)
 533       implicit none
 534       include 'DIMENSIONS'
 535       include 'sizesclu.dat'
 536       include 'COMMON.IOUNITS'
 537       include 'COMMON.CONTROL'
 538       include 'COMMON.CLUSTER'
 539       include 'COMMON.CHAIN'
 540       include 'COMMON.INTERACT'
 541       include 'COMMON.SBRIDGE'
 542       include 'COMMON.VAR'
 543       include 'COMMON.FFIELD'
 544       include 'COMMON.TORCNSTR'
 545       include 'COMMON.SAXS'
 546       logical non_conv
 547       double precision przes(3),obrot(3,3)
 548       double precision xx(3,maxres2),yy(3,maxres2)
 549       integer i,ii,j,k,icon,jcon,jconmin,igr
 550       double precision rms,rmsmin,cwork(3,maxres2)
 551       double precision ehpb,Esaxs_constr,edihcnstr
 552       rmsmin=1.0d10
 553       jconmin=nconf(igr,1)
 554       DO K=1,LICZ(IGR)
 555       jcon=nconf(igr,k)
 556       do i=1,2*nres
 557         do j=1,3
 558           c(j,i)=allcart(j,i,jcon)
 559         enddo
 560       enddo
 561       call int_from_cart1(.false.)
 562       esaxs_constr=0
 563       ehpb=0
 564       edihcnstr=0
 565       if (nsaxs.gt.0) call e_saxs(Esaxs_constr)
 566       call edis(ehpb)
 567       if (ndih_constr.gt.0)  call etor_constr(edihcnstr)
 568       rms=wsaxs*esaxs_constr+wstrain*ehpb+edihcnstr
 569 c      write (iout,*) "Esaxs_constr",esaxs_constr," Ehpb",ehpb,
 570 c     & " Edihcnstr",edihcnstr
 571       if (rms.lt.rmsmin) then
 572         jconmin=nconf(igr,k)
 573         rmsmin=rms
 574       endif
 575       ENDDO ! K
 576       write (iout,*) "fittest conformation",jconmin," penalty",rmsmin
 577       do i=1,2*nres
 578         do j=1,3
 579           c(j,i)=allcart(j,i,jconmin)
 580         enddo
 581       enddo
 582       nss=nss_all(jconmin)
 583       do k=1,nss
 584         ihpb(k)=ihpb_all(k,jconmin)
 585         jhpb(k)=jhpb_all(k,jconmin)
 586       enddo
 587       return
 588       end
 589 c------------------------------------------------------------------------------
 590       subroutine closest_coord(igr)
 591       implicit none
 592       include 'DIMENSIONS'
 593       include 'sizesclu.dat'
 594       include 'COMMON.IOUNITS'
 595       include 'COMMON.CONTROL'
 596       include 'COMMON.CLUSTER'
 597       include 'COMMON.CHAIN'
 598       include 'COMMON.INTERACT'
 599       include 'COMMON.VAR'
 600       include 'COMMON.SBRIDGE'
 601       logical non_conv
 602       double precision przes(3),obrot(3,3)
 603       integer i,ii,j,k,icon,jcon,jconmin,igr,ipermmin
 604       double precision rms,rmsmin,cwork(3,maxres2)
 605       double precision xx(3,maxres2),yy(3,maxres2)
 606       double precision rmscalc
 607       rmsmin=1.0d10
 608       jconmin=nconf(igr,1)
 609 c      write (iout,*) "CLOSEST_COORD: Average coords"
 610 c      call cartprint
 611       DO K=1,LICZ(IGR)
 612         jcon=nconf(igr,k)
 613         do i=1,2*nres
 614           do j=1,3
 615             xx(j,i)=c(j,i)
 616             yy(j,i)=allcart(j,i,jcon)
 617           enddo
 618         enddo
 619         rms=rmscalc(xx(1,1),yy(1,1),przes,obrot,ipermmin)
 620 c        write (iout,*) "jcon",jcon," rms",rms," rmsmin",rmsmin
 621         if (non_conv) print *,non_conv,icon,jcon
 622         if (rms.lt.rmsmin) then
 623           rmsmin=rms
 624           jconmin=jcon
 625         endif
 626       ENDDO ! K
 627 c      write (iout,*) "rmsmin",rmsmin," rms",rms
 628       call flush(iout)
 629       do i=1,2*nres
 630         do j=1,3
 631           c(j,i)=allcart(j,i,jconmin)
 632         enddo
 633       enddo
 634       nss=nss_all(jconmin)
 635 c      write (iout,*) "jconmin",jconmin," nss",nss
 636       call flush(iout)
 637       do k=1,nss
 638         ihpb(k)=ihpb_all(k,jconmin)
 639         jhpb(k)=jhpb_all(k,jconmin)
 640 c        write (iout,*) "k",k," ihpb",ihpb(k)," jhpb",jhpb(k)
 641       enddo
 642       call flush(iout)
 643       return
 644       end
 645 c------------------------------------------------------------------------------
 646       subroutine center
 647       implicit none
 648       include 'DIMENSIONS'
 649       include 'sizesclu.dat'
 650       include 'COMMON.IOUNITS'
 651       include 'COMMON.CONTROL'
 652       include 'COMMON.CLUSTER'
 653       include 'COMMON.CHAIN'
 654       include 'COMMON.INTERACT'
 655       include 'COMMON.VAR'
 656       double precision przes(3)
 657       integer i,ii,j,k,icon,jcon,jconmin,igr
 658       przes=0.0d0
 659       do j=1,3
 660         do i=1,nres
 661           przes(j)=przes(j)+c(j,i)
 662         enddo
 663       enddo
 664       do j=1,3
 665         przes(j)=przes(j)/nres
 666       enddo
 667       do i=1,2*nres
 668         do j=1,3
 669           c(j,i)=c(j,i)-przes(j)
 670         enddo
 671       enddo
 672       return
 673       end