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