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