      Program LOCSYM                                                    LOCS001
c --------------------------------------------------------------------- LOCS002
c     VARIABLES:                                                        LOCS003
c       cel    =  cell constants                                        LOCS004
c       t      =  coordinate transformation matrix                      LOCS005
c       x,y,z  =  cartesion coordinates                                 LOCS006
c       charge =  charge of ligand atom                                 LOCS007
c       dis    =  distance from central atom                            LOCS008
c       el     =  element name                                          LOCS009
c       rmax   =  largest permitted radius                              LOCS010
c                                                                       LOCS011
c       Compatible with SEXIE 4.0 - B.Rupp (1994)                       LOCS012
c       Compatible with SEXIE 5.0 - B.Rupp (1995)                       LOCS013
c                                                                       LOCS014
c     CONTROL VARIABLES                                                 LOCS015
c       iatoms = total number of atoms found in inputfile               LOCS016
c       ibeg   = first atom to be read                                  LOCS017
c       iend   = last atom to be read                                   LOCS018
c       nshells= highest shell to be read                               LOCS019
c       isumat = numer of atoms actually used                           LOCS020
c --------------------------------------------------------------------- LOCS021
      parameter (ims=100, maxats=48)                                    LOCS022
c --- maximum number of atoms in shell is maxats = 48   (SEXIE limit)   LOCS023
c     maximum number of shells to be read in ims = 100  (SEXIE limit)   LOCS024
c     NOTE : SEXIE limits need to be changed in all subroutines !       LOCS025
c     (limits should be higher for low symmetry problems)               LOCS026
                                                                        LOCS027
      character infile*25,outnam*10,el*2,revdat*8,icomp*4,creatr*50     LOCS028
      character title*60,spcgrp*60,ide*4,label*7,answ,inname*10,ln*70   LOCS029
      character lv,de,an                                                LOCS030
      integer sltodo(10),ishell(ims*maxats)                             LOCS031
      logical exists, isnew                                             LOCS032
      real coords(6,ims*maxats)                                         LOCS033
      real m(3,3,ims*maxats),a(3,3), d(3), v(3,3)                       LOCS034
      dimension x(ims*maxats), y(ims*maxats), z(ims*maxats)             LOCS035
      dimension dis(ims*maxats),ide(ims*maxats),foc(ims*maxats)         LOCS036
      dimension charge(ims*maxats),t(3,3),cel(6), el(ims*maxats)        LOCS037
      inname='charge.inp'                                               LOCS038
      outnam='charge.out'                                               LOCS039
                                                                        LOCS040
c --- revision comments and date here ----                              LOCS041
c     version  5.0  by B.Rupp                                           LOCS042
      revdat='08/21/95'                                                 LOCS043
                                                                        LOCS044
c --- look for SEXIE configuration file and read the computer type --   LOCS045
      inquire (file='SEXIE.CFG',exist=exists)                           LOCS046
      if (exists) then                                                  LOCS047
         open (4,file='SEXIE.CFG',status='old')                         LOCS048
         do 0001 i=1,7                                                  LOCS049
 0001       read (4,'(8x,a4)') icomp                                    LOCS050
         read(4,'(8x,a4)')   icomp                                      LOCS051
         call upstrg (icomp,4)                                          LOCS052
         close(4,status='keep')                                         LOCS053
      else                                                              LOCS054
         write (*,'(/a)') '  Cannot find configuration file SEXIE.CFG'  LOCS055
         write (*,*) ' Resume with 7-bit TTY terminal character set'    LOCS056
         icomp  = '7BIT'                                                LOCS057
      end if                                                            LOCS058
                                                                        LOCS059
c --- initialize special characters ---                                 LOCS060
      if (icomp.eq.'IBM ') then                                         LOCS061
         lv=char(179)                                                   LOCS062
         an=char(143)                                                   LOCS063
         de=char(248)                                                   LOCS064
      elseif (icomp.eq.'VAX ') then                                     LOCS065
         lv='|'                                                         LOCS066
         an=char(197)                                                   LOCS067
         de=char(176)                                                   LOCS068
      else                                                              LOCS069
         lv='|'                                                         LOCS070
         an='A'                                                         LOCS071
         de=' '                                                         LOCS072
      end if                                                            LOCS073
                                                                        LOCS074
c --- zero out x,y,z and distance ----                                  LOCS075
      do 0002 i=1,ims*maxats                                            LOCS076
          dis(i)=0.                                                     LOCS077
            x(i)=0.                                                     LOCS078
            y(i)=0.                                                     LOCS079
            z(i)=0.                                                     LOCS080
       charge(i)=0.                                                     LOCS081
          foc(i)=0.                                                     LOCS082
 0002 continue                                                          LOCS083
                                                                        LOCS084
c --- read file charge.inp and open charge.out ---------------          LOCS085
      open (5,file=inname,status='old')                                 LOCS086
      open (7,file=outnam,status='unknown')                             LOCS087
c --- junk line --                                                      LOCS088
      read (5,'(1x,a7,a60)') label,title                                LOCS089
c --- creator line --                                                   LOCS090
      read (5,'(2x,a50)') creatr                                        LOCS091
c --- junk line --                                                      LOCS092
      read (5,'(1x,a7,a60)') label,title                                LOCS093
c --- TITLE line --                                                     LOCS094
      read (5,'(1x,a7,a60)') label,title                                LOCS095
c --- SPCGRP line --                                                    LOCS096
      read (5,'(1x,a7,a60)') label,spcgrp                               LOCS097
                                                                        LOCS098
c --- CELL line --                                                      LOCS099
      read (5,'(1x,a7,3f8.4,3f7.2)') label,(cel(i), i=1,6)              LOCS100
                                                                        LOCS101
c --- CORMAT line (coordinate transformation matrix)                    LOCS102
      do 0004 j=1,3                                                     LOCS103
         read (5,'(1x,a7,3f12.6)')label,  (t(j,i),i=1,3)                LOCS104
 0004 continue                                                          LOCS105
                                                                        LOCS106
c --- SHELLS line                                                       LOCS107
      read (5,'(1x,a7,10i3)') label,(sltodo(i),i=1,10)                  LOCS108
c      write(*,'(1x,a7,10i3)') label,(sltodo(i),i=1,10)                 LOCS109
      if (sltodo(1).lt.1) then                                          LOCS110
         isnew=.true.                                                   LOCS111
      else                                                              LOCS112
         isnew=.false.                                                  LOCS113
      end if                                                            LOCS114
c --- RMAX                                                              LOCS115
      read (5,'(1x,a7,f8.4)') label,rmax                                LOCS116
c --- WEIGHT                                                            LOCS117
      read (5,'(1x,a7,i8)') label,iweigt                                LOCS118
c --- USECHG                                                            LOCS119
      read (5,'(1x,a7,i8)') label,iusech                                LOCS120
c --- CENTRE                                                            LOCS121
      read (5,'(1x,a7)') label                                          LOCS122
      i=1                                                               LOCS123
      read (5,'(i4,1x,a4,a4,6f8.4)') ishell(i),el(i),                   LOCS124
     &      ide(i),x(i),y(i),z(i),dis(i),charge(i),foc(i)               LOCS125
c --- ATOMS                                                             LOCS126
      read (5,'(1x,a7)') label                                          LOCS127
      read (5,'(1x,a7)') label                                          LOCS128
c --- loop over all atoms listed in all shells ---                      LOCS129
      do 0005 i=2,maxats*ims                                            LOCS130
         read (5,'(i4,1x,a4,a4,6f8.4)',err=0106,end=0106) ishell(i),    LOCS131
     &   el(i),ide(i),x(i),y(i),z(i),dis(i),charge(i),foc(i)            LOCS132
 0005 continue                                                          LOCS133
 0106 iatoms=i-1                                                        LOCS134
      close(5)                                                          LOCS135
                                                                        LOCS136
c --- header ------                                                     LOCS137
      call chrhed (revdat,icomp,isnew,creatr,title,spcgrp,cel,el(1),    LOCS138
     &      ide(1),x(1),y(1),z(1),lv)                                   LOCS139
                                                                        LOCS140
c --- notify if new run --                                              LOCS141
      if (isnew) then                                                   LOCS142
         write (*,'(/a,a/)')'  This appears to be a new run.' ,         LOCS143
     &                     ' You need to answer a few questions.'       LOCS144
c --- select shell or radius ---                                        LOCS145
 9000 write(*,'(a/,a$)')      '  You can either select a maximum radius LOCS146
     &or specific shells to calculate.','  Radius (R) or Shells (S) : ' LOCS147
      call inkey(answ)                                                  LOCS148
c --- get the radius ---                                                LOCS149
      if (answ.eq.'R') then                                             LOCS150
 9001    write (*,'(a$)')                                               LOCS151
     &   '  Maximum distance of atoms (RMAX) to include : '             LOCS152
         read (*,*,err=9001) rmax                                       LOCS153
         write(*,'(a,f6.3,a)') '  Max radius = ',rmax, 'A'              LOCS154
c ---    see how many atoms and shells lie within rmax ---              LOCS155
         ibeg=2                                                         LOCS156
         do 0006 i=ibeg,iatoms-1                                        LOCS157
            if (dis(i).le.rmax) then                                    LOCS158
               nshells=ishell(i)                                        LOCS159
            else                                                        LOCS160
               goto 0007                                                LOCS161
            end if                                                      LOCS162
 0006    continue                                                       LOCS163
 0007    continue                                                       LOCS164
         rmax=dis(i-1)                                                  LOCS165
         iend=i-1                                                       LOCS166
c        print*,'ibeg,iend=',ibeg,iend                                  LOCS167
                                                                        LOCS168
c ---    tell how many shells ---                                       LOCS169
         write(*,'(/a$)')'  Shells included : '                         LOCS170
         do 0008 i=1,nshells                                            LOCS171
            write(*,'(i3$)') i                                          LOCS172
            sltodo(i)=i                                                 LOCS173
 0008    continue                                                       LOCS174
         write(*,*)                                                     LOCS175
c ---    ask if list desired ---                                        LOCS176
         write(*,'(a$)')                                                LOCS177
     &      '  Do you want to see a list of the selected atoms (Y/N) : 'LOCS178
         call inkey(answ)                                               LOCS179
         if (answ.eq.'Y') then                                          LOCS180
            write(*,'(/a/)') '  List of atoms follows'                  LOCS181
            do 0011 i=ibeg,iend                                         LOCS182
               write(*,'(i4,1x,a4,a4,6f8.4)') ishell(i),el(i),          LOCS183
     &              ide(i),x(i),y(i),z(i),dis(i),charge(i),foc(i)       LOCS184
 0011       continue                                                    LOCS185
         end if                                                         LOCS186
c --- get the shell numbers ---                                         LOCS187
      else if (answ.eq.'S') then                                        LOCS188
         i=1                                                            LOCS189
 9002    write(*,'(a$)')'  Shell Number to calculate (#, E to end) : '  LOCS190
         read(*,*,err=9003) sltodo(i)                                   LOCS191
         if(sltodo(i).gt.0) i=i+1                                       LOCS192
         goto 9002                                                      LOCS193
 9003    continue                                                       LOCS194
         nshells=i-1                                                    LOCS195
c         write(*,'(10i5)') (sltodo(i),i=1,10)                          LOCS196
c ---    sort the shells properly and strip zeros ---                   LOCS197
         call nulsrt(10,sltodo)                                         LOCS198
c ---    tell which shells ---                                          LOCS199
         write(*,'(/a$)')'  Shells included : '                         LOCS200
         do 0012 j=1,nshells                                            LOCS201
            write(*,'(i3$)') sltodo(j)                                  LOCS202
 0012    continue                                                       LOCS203
         write(*,*)                                                     LOCS204
c ---    determine ibeg and iend ---                                    LOCS205
         ibeg=0                                                         LOCS206
         do 1000 i=2,1+ims*maxats                                       LOCS207
            if (ishell(i).eq.sltodo(1)) then                            LOCS208
               ibeg=i                                                   LOCS209
               goto 1001                                                LOCS210
            end if                                                      LOCS211
 1000    continue                                                       LOCS212
 1001    continue                                                       LOCS213
         iend=0                                                         LOCS214
         do 1002 i=2,1+ims*maxats                                       LOCS215
            if (ishell(i).gt.sltodo(nshells)) then                      LOCS216
               iend=i-1                                                 LOCS217
               goto 1003                                                LOCS218
            end if                                                      LOCS219
 1002    continue                                                       LOCS220
 1003    continue                                                       LOCS221
c        print*,'ibeg,iend=',ibeg,iend                                  LOCS222
c ---    ask if list desired ---                                        LOCS223
         write(*,'(/a$)')                                               LOCS224
     &      '  Do you want to see a list of the selected atoms (Y/N) : 'LOCS225
         call inkey(answ)                                               LOCS226
         if (answ.eq.'Y') then                                          LOCS227
            write(*,'(/a/)') '  List of atoms follows'                  LOCS228
            rmax=0.01                                                   LOCS229
            do 0013 i=ibeg, iend                                        LOCS230
               do 0014 j=1,nshells                                      LOCS231
               if (ishell(i).eq.sltodo(j)) then                         LOCS232
                  write(*,'(i4,1x,a4,a4,6f8.4)') ishell(i),el(i),       LOCS233
     &            ide(i),x(i),y(i),z(i),dis(i),charge(i),foc(i)         LOCS234
                  if (dis(i).gt.rmax) rmax=dis(i)                       LOCS235
               end if                                                   LOCS236
 0014          continue                                                 LOCS237
 0013       continue                                                    LOCS238
         end if                                                         LOCS239
      else                                                              LOCS240
         goto 9000                                                      LOCS241
      end if                                                            LOCS242
                                                                        LOCS243
c --- now, for weight ---                                               LOCS244
 9014 write(*,'(/a)')'  You have to select a weight W for 1/R**W.'      LOCS245
      write(*,'(a)') '  0 means no weighting in r (geometry only).'     LOCS246
      write(*,'(a$)')'  Your choice of weight (0 to 9) : '              LOCS247
      read(*,*,err=9014) iweigt                                         LOCS248
                                                                        LOCS249
c --- shall we include the charge ---                                   LOCS250
 9914 write(*,'(/a$)')'  How do you want to use the individual charges 'LOCS251
      write(*,'(a)')  '  of the ligands (charge**n) ?'                  LOCS252
      write(*,'(a)')  '  No charge          (0)'                        LOCS253
      write(*,'(a)')  '  Linear charge      (1)'                        LOCS254
      write(*,'(a$)') '  Squared charge     (2) : '                     LOCS255
                                                                        LOCS256
      read(*,'(i1)',err=9914) iusech                                    LOCS257
      if ((iusech.gt.2).or.(iusech.lt.0)) goto 9914                     LOCS258
                                                                        LOCS259
c --- verify selection of parameters ---                                LOCS260
      write(*,'(/a/)') '  Your selection of parameters : '              LOCS261
      write(*,'(1x,a7,10i3)')' SHELLS ',(sltodo(i),i=1,10)              LOCS262
      write(*,'(1x,a7,f8.4,a)')' RMAX   ',rmax,'   Max. radial distance'LOCS263
      write(*,'(1x,a7,i8,a)')' WEIGHT ',iweigt,'   Weight (1/r**w)'     LOCS264
      write(*,'(1x,a7,i8,a)')' USECHG ',iusech,'   Use charge '         LOCS265
                                                                        LOCS266
c --- get a file name --                                                LOCS267
      write(*,'(a$)') '  Create input file (Y) or redo (N) : '          LOCS268
      call inkey(answ)                                                  LOCS269
      if (answ.eq.'Y') then                                             LOCS270
 9200  write(*,'(a$)')'  Filename (max 8 chars, no extension,please) : 'LOCS271
         read(*,'(a25)',err=9200) infile                                LOCS272
         call filnam(infile,'.CHG',0)                                   LOCS273
c ---    write new input file ---                                       LOCS274
         open(8,file=infile,status='unknown',err=9200)                  LOCS275
         open (5,file=inname,status='old')                              LOCS276
         do 0015 i=1,9                                                  LOCS277
            read(5,'(a)')ln                                             LOCS278
            write(8,'(a)')ln                                            LOCS279
 0015    continue                                                       LOCS280
      write(8,'(1x,a7,10i3)')'SHELLS ',(sltodo(i),i=1,10)               LOCS281
      write(8,'(1x,a7,f8.4,a)')'RMAX   ',rmax,'   Max. radial distance' LOCS282
      write(8,'(1x,a7,i8,a)')'WEIGHT ',iweigt,'   Weight (1/r**w)'      LOCS283
      write(8,'(1x,a7,i8,a)')'USECHG ',iusech,'   Use charge '          LOCS284
         do 0016 i=1,4                                                  LOCS285
            read(5,'(a)',end=9300,err=9300)ln                           LOCS286
 0016    continue                                                       LOCS287
         do 0017 i=1,4                                                  LOCS288
            read(5,'(a)',end=9300,err=9300)ln                           LOCS289
            write(8,'(a)')ln                                            LOCS290
 0017    continue                                                       LOCS291
c ---    print out and sum up the atoms ---                             LOCS292
         isumat=0                                                       LOCS293
         do 0019 i=ibeg, iend                                           LOCS294
            do 0020 j=1,nshells                                         LOCS295
            if (ishell(i).eq.sltodo(j)) then                            LOCS296
               write(8,'(i4,1x,a4,a4,6f8.4)') ishell(i),el(i),          LOCS297
     &         ide(i),x(i),y(i),z(i),dis(i),charge(i),foc(i)            LOCS298
               isumat=isumat+1                                          LOCS299
            end if                                                      LOCS300
 0020         continue                                                  LOCS301
 0019    continue                                                       LOCS302
c ---     this writes all positions ---                                 LOCS303
c         do 0018 i=1,maxats*ims                                        LOCS304
c            read(5,'(a)',end=9300,err=9300)ln                          LOCS305
c            write(8,'(a)')ln                                           LOCS306
c 0018    continue                                                      LOCS307
 9300    continue                                                       LOCS308
         write(*,'(a,a/)') '  Created input file ',infile               LOCS309
         close(8)                                                       LOCS310
         close(5)                                                       LOCS311
      else                                                              LOCS312
         goto 9000                                                      LOCS313
      end if                                                            LOCS314
      else                                                              LOCS315
c ---    batch input sort the shells properly and strip zeros ---       LOCS316
         call nulsrt(10,sltodo)                                         LOCS317
c ---    tell which shells ---                                          LOCS318
         nshells=0                                                      LOCS319
         do 9204 i=1,10                                                 LOCS320
            if(sltodo(i).gt.0) nshells=nshells+1                        LOCS321
 9204    continue                                                       LOCS322
c ---    determine ibeg and iend from batch input                       LOCS323
         ibeg=0                                                         LOCS324
         do 1200 i=2,1+ims*maxats                                       LOCS325
            if (ishell(i).eq.sltodo(1)) then                            LOCS326
               ibeg=i                                                   LOCS327
               goto 1201                                                LOCS328
            end if                                                      LOCS329
 1200    continue                                                       LOCS330
 1201    continue                                                       LOCS331
         iend=0                                                         LOCS332
         do 1202 i=2,1+ims*maxats                                       LOCS333
            if (ishell(i).gt.sltodo(nshells)) then                      LOCS334
               iend=i-1                                                 LOCS335
               goto 1203                                                LOCS336
            end if                                                      LOCS337
 1202    continue                                                       LOCS338
 1203    continue                                                       LOCS339
         isumat=0                                                       LOCS340
         do 1219 i=ibeg, iend                                           LOCS341
            do 1220 j=1,nshells                                         LOCS342
               if (ishell(i).eq.sltodo(j)) isumat=isumat+1              LOCS343
 1220       continue                                                    LOCS344
 1219    continue                                                       LOCS345
      end if                                                            LOCS346
c --- end if isnew --- BATCH execution follows --- output first ---     LOCS347
                                                                        LOCS348
      call parout                                                       LOCS349
     & (icomp,sltodo,rmax,iweigt,iusech,an,lv,isumat)                   LOCS350
c                    ishell,el,ide,x,y,x,dis,charge                     LOCS351
      do 0029 ipr=6,7                                                   LOCS352
         do 0031 i=ibeg,iend                                            LOCS353
            do 0030 j=1,nshells                                         LOCS354
            if (ishell(i).eq.sltodo(j)) then                            LOCS355
              write(ipr,'(1x,a,i4,3x,a4,5x,a4,5(2x,f8.4),f6.3,1x,a)')lv,LOCS356
     &  ishell(i),el(i),ide(i),x(i),y(i),z(i),dis(i),charge(i),foc(i),lvLOCS357
            end if                                                      LOCS358
 0030       continue                                                    LOCS359
 0031    continue                                                       LOCS360
         call frames (ipr,3,icomp)                                      LOCS361
 0029 continue                                                          LOCS362
                                                                        LOCS363
c --- set upper limit of shells to rmax ---                             LOCS364
c      ibeg=2                                                           LOCS365
c      do 0035 i=ibeg,iatoms-1                                          LOCS366
c         if (dis(i).le.rmax) then                                      LOCS367
c            nshells=ishell(i)                                          LOCS368
c         else                                                          LOCS369
c            goto 0036                                                  LOCS370
c         end if                                                        LOCS371
c 0035 continue                                                         LOCS372
c 0036 continue                                                         LOCS373
                                                                        LOCS374
c --- rearrange input x,y,x, distance into coords(i,j) ---              LOCS375
      isumat=1                                                          LOCS376
      do 0041 i=ibeg,iend                                               LOCS377
         do 0042 j=1,nshells                                            LOCS378
            if (ishell(i).eq.sltodo(j)) then                            LOCS379
               coords(1,isumat)=x(i)                                    LOCS380
               coords(2,isumat)=y(i)                                    LOCS381
               coords(3,isumat)=z(i)                                    LOCS382
               coords(4,isumat)=dis(i)                                  LOCS383
               coords(5,isumat)=charge(i)                               LOCS384
               coords(6,isumat)=foc(i)                                  LOCS385
               isumat=isumat+1                                          LOCS386
            end if                                                      LOCS387
 0042    continue                                                       LOCS388
 0041 continue                                                          LOCS389
      isumat=isumat-1                                                   LOCS390
                                                                        LOCS391
c      do 201 ipr=6,7                                                   LOCS392
c      do 0200 i=1,isumat                                               LOCS393
c         write(ipr,'(4(1x,f10.5))') (coords(j,i),j=1,5)                LOCS394
c 0200 continue                                                         LOCS395
c 201  continue                                                         LOCS396
c                                                                       LOCS397
                                                                        LOCS398
c --- fill matrics mij for each atom - geometry only --                 LOCS399
                                                                        LOCS400
      do 0053 k=1,isumat                                                LOCS401
         do 0054 i=1,3                                                  LOCS402
            do 0055 j=1,3                                               LOCS403
               a(i,j)=0                                                 LOCS404
               m(i,j,k)=coords(j,k)*coords(i,k)                         LOCS405
 0055       continue                                                    LOCS406
c            write(7,'(3f10.4)') (m(i,j,k),j=1,3)                       LOCS407
 0054    continue                                                       LOCS408
 0053 continue                                                          LOCS409
                                                                        LOCS410
c --- sum them up into the sum matrix A ---                             LOCS411
      do 0060 k=1,isumat                                                LOCS412
         do 0061 i=1,3                                                  LOCS413
            do 0062 j=1,3                                               LOCS414
               a(i,j)=a(i,j)+m(i,j,k)                                   LOCS415
 0062       continue                                                    LOCS416
 0061    continue                                                       LOCS417
 0060 continue                                                          LOCS418
                                                                        LOCS419
      do 0071 ipr=6,7                                                   LOCS420
         write(ipr,'(/a/)') ' PRINCIPAL AXIS TRANSFORMATION FOLLOWS '   LOCS421
         write(ipr,*) 'SUMMED GEOMETRY MATRIX Mij=SUM(mij)'             LOCS422
         do 0070 i=1,3                                                  LOCS423
            write(ipr,'(3(1x,f10.5))') (a(i,j),j=1,3)                   LOCS424
 0070    continue                                                       LOCS425
 0071 continue                                                          LOCS426
                                                                        LOCS427
      n=3                                                               LOCS428
      ns=3                                                              LOCS429
                                                                        LOCS430
      call jacobi(a,n,ns,d,v,nrot)                                      LOCS431
                                                                        LOCS432
      do 0081 ipr=6,7                                                   LOCS433
      write (ipr,'(a,i3)')                                              LOCS434
     &   ' NUMBER OF ROTATIONS=',nrot                                   LOCS435
c      do 0080 i=1,3                                                    LOCS436
c         write(ipr,'(3f10.5)') (a(i,j),j=1,3)                          LOCS437
c 0080 continue                                                         LOCS438
      write (ipr,'(a)') ' NORMALIZED EIGENVECTORS '                     LOCS439
      do 0090 i=1,3                                                     LOCS440
         write(ipr,'(3f11.5)') (v(i,j),j=1,3)                           LOCS441
 0090 continue                                                          LOCS442
                                                                        LOCS443
      write (ipr,'(a)') ' EIGENVALUES'                                  LOCS444
      write(ipr,'(3f10.5)') (d(j),j=1,3)                                LOCS445
                                                                        LOCS446
 0081 continue                                                          LOCS447
                                                                        LOCS448
c --- full calculation ---                                              LOCS449
                                                                        LOCS450
      do 0153 k=1,isumat                                                LOCS451
         do 0154 i=1,3                                                  LOCS452
            do 0155 j=1,3                                               LOCS453
               a(i,j)=0                                                 LOCS454
               if (iusech.eq.1) then                                    LOCS455
c ----------------------------------------------------------------------
c                 m(i,j,k)=coords(4,k)*coords(j,k)*coords(i,k)/         LOCS456
c    &                          (coords(4,k)**iweigt)                   LOCS457
c ----------------------------------------------------------------------
               m(i,j,k)=coords(6,k)*coords(5,k)*coords(j,k)*coords(i,k)/LOCS458
     &                          (coords(4,k)**iweigt)                   LOCS459
               elseif (iusech.eq.2) then                                LOCS460
               m(i,j,k)=coords(6,k)*coords(5,k)*coords(5,k)*coords(j,k)*LOCS461
     &                          coords(i,k)/(coords(4,k)**iweigt)       LOCS462
                  if (coords(5,k).eq.0.0) then                          LOCS464
                     write (6,*)'  *** WARNING ****** ZERO CHARGE ****' LOCS465
                     write (7,*)'  *** WARNING ****** ZERO CHARGE ****' LOCS466
                  end if                                                LOCS467
               else                                                     LOCS468
                  m(i,j,k)=coords(j,k)*coords(i,k)/(coords(4,k)**iweigt)LOCS469
               end if                                                   LOCS470
 0155       continue                                                    LOCS471
c           write(7,'(3f10.4)') (m(i,j,k),j=1,3)                        LOCS472
 0154    continue                                                       LOCS473
c       write(*,*) (coords(ik,k),ik=1,6)
c       pause
 0153 continue                                                          LOCS474
                                                                        LOCS475
c --- sum them up into the sum matrix A ---                             LOCS476
      do 0160 k=1,isumat                                                LOCS477
         do 0161 i=1,3                                                  LOCS478
            do 0162 j=1,3                                               LOCS479
               a(i,j)=a(i,j)+m(i,j,k)                                   LOCS480
 0162       continue                                                    LOCS481
 0161    continue                                                       LOCS482
 0160 continue                                                          LOCS483
                                                                        LOCS484
      do 0171 ipr=6,7                                                   LOCS485
      write(ipr,*)                                                      LOCS486
     &  'SUMMED MATRIX Mij=SUM(foc(k)*charge(k)*mij(k)/r(k)**weight)'   LOCS487
      do 0170 i=1,3                                                     LOCS488
         write(ipr,'(3(1x,f10.5))') (a(i,j),j=1,3)                      LOCS489
 0170 continue                                                          LOCS490
 0171 continue                                                          LOCS491
                                                                        LOCS492
      call jacobi(a,3,3,d,v,nrot)                                       LOCS493
                                                                        LOCS494
      do 0181 ipr=6,7                                                   LOCS495
      write (ipr,'(a,i3)')                                              LOCS496
     &   ' NUMBER OF ROTATIONS=',nrot                                   LOCS497
c      do 0180 i=1,3                                                    LOCS498
c         write(ipr,'(3f10.5)') (a(i,j),j=1,3)                          LOCS499
c 0180 continue                                                         LOCS500
      write (ipr,'(a)') ' NORMALIZED EIGENVECTORS '                     LOCS501
      do 0190 i=1,3                                                     LOCS502
         write(ipr,'(3f11.5)') (v(i,j),j=1,3)                           LOCS503
 0190 continue                                                          LOCS504
                                                                        LOCS505
      write (ipr,'(a)') ' EIGENVALUES'                                  LOCS506
      write(ipr,'(3f10.5)') (d(j),j=1,3)                                LOCS507
                                                                        LOCS508
      write (ipr,'(a)') ' NORMALIZED EIGENVALUES '                      LOCS509
      renorm=sqrt(d(1)**2 + d(2)**2 + d(3)**2)                          LOCS510
      write(ipr,'(3f10.5)') d(1)/renorm, d(2)/renorm, d(3)/renorm       LOCS511
      write (ipr,'(a)') ' NORMALIZED EIGENVALUES *sqrt(3) '             LOCS512
      renorm=renorm/sqrt(3.0)                                           LOCS513
      write(ipr,'(3f10.5)') d(1)/renorm, d(2)/renorm, d(3)/renorm       LOCS514
                                                                        LOCS515
                                                                        LOCS516
 0181 continue                                                          LOCS517
                                                                        LOCS518
                                                                        LOCS519
      stop ' Program terminated ok'                                     LOCS520
                                                                        LOCS521
      end                                                               LOCS522
