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