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