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)
  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         write (iout,*) "cfname",cfname
 190         OPEN(ipdb,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
 191         write (ipdb,'(a,f8.2)')
 192      &    "REMAR 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             call center
 231             call pdbout(totfree(icon)/beta_h(ib),rmstb(icon),titel)
 232             write (ipdb,'("TER")')
 233           enddo
 234           close(ipdb)
 235 c Average structures and structures closest to average
 236           cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
 237      &    //"K_"//'ave'//exten
 238           OPEN(ipdb,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED',
 239      &     position="APPEND")
 240           call ave_coord(i)
 241           write (ipdb,'(a,i5)') "REMARK CLUSTER",i
 242           call center
 243           call pdbout(totfree_gr(i)/beta_h(ib),rmsave(i),titel)
 244           write (ipdb,'("TER")')
 245           if (print_fittest.and.(nsaxs.gt.0 .or. nhpb.gt.0
 246      &     .or.npeak.gt.0)) then
 247             call fittest_coord(i)
 248           else
 249             call closest_coord(i)
 250           endif
 251 c            write (iout,*) "Calling rmsnat"
 252           rms_closest(i) = rmsnat(i)
 253
 254           write (iout,*) "Cluster",i
 255           call TMscore_sub(rmsd,gdt_ts_closest(i),gdt_ha_closest(i),
 256      &      tmscore_closest(i),cfname,.true.)
 257 c          write (iout,*) "WRTCLUST: nsaxs",nsaxs," i",i
 258           if (nsaxs.gt.0 .and. saxs_mode.eq.0) then
 259             call e_saxs(Esaxs_constr)
 260             Cnorm=0.0d0
 261             do j=1,nsaxs-1
 262               Cnorm=Cnorm+(distsaxs(j+1)-distsaxs(j))*
 263      &             (Pcalc(j+1)+Pcalc(j))/2
 264             enddo
 265             do j=1,nsaxs
 266               Pcalc_all(j,i)=Pcalc(j)/Cnorm
 267             enddo
 268 c            write (iout,*) "Pcalc"
 269 c            write (iout,'(f6.2,f10.5)') (distsaxs(j),Pcalc(j),j=1,nsaxs)
 270             Esaxs_all(i)=Esaxs_constr
 271             write (iout,*) "Esaxs",Esaxs_constr
 272           endif
 273           nviolxlink=0
 274           if (link_start.gt.0) then
 275           do j=link_start,link_end
 276             if (irestr_type(j).eq.10 .or. irestr_type(j).eq. 11) then
 277               dxlink=dist(ihpb(j),jhpb(j))
 278               if (dxlink.le.25.0d0) then
 279               write (iout,'(a,i2,2i5,f8.2)') "XLINK-",
 280      &          irestr_type(j),ihpb(j),jhpb(j),
 281      &          dxlink
 282               else
 283               nviolxlink=nviolxlink+1
 284               write (iout,'(a,i2,2i5,f8.2,2h *)') "XLINK-",
 285      &          irestr_type(j),ihpb(j),jhpb(j),
 286      &          dxlink
 287               endif
 288             endif
 289           enddo
 290           if (nviolxlink.gt.0)
 291      &      write (iout,*) nviolxlink," crosslink violations."
 292 c          write (iout,*) "Family",i," rmsd",rmsd,"gdt_ts",
 293 c     &      gdt_ts_closest(i)," gdt_ha",gdt_ha_closest(i),
 294 c     &      "tmscore",tmscore_closest(i)
 295           endif
 296 c Determine # violated NMR restraints
 297           if (link_end_peak.gt.0) then
 298           nviolpeak=0
 299           write (NUMM,'(bz,i4.4)') i
 300           cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
 301      &    //"K_"//NUMM(:ilen(NUMM))//'.nmr'
 302           open(jrms,file=cfname)
 303           do j=link_start_peak,link_end_peak
 304             viol_nmr=.true.
 305             do ip=ipeak(1,j),ipeak(2,j)
 306               ii=ihpb_peak(ip)
 307               jj=jhpb_peak(ip)
 308               dd=dist(ii,jj)
 309 c              iip=ip-ipeak(1,j)+1
 310 C iii and jjj point to the residues for which the distance is assigned.
 311               if (ii.gt.nres) then
 312                 iii=ii-nres
 313                 jjj=jj-nres
 314                 iiib=1
 315               else
 316                 iii=ii
 317                 jjj=jj
 318                 iiib=0
 319               endif
 320               if (dd.lt.dhpb1_peak(ip)) then
 321                 viol_nmr=.false.
 322 c                write (iout,*) j,iii,jjj,iiib
 323                 write (jrms,'(4i6)') j,iii,jjj,iiib
 324               endif
 325             enddo
 326             if (viol_nmr) then
 327               nviolpeak=nviolpeak+1
 328               list_peak_viol(nviolpeak)=j
 329             endif
 330           enddo
 331           if (nviolpeak.gt.0) then
 332            write (iout,'(a,i5,2h (f8.4,2h%))')
 333      &      "Number of violated NMR restraints:",
 334      &      nviolpeak,100*(nviolpeak+0.)/npeak
 335            write (iout,'(a)')"List of violated restraints:"
 336            write (iout,'(16i5)') (list_peak_viol(j),j=1,nviolpeak)
 337           endif
 338           close(jrms)
 339           endif
 340           if (.not.raw_psipred .and. idihconstr_end.gt.0) then
 341             cfname=prefixp(:ilen(prefixp))//"_T"
 342      &      //ctemper(:ilen(ctemper))
 343      &      //"K_"//NUMM(:ilen(NUMM))//'.angle'
 344           open(jrms,file=cfname)
 345             call int_from_cart1(.false.)
 346             nangviol=0
 347             do j=idihconstr_start,idihconstr_end
 348               itori=idih_constr(j)
 349               phii=phi(itori)
 350               difi=pinorm(phii-phi0(j))
 351               if (difi.gt.drange(j) .or. difi.lt.-drange(j))
 352      &          nangviol=nangviol+1
 353               write (jrms,'(i5,3f10.3)') itori,phii*rad2deg,
 354      &          phi0(j)*rad2deg,rad2deg*drange(j)
 355             enddo
 356             write (iout,'(a,i5)')"Number of angle-restraint violations:"
 357      &           ,nangviol
 358             close(jrms)
 359           endif
 360           call center
 361           call pdbout(totfree_gr(i)/beta_h(ib),rms_closest(i),titel)
 362           write (ipdb,'("TER")')
 363           close (ipdb)
 364           I=I+1
 365           ICON=NCONF(I,1)
 366           emin1=totfree(icon)
 367         ENDDO
 368         ngr_print=i-1
 369         if (nsaxs.gt.0 .and. saxs_mode.eq.0) then
 370           cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
 371      &    //"K_"//'ave'//'.dist'
 372           OPEN(99,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
 373           write (99,'(5h#     ,10f10.5)')
 374      &      (Esaxs_all(i)*wsaxs,i=1,ngr_print)
 375           do j=1,nsaxs
 376             write (99,'(f6.2,10f10.5)') distsaxs(j),
 377      &        (Pcalc_all(j,i),i=1,ngr_print)
 378           enddo
 379           close(99)
 380         endif
 381       ENDIF
 382       IF (printmol2(icut).gt.0) THEN
 383 c Write out a number of conformations from each family in PDB format and
 384 c create InsightII command file for their displaying in different colors
 385         I=1
 386         ICON=NCONF(1,1)
 387         EMIN=ENERGY(ICON)
 388         emin1=emin
 389         DO WHILE(I.LE.NGR .AND. totfree_gr(i)-EMIN.LE.ECUT)
 390           write (NUMM,'(bz,i4.4)') i
 391           cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
 392      &      //"K_"//numm(:ilen(numm))//extmol
 393           OPEN(imol2,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
 394           ncon_out=1
 395           do while (ncon_out.lt.printmol2(icut) .and.
 396      &     ncon_out.lt.licz(i).and.
 397      &     totfree(nconf(i,ncon_out+1))-EMIN1.LE.ECUT)
 398             ncon_out=ncon_out+1
 399           enddo
 400           do j=1,ncon_out
 401             icon=nconf(i,j)
 402             do ii=1,2*nres
 403               do k=1,3
 404                 c(k,ii)=allcart(k,ii,icon)
 405               enddo
 406             enddo
 407             CALL MOL2OUT(totfree(icon)/beta_h(ib),'STRUCTURE'//numm)
 408           enddo
 409           CLOSE(imol2)
 410           I=I+1
 411           ICON=NCONF(I,1)
 412           emin1=totfree(icon)
 413         ENDDO
 414       ENDIF
 415       call WRITE_STATS(ICUT,NCON,IB)
 416   100 FORMAT (//'THERE ARE ',I4,' FAMILIES OF CONFORMATIONS')
 417   200 FORMAT (/'FAMILY ',I4,' WITH TOTAL FREE ENERGY',1pE15.5,
 418      & ' CONTAINS ',I4,' CONFORMATION(S): ')
 419 c 300 FORMAT ( 8(I4,F6.1))
 420   300 FORMAT (5(I4,1pe12.3))
 421   400 FORMAT (//'ASSIGNMENT OF CONSECUTIVE CONFORMATIONS TO FAMILIES:')
 422   500 FORMAT (8(2I4,2X))
 423   600 FORMAT ('REMARK FAMILY',I4,' CONFORMATION',I4,' ENERGY ',E15.6)
 424       RETURN
 425       END
 426 c------------------------------------------------------------------------------
 427       subroutine ave_coord(igr)
 428       implicit none
 429       include 'DIMENSIONS'
 430       include 'sizesclu.dat'
 431       include 'COMMON.CONTROL'
 432       include 'COMMON.CLUSTER'
 433       include 'COMMON.CHAIN'
 434       include 'COMMON.INTERACT'
 435       include 'COMMON.VAR'
 436       include 'COMMON.TEMPFAC'
 437       include 'COMMON.IOUNITS'
 438       logical non_conv
 439       double precision przes(3),obrot(3,3)
 440       double precision xx(3,maxres2),csq(3,maxres2)
 441       double precision eref
 442       double precision rmscalc
 443 c      double precision rmscheck
 444       integer i,ii,j,k,icon,jcon,igr,ipermmin
 445       double precision rms,boltz,qpart,cwork(3,maxres2),cref1(3,maxres2)
 446 c      write (iout,*) "AVE_COORD: igr",igr
 447       jcon=nconf(igr,1)
 448       eref=totfree(jcon)
 449       boltz = dexp(-totfree(jcon)+eref)
 450       qpart=boltz
 451       do i=1,2*nres
 452         do j=1,3
 453           c(j,i)=allcart(j,i,jcon)*boltz
 454           cref1(j,i)=allcart(j,i,jcon)
 455           csq(j,i)=allcart(j,i,jcon)**2*boltz
 456         enddo
 457       enddo
 458       DO K=2,LICZ(IGR)
 459         jcon=nconf(igr,k)
 460 c        write (iout,*) "k",k," jcon",jcon
 461         do i=1,2*nres
 462           do j=1,3
 463             cwork(j,i)=allcart(j,i,jcon)
 464           enddo
 465         enddo
 466         rms=rmscalc(cwork(1,1),cref1(1,1),przes,obrot,ipermmin)
 467 c        write (iout,*) "rms",rms," ipermmin",ipermmin
 468 c        do i=1,3
 469 c          write (iout,'(i3,f10.5,5x,3f10.5)')i,przes(i),
 470 c     &      (obrot(i,j),j=1,3)
 471 c        enddo
 472 c        if (rms.lt.0.0) then
 473 c          print *,'error, rms^2 = ',rms,icon,jcon
 474 c          stop
 475 c        endif
 476 c        if (non_conv) print *,non_conv,icon,jcon
 477         boltz=dexp(-totfree(jcon)+eref)
 478         qpart = qpart + boltz
 479         do i=1,2*nres
 480           do j=1,3
 481             xx(j,i)=allcart(j,i,jcon)
 482           enddo
 483         enddo
 484         call matvec(cwork,obrot,xx,2*nres)
 485         do i=1,2*nres
 486 c          write (iout,'(i5,2(3f10.5,5x))') i,(cwork(j,i),j=1,3),
 487 c     &    (allcart(j,i,jcon),j=1,3)
 488           do j=1,3
 489             cwork(j,i)=cwork(j,i)+przes(j)
 490             c(j,i)=c(j,i)+cwork(j,i)*boltz
 491             csq(j,i)=csq(j,i)+cwork(j,i)**2*boltz
 492           enddo
 493         enddo
 494 c rms check
 495 c        rmscheck=0.0d0
 496 c        do i=nnt,nct
 497 c          do j=1,3
 498 c            rmscheck=rmscheck+(cwork(j,i)-cref1(j,i))**2
 499 c          enddo
 500 c        enddo
 501 c        write (iout,*) "rmscheck",dsqrt(rmscheck/(nct-nnt+1)),rms
 502       ENDDO ! K
 503       do i=1,2*nres
 504         do j=1,3
 505           c(j,i)=c(j,i)/qpart
 506           csq(j,i)=csq(j,i)/qpart-c(j,i)**2
 507         enddo
 508 c        write (iout,'(i5,3f10.5)') i,(csq(j,i),j=1,3)
 509       enddo
 510       do i=nnt,nct
 511         tempfac(1,i)=0.0d0
 512         tempfac(2,i)=0.0d0
 513         do j=1,3
 514           tempfac(1,i)=tempfac(1,i)+csq(j,i)
 515           tempfac(2,i)=tempfac(2,i)+csq(j,i+nres)
 516         enddo
 517         tempfac(1,i)=dsqrt(tempfac(1,i))
 518         tempfac(2,i)=dsqrt(tempfac(2,i))
 519       enddo
 520       return
 521       end
 522 c------------------------------------------------------------------------------
 523       subroutine fittest_coord(igr)
 524       implicit none
 525       include 'DIMENSIONS'
 526       include 'sizesclu.dat'
 527       include 'COMMON.IOUNITS'
 528       include 'COMMON.CONTROL'
 529       include 'COMMON.CLUSTER'
 530       include 'COMMON.CHAIN'
 531       include 'COMMON.INTERACT'
 532       include 'COMMON.VAR'
 533       include 'COMMON.FFIELD'
 534       include 'COMMON.TORCNSTR'
 535       include 'COMMON.SAXS'
 536       logical non_conv
 537       double precision przes(3),obrot(3,3)
 538       double precision xx(3,maxres2),yy(3,maxres2)
 539       integer i,ii,j,k,icon,jcon,jconmin,igr
 540       double precision rms,rmsmin,cwork(3,maxres2)
 541       double precision ehpb,Esaxs_constr,edihcnstr
 542       rmsmin=1.0d10
 543       jconmin=nconf(igr,1)
 544       DO K=1,LICZ(IGR)
 545       jcon=nconf(igr,k)
 546       do i=1,2*nres
 547         do j=1,3
 548           c(j,i)=allcart(j,i,jcon)
 549         enddo
 550       enddo
 551       call int_from_cart1(.false.)
 552       esaxs_constr=0
 553       ehpb=0
 554       edihcnstr=0
 555       if (nsaxs.gt.0) call e_saxs(Esaxs_constr)
 556       call edis(ehpb)
 557       if (ndih_constr.gt.0)  call etor_constr(edihcnstr)
 558       rms=wsaxs*esaxs_constr+wstrain*ehpb+edihcnstr
 559 c      write (iout,*) "Esaxs_constr",esaxs_constr," Ehpb",ehpb,
 560 c     & " Edihcnstr",edihcnstr
 561       if (rms.lt.rmsmin) then
 562         jconmin=nconf(igr,k)
 563         rmsmin=rms
 564       endif
 565       ENDDO ! K
 566       write (iout,*) "fittest conformation",jconmin," penalty",rmsmin
 567       do i=1,2*nres
 568         do j=1,3
 569           c(j,i)=allcart(j,i,jconmin)
 570         enddo
 571       enddo
 572       return
 573       end
 574 c------------------------------------------------------------------------------
 575       subroutine closest_coord(igr)
 576       implicit none
 577       include 'DIMENSIONS'
 578       include 'sizesclu.dat'
 579       include 'COMMON.IOUNITS'
 580       include 'COMMON.CONTROL'
 581       include 'COMMON.CLUSTER'
 582       include 'COMMON.CHAIN'
 583       include 'COMMON.INTERACT'
 584       include 'COMMON.VAR'
 585       logical non_conv
 586       double precision przes(3),obrot(3,3)
 587       integer i,ii,j,k,icon,jcon,jconmin,igr,ipermmin
 588       double precision rms,rmsmin,cwork(3,maxres2)
 589       double precision xx(3,maxres2),yy(3,maxres2)
 590       double precision rmscalc
 591       rmsmin=1.0d10
 592       jconmin=nconf(igr,1)
 593       DO K=1,LICZ(IGR)
 594         jcon=nconf(igr,k)
 595         do i=1,2*nres
 596           do j=1,3
 597             xx(j,i)=c(j,i)
 598             yy(j,i)=allcart(j,i,jcon)
 599           enddo
 600         enddo
 601         rms=rmscalc(xx(1,1),yy(1,1),przes,obrot,ipermmin)
 602 c        write (iout,*) "jcon",jcon," rms",rms," rmsmin",rmsmin
 603         if (non_conv) print *,non_conv,icon,jcon
 604         if (rms.lt.rmsmin) then
 605           rmsmin=rms
 606           jconmin=jcon
 607         endif
 608       ENDDO ! K
 609 c      write (iout,*) "rmsmin",rmsmin," rms",rms
 610 c      call flush(iout)
 611       do i=1,2*nres
 612         do j=1,3
 613           c(j,i)=allcart(j,i,jconmin)
 614         enddo
 615       enddo
 616       return
 617       end
 618 c------------------------------------------------------------------------------
 619       subroutine center
 620       implicit none
 621       include 'DIMENSIONS'
 622       include 'sizesclu.dat'
 623       include 'COMMON.IOUNITS'
 624       include 'COMMON.CONTROL'
 625       include 'COMMON.CLUSTER'
 626       include 'COMMON.CHAIN'
 627       include 'COMMON.INTERACT'
 628       include 'COMMON.VAR'
 629       double precision przes(3)
 630       integer i,ii,j,k,icon,jcon,jconmin,igr
 631       przes=0.0d0
 632       do j=1,3
 633         do i=1,nres
 634           przes(j)=przes(j)+c(j,i)
 635         enddo
 636       enddo
 637       do j=1,3
 638         przes(j)=przes(j)/nres
 639       enddo
 640       do i=1,2*nres
 641         do j=1,3
 642           c(j,i)=c(j,i)-przes(j)
 643         enddo
 644       enddo
 645       return
 646       end