program patsol c c PROGRAM PATSOL c Patterson Solution Program c c Liang Tong c c Department of Chemistry c University of California c Berkeley, CA 94720 c c Department of Biological Sciences c Purdue University c West Lafayette, IN 47907 c c ------------------ c VAX/Ultrix Version c ------------------ c include 'main.cmn' c character*80 record, rekord character*11 tmpfil (5) data tmpfil/'patsol.tmp', 'patsol1.tmp', 'patsol2.tmp', . 'patsol3.tmp', 'patsol4.tmp'/ c c Disable buffering to the log file c call outbuf c c Initializations c call initss c c copyright stamp c write (log, 30) vershn (1:21) c c Open a scratch file c c open (unit=lin, form='formatted', status='unknown') c open (unit=log, form='formatted', status='unknown') do 1100, nch = 1, 5 open (unit=lscrch, form='FORMATTED', status='NEW', . file=tmpfil(nch), err=1100) go to 100 1100 continue write (log, 20) call exit c 100 read (lin, 10, end=99) record nch = 80 rekord = record call parser (record, rekord, nch) call intprt (record, rekord, nch) if (nch .lt. 0) go to 99 go to 100 c 99 continue c open (unit=lprt, status='NEW', form='FORMATTED', file=prtfil) c if(map) call gthead (0) c call spgsym call patsym call sglvec call ckoren if (nlcaxs .gt. 0) call ncrsym c call chekin c c Patterson map input c if (map) then call patcel call mapinp endif c c Convert all positions and translation components of symmetry c operators to DOUBLED cell grids. c call cvgrid c if (pat) then call patpik if (qmulti) call symdiv call outpat else if (qmulti) call symdiv endif c if (qmulti) then do 1000, i = 1, nnxyz(1) * nnxyz(2) * nnxyz(3) if (ipatmp(i) .lt. 0) ipatmp(i) = - ipatmp(i) 1000 continue endif c c Remove a box around the origin if sum function is being used c if (qsum .and. map) call rmorgn c c NCRshow c if (locsho.gt.0 .and. nncr.gt.0) then call ncrsho write (log, 40) call endgam endif c c NCRSite c if (locsit.gt.0 .and. nncr.gt.0) then call ncrsit write (log, 40) call endgam endif c c SGLshow c if (nsglsh .ne. 0) call sglsho c c DBLshow c if (ndblsh .ne. 0) call dblsho c c SYMshow c if (nsymsh .ne. 0) call symsho c c SGLfit c if (nsglft .ne. 0) call sglfit c c DBLfit c if (ndblft .ne. 0) call dblfit c c Single site solution c if (sgl) then if (nncr .le. 0) call single if (nncr .gt. 0) call ncrsgl endif c c Look for more sites based on input solution lists c if (nsoln .ne. 0) then call fillup if (sol) call solutn endif c c Two atom search - use single site solution lists c if (dbl) then if (nncr .le. 0) then call double if (sol) call triple else call ncrdbl endif endif c c Two atom search - use Patterson peaks as cross vectors c if (patx) then call patdbl if (sol) call triple endif c c Two atom search - use input cross vectors c if (ncross .ne. 0) then call crsdbl if (sol) call triple endif c c All done c write (log, 40) call endgam c 10 format (a) 20 format (' PATSOL> *Fatal* ', 'Unable to open a scratch file', /, . ' PROGRAM ABORT') 30 format (' PATSOL> ', /, ' PATSOL> ', 'Program ', a, . ', by Liang Tong.', /, ' PATSOL> ', 'This program may ', . 'not be re-distributed in any form.', /, ' PATSOL> ') 40 format (' PATSOL> ', /, ' PATSOL> ', 'Program finished ', . 'successfully.') c end c c----------------------------------------------------------------------- c subroutine ccpinp c include 'main.cmn' c write (log, 10) call endgam c 10 format (' CCPINP> *Fatal* ', 'CCP4 input routine is not present', . /, ' CCPINP> ', 'Re-compile program with the correct ', . 'routine') c return end c c----------------------------------------------------------------------- c subroutine chekin c include 'main.cmn' c if (ncell(1)*ncell(2)*ncell(3) .eq. 0) then write (log, 10) call endgam endif c if (map .and. nxyz(1)*nxyz(2)*nxyz(3).eq.0) then write (log, 40) map = .false. endif c if (isym .eq. 0) then write (log, 20) isym = 1 do 1100, i = 1, 3 its (i, 1) = 0 do 1200, j = 1, 3 ifs (i, j, 1) = 0 1200 continue ifs (i, i, 1) = 1 1100 continue endif c if (.not. map) then pat = .false. patx = .false. sol = .false. nsoln = 0 ncross = 0 nsglft = 0 ndblft = 0 if (nncr .le. 0) then sgl = .false. dbl = .false. endif endif c if (nncr .ne. 0) then if (sgl) dbl = .false. if (dbl) sgl = .false. if (dbl) qspher = .false. qsum = .true. sol = .false. patx = .false. nsoln = 0 ncross = 0 nsglsh = 0 ndblsh = 0 nsymsh = 0 nsglft = 0 ndblft = 0 write (log, 70) endif c if (nncr .le. 0) then if (dbl) sgl = .true. if (patx) sgl = .true. if (sol) sgl = .true. if (ncross .ne. 0) sgl = .true. if (nsoln .ne. 0) sgl = .true. if (.not. qsum) write (log, 50) if (qsum) write (log, 60) c locsho = 0 locsit = 0 c if (isym .eq. 1) then sgl = .false. nsglpk = 1 do 1000, i = 1, 3 isgxyz (i, 1) = 1 isglpk (i, 1) = 1 1000 continue isgxyz (4, 1) = 999 isglpk (4, 1) = 999 write (log, 30) endif endif c c Convert the input map region into grid region. Array nnxyz c is used for addressing the Patterson map grid points later. c do 1300, i = 1, 3 if (nxyz(i) .gt. 1000) then x = float(nxyz (i)) / 1000000.0 x = x * float (ncell (i)) nxyz (i) = x + 1.5 endif nnxyz (i) = nxyz (i) 1300 continue c 10 format (' CHEKIN> *Fatal* ', 'Cell grid can not be zero', . /, ' PROGRAM ABORT') 20 format (' CHEKIN> ', 'No symmetry operator input, P1 is assumed') 30 format (' CHEKIN> ', 'Space group is P1, no single-site search ', . 'will be carried out') 40 format (' CHEKIN> *Warning* ', 'Map region is not defined', /, . ' CHEKIN> Map will not be read in') 50 format (' CHEKIN> ', 'Minimum function will be used.') 60 format (' CHEKIN> ', 'Sum function will be used.') 70 format (' CHEKIN> ', 'Patterson interpretation with local ', . 'symmetry', /, ' CHEKIN> ', 'Sum function will be used') c return end c c----------------------------------------------------------------------- c subroutine ckaxis (is, th1, th2, th3) c include 'main.cmn' c dimension rotmat (3, 3) c do 1000, i = 1, 3 do 1100, j = 1, 3 rotmat (i, j) = rlcrot (i, j, is) 1100 continue 1000 continue c th1 = -1000.0 if (polar .eq. 'XZK') then call polar1 (th1, th2, th3, rotmat) call ckpol1 (th1, th2, th3, rotmat) else if (polar .eq. 'XYK') then call polar2 (th1, th2, th3, rotmat) call ckpol2 (th1, th2, th3, rotmat) else write (log, 10) 'Polar', polar endif c if (th3 .lt. 0.5) then mlcaxs (is) = 0 return endif c if (is .lt. nlcaxs) then n = 360.0 / th3 + 0.5 if (.not. qlcxpd) then mlcaxs (is) = n return endif dth = abs (360.0 - float(n) * th3) if (dth .gt. 0.5) then write (log, 20) th3 call endgam endif c th3 = 360.0 / float (n) if (mlcaxs(is).ne.0 .and. n.ne.mlcaxs (is)) then write (log, 30) n, mlcaxs (is), n mlcaxs (is) = n else if (mlcaxs (is) .eq. 0) then mlcaxs (is) = n endif endif c 10 format (/, ' CKAXIS> *Fatal* ', 'Unknown ', a, ' Angle Type <', . a, '>', /, ' PROGRAM ABORT') 20 format (/, ' CKAXIS> *Fatal* ', 'Input axis has rotation of', . f6.1, . ' degrees -- Not divisible to 360.0', /, ' PROGRAM ABORT') 30 format (' CKAXIS> *Warning* ', 'Input rotation has ', i2, . ' fold symmetry, rather than that of ', i2, ' as input.', . /, ' CKAXIS> ', i2, ' fold symmetry will be used') c return end c c----------------------------------------------------------------------- c subroutine cklmts (iopt) c include 'main.cmn' c call nupage c if (iopt .eq. 1) then do 1000, n = 1, 3 if (slimts(3, n) .lt. 0.0005) then slimts (3, n) = 0.0 slimts (2, n) = slimts (1, n) nlimts (n) = 1 else nlimts (n) = (slimts (2, n) - slimts (1, n)) . / slimts (3, n) + 1.5 slimts (2, n) = slimts (1, n) + float (nlimts (n) - 1) . * slimts (3, n) endif 1000 continue c if (sgl) write (lprt, 10) if (dbl) write (lprt, 20) do 1100, n = 1, 3 write (lprt, 30) n, (slimts(i, n), i=1, 3), nlimts(n) 1100 continue endif c if (dbl) qspher = .false. if (qspher) write (lprt, 60) polar if (.not.qspher .and. sgl) write (lprt, 70) if (dbl) write (lprt, 80) if (qckall) write (lprt, 40) if (.not. qckall) write (lprt, 50) if (qmuncr) write (lprt, 90) if (.not. qmuncr) write (lprt, 100) if (sgl .and. .not.qspher) then write (lprt, 210) sizmol sizmol = sizmol * sizmol else if (dbl) then write (lprt, 220) sizmol sizmol = 4 * sizmol * sizmol endif if (qdump .and. nlimts(3).gt.40) qdump = .false. if (qdump) write (lprt, 180) if (.not. qdump) write (lprt, 190) write (lprt, 160) icutlo - 2000 write (lprt, 200) iorign c if (sgl) return c c Find the input crystallographic symmetry operator c if (isym .le. 1) then write (log, 110) call endgam endif if (isymop .gt. nsym) then write (log, 120) isymop isymop = 2 endif c do 1200, i = 1, 3 dsymk (i) = float(its (i, isymop)) / float (ncell (i)) do 1300, j = 1, 3 tsymk (i, j) = ifs (i, j, isymop) tsymk1 (i, j) = ifs (i, j, isymop) 1300 continue tsymk1 (i, i) = tsymk1 (i, i) - 1.0 1200 continue c write (lprt, 170) pinput write (lprt, 130) isymop write (lprt, 140) dsymk write (lprt, 150) ((tsymk(i, j), j=1, 3), i=1, 3), . ((tsymk1(i, j), j=1, 3), i=1, 3) c 10 format (5x, 'Self-vector search with the specified local ', . 'symmetry', //, 5x, 'Search Region -- ', //) 20 format (5x, 'Cross-vector search with the specified local ', . 'symmetry', //, 5x, 'Search Region -- ', //) 30 format (10x, i2, 3x, 'Limits : ', 2x, 3f10.4, . 3x, 'Number of points : ', 2x, i4) 40 format (/, 5x, 'All vectors will be checked') 50 format (/, 5x, 'Only those vectors between a reference ', . 'position and other positions will be checked') 60 format (/, 5x, 'Search will be done in spherical polar ', . 'coordinates (1=Phi, 2=Psi, 3=R).', /, . 5x, 'Polar angle convention is ', a) 70 format (/, 5x, 'Search will be done in Cartesian coordinates ', . '(1=P, 2=Q, 3=R).') 80 format (/, 5x, 'Search will be done in fractional coordinates ', . '(1=a, 2=b, 3=c).') 90 format (/, 5x, 'Multiplicity correction will be carried out.') 100 format (/, 5x, 'Multiplicity correction will not be carried ', . 'out.') 110 format (' CKLMTS> *Fatal* ', 'Space group is P1. Cross-vector ', . 'search is not necessary.', /, ' PROGRAM ABORT') 120 format (' CKLMTS> *Warning* ', 'Crystallographic symmetry ', . 'operator #', i2, ' is not found.', /, ' CKLMTS> ', . 'Symmetry operator #2 will be used.') 130 format (/, 5x, 'Crystallographic symmetry operator #', i2, . ' will be used in the cross-vector search.') 140 format (/, 10x, 'Translation Elements : ', 3f8.5) 150 format (10x, 'Rotation Elements : ', 3(3f8.3, 3x), /, . 10x, 'Rotation - Identity : ', 3(3f8.3, 3x)) 160 format (/, 5x, 'Individual vector cut-off -- Cutlo = ', i4) 170 format (/, 5x, 'Position from self-vector search (P input) : ', . 3f10.4) 180 format (/, 5x, 'The vector search map will be dumped to the ', . 'print file.') 190 format (/, 5x, 'The vector search map will not be dumped to the ', . 'print file.') 200 format (/, 5x, 'Size of box around the origin ', 3i4, /, 5x, . 'Vectors within this box will not be considered in the ', . 'calculation.') 210 format (/, 5x, 'Maximum distance of search point to origin ', . f10.2, ' Angstroms.') 220 format (/, 5x, 'Radius of molecular assembly for packing ', . 'calculations ', f10.2, ' Angstroms.') c return end c c----------------------------------------------------------------------- c subroutine ckoren c c This routine locates the possible alternative origins of the space c group. Positions that are related by these origin translations c share the same self-vectors and hence can not be distinguished from c each other. The routine also determines whether (x, y, z) and c (-x, -y, -z) share the same Harker vectors. c c All these are carried out by checking the identity of the single-site c vector matrices. c c SUBROUTINE CALLED -- SMPATT c c IWORKS -- c 13 - 16 : position (x, y, z), or (-x, -y, -z), with OR. c 21 - 28 : the 8 possible origins to check. c 1 - 4 : self-vector position corresponding to 13-16. c 5 - 8 : self-vector position corresponding to (x, y, z) c include 'main.cmn' c c First -- alternate origins? c do 1000, i = 1, 3 do 1100, j = 1, 4 if (i .eq. j) then iworks (i, j+12) = 1 else iworks (i, j+12) = 0 endif 1100 continue 1000 continue c do 1200, i = 21, 28 do 1300, k = 1, 3 iworks (k, i) = 0 1300 continue 1200 continue c c 22 - (1/2, 0, 0) iworks (1, 22) = 6 c c 23 - (1/2, 1/2, 0) iworks (1, 23) = 6 iworks (2, 23) = 6 c c 24 - (0, 1/2, 0) iworks (2, 24) = 6 c c 25 - (0, 0, 1/2) iworks (3, 25) = 6 c c 26 - (1/2, 0, 1/2) iworks (1, 26) = 6 iworks (3, 26) = 6 c c 27 - (1/2, 1/2, 1/2) iworks (1, 27) = 6 iworks (2, 27) = 6 iworks (3, 27) = 6 c c 28 - (0, 1/2, 1/2) iworks (2, 28) = 6 iworks (3, 28) = 6 c write (lprt, 10) m = 0 do 100, i = 1, 8 c if (icent.eq.2 .and. i.ge.5) go to 100 if (icent.eq.3 .and. i.ge.5) go to 100 if (icent.eq.4 .and. (i.eq.3.or.i.eq.4.or.i.ge.7)) go to 100 if (icent.eq.5 .and. i.ge.3) go to 100 if (icent.eq.6 .and. i.ge.5) go to 100 c do 1400, j = 1, 3 iworks (j, 16) = iworks (j, i+20) if (fixxyz(j) .and. iworks(j, 16).ne.0) go to 100 1400 continue c do 1500, j = 1, nsgl do 1600, k = 1, 3 do 1700, L = 1, 4 iworks (k, L) = 0 if (L .eq. 4) iworks(k, L) = isglt(k, j) do 1800, n = 1, 3 iworks(k, L) = iworks(k, L) + . isglr(k, n, j) * iworks(n, L+12) 1800 continue 1700 continue 1600 continue do 1900, k = 1, 3 iworks (k, 8) = isglt (k, j) do 2000, L = 1, 3 iworks (k, L+4) = isglr (k, L, j) 2000 continue 1900 continue call smpatt (isame, 1) if (isame .eq. 0) go to 100 1500 continue m = m + 1 do 2100, j = 1, 3 isglor (j, m) = iworks (j, i+20) 2100 continue write (lprt, 20) m, (isglor(j, m), j=1, 3) 100 continue c c Now enantiomorph c write (lprt, *) do 2200, j = 1, 3 do 2300, k = 1, 4 if (j .eq. k) then iworks (j, k+12) = -1 else iworks (j, k+12) = 0 endif 2300 continue 2200 continue do 2400, j = 1, nsgl do 2500, k = 1, 3 do 2600, L = 1, 4 iworks (k, L) = 0 if (L .eq. 4) iworks(k, L) = isglt(k, j) do 2700, n = 1, 3 iworks(k, L) = iworks(k, L) + . isglr(k, n, j) * iworks(n, L+12) 2700 continue 2600 continue 2500 continue do 2800, k = 1, 3 iworks (k, 8) = isglt (k, j) do 2900, L = 1, 3 iworks (k, L+4) = isglr (k, L, j) 2900 continue 2800 continue call smpatt (isame, 1) if (isame .eq. 0) go to 200 2400 continue write (lprt, 30) return c 200 write (lprt, 40) enan = .false. c 10 format (5x, 'The positions related by the following alternate ', . 'origins share the same self-vectors ', /) 20 format (10x, i4, 5x, 3i3) 30 format (5x, 'Positions (x, y, z) and (-x, -y, -z) share the ', . 'same self-vectors') 40 format (5x, 'Positions (x, y, z) and (-x, -y, -z) do not ', . 'share the same self-vectors') c return end c c----------------------------------------------------------------------- c subroutine ckpack (x, y, z, ipack) c include 'main.cmn' c dimension xyz(3), xyzn(3) c ipack = 0 c xyz (1) = x xyz (2) = y xyz (3) = z c do 1000, i = 2, nsym c txyz = 0.0 do 1100, j = 1, 3 c xyzn (j) = float (its (j, i)) / float (ncell (j)) do 1200, k = 1, 3 xyzn (j) = xyzn (j) + float (ifs (j, k, i)) * xyz (k) 1200 continue c xyzn (j) = mod (xyzn (j) - xyz (j), 1.0) if (xyzn(j) .lt. 0.0) xyzn (j) = xyzn (j) + 1.0 if (xyzn(j) .gt. 0.5) xyzn (j) = xyzn (j) - 1.0 txyz = txyz + xyzn (j) * xyzn (j) * cell (j, 2) 1100 continue c txyz = txyz + cell (4, 2) * xyzn (2) * xyzn (3) . + cell (5, 2) * xyzn (1) * xyzn (3) . + cell (6, 2) * xyzn (1) * xyzn (2) c if (txyz .lt. sizmol) then ipack = 1 return endif c 1000 continue c return end c c----------------------------------------------------------------------- c subroutine ckpol1 (th1, th2, th3, rotmat) c include 'main.cmn' c dimension rotmat(3, 3), rotnew(3, 3) c c This checks to make sure the extraction of angles from the matrix c has been done correctly. c c call polar1 (th1, th2, th3, rotnew) c do 1000, i = 1, 3 do 1100, j = 1, 3 s1 = abs (rotmat (i, j) - rotnew (i, j)) if (s1 .gt. 0.003) then write (log, 10) ((rotmat(m, n), n=1, 3), m=1, 3) write (log, 20) th1, th2, th3 write (log, 30) ((rotnew(m, n), n=1, 3), m=1, 3) return endif 1100 continue 1000 continue c 10 format (' CKPOL1> *Warning* ', 'Error extracting angles from ', . 'this rotation matrix', /, ' CKPOL1> ', 3(2x, 3f7.3)) 20 format (' CKPOL1> ', 'Extracted angles are : ', 3f10.3) 30 format (' CKPOL1> ', 'New rotation matrix is', /, . ' CKPOL1> ', 3(2x, 3f7.3)) c return end c c----------------------------------------------------------------------- c subroutine ckpol2 (th1, th2, th3, rotmat) c include 'main.cmn' c dimension rotmat(3, 3), rotnew(3, 3) c c This checks to make sure the extraction of angles from the matrix c has been done correctly. c c call polar2 (th1, th2, th3, rotnew) c do 1000, i = 1, 3 do 1100, j = 1, 3 s1 = abs (rotmat (i, j) - rotnew (i, j)) if (s1 .gt. 0.003) then write (log, 10) ((rotmat(m, n), n=1, 3), m=1, 3) write (log, 20) th1, th2, th3 write (log, 30) ((rotnew(m, n), n=1, 3), m=1, 3) return endif 1100 continue 1000 continue c 10 format (' CKPOL2> *Warning* ', 'Error extracting angles from ', . 'this rotation matrix', /, ' CKPOL2> ', 3(2x, 3f7.3)) 20 format (' CKPOL2> ', 'Extracted angles are : ', 3f10.3) 30 format (' CKPOL2> ', 'New rotation matrix is', /, . ' CKPOL2> ', 3(2x, 3f7.3)) c return end c c----------------------------------------------------------------------- c subroutine crsdbl c c This routine carries out a two-atom search based on a list of input c cross vectors. c c SUBROUTINES CALLED : GTSGLE, XVECTR, SAVDBL, DBLREF, OUTDBL c c IWORKS -- c 1 : site #1 c 2 : site #1 + cross-vector c include 'main.cmn' c write (log, 10) ndbl = 0 jdblct = idblct - 2000 if (jdblct .gt. 0) then jdblct = jdblct / 3 else if (jdblct .lt. 0) then jdblct = max (-500, jdblct * 2) else jdblct = -200 endif jdblct = jdblct + 2000 c do 1000, i = 1, maxsgl if (isgxyz(4, i) .lt. 0) go to 200 do 1100, j = 1, 3 iworks (j, 1) = isgxyz (j, i) 1100 continue c do 1400, j = 1, ncross do 1200, k = 1, 3 iworks (k, 2) = mod(iworks(k, 1) + icross(k, j) - 2, . ncell(k)) iworks (k, 2) = iworks (k, 2) + 1 1200 continue c call xvectr (mn, minx, jdblct) if (mn .lt. jdblct) go to 1400 c call spgasu (iasu, 1) if (iasu .eq. 1) go to 1400 c call gtsgle (isame) if (isame .eq. 0) go to 1400 c do 1300, ii = 1, 3 jworks (ii, 1) = iworks (ii, 1) jworks (ii, 2) = iworks (ii, 2) 1300 continue jworks (4, 1) = isgxyz (4, i) jworks (4, 2) = isgxyz (4, isame) call savdbl (mn) c 1400 continue 1000 continue c 200 continue c call dblref call outdbl (2) c 10 format (' CRSDBL> ', 'Two-atom search with input cross-vectors') c return end c c----------------------------------------------------------------------- c subroutine cvgrid c c This routine converts the translational components in the symmetry c operators from the 12 base to the base of the actual number of grid c points along the cell edges. One important thing is to check to c make sure that the number of grid points is divisible by the c fractional translation to insure accurate representation of the c symmetry. c include 'main.cmn' c c Check NCELL is divisible by the translation components. c Note the ITS and IPTS are based on doubled NCELL. c do 1000, i = 1, nsym do 1100, j = 1, 3 if (its (j, i).eq.2 .and. mod(ncell(j), 6).ne.0) then write (log, 10) ncell(j), 6 call endgam else if (its (j, i).eq.3 .and. mod(ncell(j), 4).ne.0) then write (log, 10) ncell(j), 4 call endgam else if (its (j, i).eq.4 .and. mod(ncell(j), 3).ne.0) then write (log, 10) ncell(j), 3 call endgam else if (its (j, i).eq.6 .and. mod(ncell(j), 2).ne.0) then write (log, 10) ncell(j), 2 call endgam else if (its (j, i).eq.8 .and. mod(ncell(j), 3).ne.0) then write (log, 10) ncell(j), 3 call endgam else if (its (j, i).eq.9 .and. mod(ncell(j), 4).ne.0) then write (log, 10) ncell(j), 4 call endgam else if (its (j, i).eq.10 .and. mod(ncell(j), 6).ne.0) then write (log, 10) ncell(j), 6 call endgam endif c its (j, i) = ncell (j) * its (j, i) / 6 c 1100 continue 1000 continue c do 1200, i = 1, npsym do 1300, j = 1, 3 if (ipts (j, i).eq.4 .and. mod(ncell(j), 3).ne.0) then write (log, 10) ncell(j), 3 call endgam else if (ipts (j, i).eq.8 .and. mod(ncell(j), 3).ne.0) then write (log, 10) ncell(j), 3 call endgam else if (ipts (j, i).eq.6 .and. mod(ncell(j), 2).ne.0) then write (log, 10) ncell(j), 2 call endgam endif c ipts (j, i) = ncell (j) * ipts (j, i) / 6 c 1300 continue 1200 continue c do 1400, i = 1, nsgl do 1500, j = 1, 3 isglt (j, i) = ncell (j) * isglt (j, i) / 6 1500 continue 1400 continue c do 1600, i = 1, 8 if (isglor(1, i) .ge. 0) then do 1700, j = 1, 3 isglor (j, i) = ncell (j) * isglor (j, i) / 6 1700 continue endif 1600 continue c c Double the cell grids c do 1900, i = 1, 3 ncell (i) = 2 * ncell (i) nxyz (i) = 2 * nxyz (i) ntol (i) = 2 * ntol (i) iorign (i) = 2 * iorign (i) 1900 continue c c Isglsh c do 2000, i = 1, nsglsh do 2100, j = 1, 3 if (abs(isglsh(j, i)).lt.1000 .and. . abs(isglsh(j, i)).gt.0) then isglsh (j, i) = 2 * isglsh (j, i) - 1 else x = float(isglsh (j, i)) / 1000000.0 x = x * float (ncell (j)) if (x .ge. 0.0) then isglsh (j, i) = x + 1.5 else isglsh (j, i) = x + 0.5 endif endif 2100 continue 2000 continue c c Idblsh c do 2200, i = 1, ndblsh do 2300, j = 1, 6 if (abs(idblsh(j, i)).lt.1000 .and. . abs(idblsh(j, i)).gt.0) then idblsh (j, i) = 2 * idblsh (j, i) - 1 else x = float(idblsh (j, i)) / 1000000.0 jj = j if (jj .gt. 3) jj = jj - 3 x = x * float (ncell (jj)) if (x .ge. 0.0) then idblsh (j, i) = x + 1.5 else idblsh (j, i) = x + 0.5 endif endif 2300 continue 2200 continue c c Isymsh c do 2900, i = 1, nsymsh do 3000, j = 1, 3 if (abs(isymsh(j, i)).lt.1000 .and. . abs(isymsh(j, i)).gt.0) then isymsh (j, i) = 2 * isymsh (j, i) - 1 else x = float(isymsh (j, i)) / 1000000.0 x = x * float (ncell (j)) if (x .ge. 0.0) then isymsh (j, i) = x + 1.5 else isymsh (j, i) = x + 0.5 endif endif 3000 continue 2900 continue c c Isglft c do 3300, i = 1, nsglft do 3400, j = 1, 3 if (abs(isglft(j, i)).lt.1000 .and. . abs(isglft(j, i)).gt.0) then isglft (j, i) = 2 * isglft (j, i) - 1 else x = float(isglft (j, i)) / 1000000.0 x = x * float (ncell (j)) if (x .ge. 0.0) then isglft (j, i) = x + 1.5 else isglft (j, i) = x + 0.5 endif endif 3400 continue 3300 continue c c Idblft c do 3100, i = 1, ndblft do 3200, j = 1, 6 if (abs(idblft(j, i)).lt.1000 .and. . abs(idblft(j, i)).gt.0) then idblft (j, i) = 2 * idblft (j, i) - 1 else x = float(idblft (j, i)) / 1000000.0 jj = j if (jj .gt. 3) jj = jj - 3 x = x * float (ncell (jj)) if (x .ge. 0.0) then idblft (j, i) = x + 1.5 else idblft (j, i) = x + 0.5 endif endif 3200 continue 3100 continue c c Icross c do 2400, i = 1, ncross do 2500, j = 1, 3 if (abs(icross(j, i)).lt.1000 .and. . abs(icross(j, i)).gt.0) then icross (j, i) = 2 * icross (j, i) - 1 else x = float(icross (j, i)) / 1000000.0 x = x * float (ncell (j)) if (x .ge. 0.0) then icross (j, i) = x + 1.5 else icross (j, i) = x + 0.5 endif endif 2500 continue 2400 continue c c Isoln c do 2600, i = 1, nsoln do 2700, j = 1, 3 if (abs(isoln(j, i)).lt.1000 .and. . abs(isoln(j, i)).gt.0) then isoln (j, i) = 2 * isoln (j, i) - 1 else x = float(isoln (j, i)) / 1000000.0 x = x * float (ncell (j)) if (x .ge. 0.0) then isoln (j, i) = x + 1.5 else isoln (j, i) = x + 0.5 endif endif 2700 continue 2600 continue c c Nsxyz c do 2800, i = 1, 3 if (nsxyz(i) .le. 0) then nsxyz (i) = 0 else if (nsxyz(i) .lt. 1000) then nsxyz (i) = nsxyz (i) * 2 else x = float(nsxyz (i)) / 1000000.0 x = x * float (ncell (i)) nsxyz (i) = x + 1.5 if (nsxyz(i) .gt. ncell(i)) nsxyz (i) = ncell (i) endif 2800 continue c c Now convert the cut-offs to sigma=100 base c txyz = float (isglct) / 10000. if (abs(txyz) .lt. 10.0) txyz = 100.0 * txyz isglct = txyz + sign (0.5, txyz) c txyz = float (idblct) / 10000. if (txyz .lt. 10.0) txyz = 100.0 * txyz idblct = txyz + sign (0.5, txyz) c txyz = float (ipatct) / 10000. if (txyz .lt. 10.0) txyz = 100.0 * txyz ipatct = txyz + sign (0.5, txyz) c txyz = float (itplsg) / 10000. if (abs(txyz) .lt. 10.0) txyz = 100.0 * txyz itplsg = txyz + sign (0.5, txyz) c txyz = float (itpldb) / 10000. if (abs(txyz) .lt. 10.0) txyz = 100.0 * txyz itpldb = txyz + sign (0.5, txyz) c txyz = float (icutlo) / 10000. if (abs(txyz) .lt. 10.0) txyz = 100.0 * txyz icutlo = txyz + sign (0.5, txyz) c isglct = isglct + 2000 idblct = idblct + 2000 ipatct = ipatct + 2000 itplsg = itplsg + 2000 itpldb = itpldb + 2000 icutlo = icutlo + 2000 c write (log, 20) ncell write (log, 30) nxyz write (log, 40) ntol c 10 format (' CVGRID> *Fatal* ', 'Number of grid points along ', . ' axis ', i4, ' is not divisible by ', i2, /, . ' PROGRAM ABORT') 20 format (' CVGRID> ', 'Doubled unit cell grid ', 3i4) 30 format (' CVGRID> ', 'Doubled input Patterson map region ', 3i4) 40 format (' CVGRID> ', 'Tolerance limits ', 3i3) c return end c c----------------------------------------------------------------------- c subroutine dblfit c c This routine attempts to maximize the cross-vector peak height c by jiggling the input sites. c c SUBROUTINES CALLED -- NUPAGE, SGLSHO, XVECTR c c IWORKS -- c 1 : single site position (i, j, k) c 2 : a single-site vector for (i, j, k) c 3 : a known single-site vector for (i, j, k) c JWORKS -- c 1 - jwork : the single-site vectors and their redundancies c include 'main.cmn' c do 1000, i = 1, ndblft c call nupage write (lprt, 10) (k, (float (idblft (j + (k - 1) * 3, i) - 1) / . float (ncell (j)), j = 1, 3), k = 1, 2) c do 2300, k = 1, 2 do 1100, j = 1, 3 iworks (j, 1) = idblft (j + (k - 1) * 3, i) 1100 continue call svectr (ival, minhk, 0) write(lprt, 70) k do 2400, j = 1, jwork write (lprt, 30) j, (jworks(m1, j), m1=1, 3), . (float(jworks(m1, j)-1)/float(ncell(m1)), m1=1, 3), . jworks(4, j) - 2000 2400 continue write (lprt, 80) ival - 2000 2300 continue c do 1200, j = 1, 3 iworks (j, 1) = idblft (j, i) iworks (j, 2) = idblft (j + 3, i) 1200 continue c call xvectr (ival, minx, 0) c write(lprt, 20) do 1300, j = 1, jwork write (lprt, 30) j, (jworks(k, j), k=1, 3), . (float(jworks(k, j)-1)/float(ncell(k)), k=1, 3), . jworks(4, j) - 2000 1300 continue write (lprt, 40) ival - 2000 c iworks (1, 100) = -10000 do 1400, L1 = - ntol (1) / 2, ntol (1) / 2 if (fixxyz(1) .and. L1.ne.0) go to 1400 do 1500, m1 = - ntol (2) / 2, ntol (2) / 2 if (fixxyz(2) .and. m1.ne.0) go to 1500 do 1600, n1 = - ntol (3) / 2, ntol (3) / 2 if (fixxyz(3) .and. n1.ne.0) go to 1600 c do 1700, L2 = - ntol (1) / 2, ntol (1) / 2 do 1800, m2 = - ntol (2) / 2, ntol (2) / 2 do 1900, n2 = - ntol (3) / 2, ntol (3) / 2 c if (L1.eq.0 .and. m1.eq.0 .and. n1.eq.0 . .and. L2.eq.0 .and. m2.eq.0 .and. n2.eq.0) . go to 1900 c iworks (1, 1) = idblft (1, i) + L1 iworks (2, 1) = idblft (2, i) + m1 iworks (3, 1) = idblft (3, i) + n1 iworks (1, 2) = idblft (4, i) + L2 iworks (2, 2) = idblft (5, i) + m2 iworks (3, 2) = idblft (6, i) + n2 c call xvectr (k, minx, 0) if (k .le. ival) go to 1900 c iworks (1, 100) = iworks (1, 1) iworks (2, 100) = iworks (2, 1) iworks (3, 100) = iworks (3, 1) iworks (1, 101) = iworks (1, 2) iworks (2, 101) = iworks (2, 2) iworks (3, 101) = iworks (3, 2) ival = k c 1900 continue 1800 continue 1700 continue 1600 continue 1500 continue 1400 continue c if (iworks(1, 100) .gt. -999) then call nupage write (lprt, 50) (k, (float (iworks (j, k + 99) - 1) / . float (ncell (j)), j = 1, 3), k = 1, 2) c do 2500, k = 1, 2 do 2600, j = 1, 3 iworks (j, 1) = iworks (j, 99 + k) 2600 continue call svectr (ival, minhk, 0) write(lprt, 70) k do 2700, j = 1, jwork write (lprt, 30) j, (jworks(m1, j), m1=1, 3), . (float(jworks(m1, j)-1)/float(ncell(m1)), m1=1, 3), . jworks(4, j) - 2000 2700 continue write (lprt, 80) ival - 2000 2500 continue c do 2100, j = 1, 3 iworks (j, 1) = iworks (j, 100) iworks (j, 2) = iworks (j, 101) 2100 continue call xvectr (ival, minx, 0) c write(lprt, 20) do 2200, j = 1, jwork write (lprt, 30) j, (jworks(k, j), k=1, 3), . (float(jworks(k, j)-1)/float(ncell(k)), k=1, 3), . jworks(4, j) - 2000 2200 continue write (lprt, 40) ival - 2000 else write (lprt, 60) endif 1000 continue c 10 format (5x, 'Maximizing the cross-vector height for the input ', . 'two-atom solution', /, 10x, i3, 3x, 3f7.3, /, . 10x, i3, 3x, 3f7.3) 20 format (/, 5x, 'Cross vectors between the two sites ', //, . 9x, 'No.', 3x, ' I ', ' J ', ' K ', 6x, 'u', 6x, 'v', . 6x, 'w', 6x, 'Height') 30 format (10x, i2, 3x, 3i4, 3x, 3f7.3, 3x, i5) 40 format (/, 10x, 'Cross vector search map value ', i5) 50 format (5x, 'The two-atom solution after refinement ', . /, 10x, i3, 3x, 3f7.3, /, 10x, i3, 3x, 3f7.3) 60 format (//, 5x, 'No improvement in cross vector height from ', . 'refinement') 70 format (/, 5x, 'Self vectors for site #', i1, //, . 9x, 'No.', 3x, ' I ', ' J ', ' K ', 6x, 'u', 6x, 'v', . 6x, 'w', 6x, 'Height') 80 format (/, 10x, 'Self vector search map value ', i5) c return end c c----------------------------------------------------------------------- c subroutine dblref c c This routine maximizes the cross-vector heights for two-atom c solutions by jiggling the individual sites by half a grid point. c The combination that gives the highest cross-vector (irrespective of c the self-vector heights) will be saved by the program. c c SUBROUTINES CALLED : SVECTR, XVECTR c c IWORKS -- c 11 : atom position #1 c 12 : atom position #2 c 1 : position #1 modified by a little c 2 : position #2 modified by a little c include 'main.cmn' c do 1000, i = 1, ndbl c ixht = idblpk (6, i) do 1100, j = 1, 3 iworks (j, 11) = idblpk (j, i) / 1000 iworks (j, 12) = mod (idblpk (j, i), 1000) 1100 continue c do 1200, L1 = -1, 1 if (fixxyz(1) .and. L1.ne.0) go to 1200 do 1300, m1 = -1, 1 if (fixxyz(2) .and. m1.ne.0) go to 1300 do 1400, n1 = -1, 1 if (fixxyz(3) .and. n1.ne.0) go to 1400 c do 1500, L2 = -1, 1 do 1600, m2 = -1, 1 do 1700, n2 = -1, 1 c if (L1.eq.0 .and. m1.eq.0 .and. n1.eq.0 . .and. L2.eq.0 .and. m2.eq.0 .and. n2.eq.0) . go to 1700 c iworks (1, 1) = iworks (1, 11) + L1 iworks (2, 1) = iworks (2, 11) + m1 iworks (3, 1) = iworks (3, 11) + n1 iworks (1, 2) = iworks (1, 12) + L2 iworks (2, 2) = iworks (2, 12) + m2 iworks (3, 2) = iworks (3, 12) + n2 c call xvectr (mn, minx, idblct) if (mn .le. ixht) go to 1700 c ixht = mn LL1 = L1 mm1 = m1 nn1 = n1 LL2 = L2 mm2 = m2 nn2 = n2 c 1700 continue 1600 continue 1500 continue 1400 continue 1300 continue 1200 continue c if (ixht .le. idblpk(6, i)) go to 1000 c iworks (1, 1) = iworks (1, 11) + LL1 iworks (2, 1) = iworks (2, 11) + mm1 iworks (3, 1) = iworks (3, 11) + nn1 call svectr (mn1, minhk, 0) c iworks (1, 1) = iworks (1, 12) + LL2 iworks (2, 1) = iworks (2, 12) + mm2 iworks (3, 1) = iworks (3, 12) + nn2 call svectr (mn2, minhk, 0) c iworks (1, 1) = iworks (1, 11) + LL1 iworks (2, 1) = iworks (2, 11) + mm1 iworks (3, 1) = iworks (3, 11) + nn1 iworks (1, 2) = iworks (1, 12) + LL2 iworks (2, 2) = iworks (2, 12) + mm2 iworks (3, 2) = iworks (3, 12) + nn2 c do 1800, j = 1, 2 do 1900, k = 1, 3 iworks (k, j) = mod (iworks (k, j) . + ncell (k), ncell (k)) if (iworks(k, j) .eq. 0) iworks (k, j) = ncell (k) jworks (k, j) = iworks (k, j) 1900 continue 1800 continue jworks (4, 1) = mn1 jworks (4, 2) = mn2 c call savdbl (ixht) c 1000 continue c c Now go over all the sites to remove those with cross-vector heights c less than Idblct. c ixht = ndbl ndbl = 0 c do 2000, i = 1, ixht if (idblpk(6, i) .lt. idblct) go to 2000 ndbl = ndbl + 1 do 2100, j = 1, 6 idblpk (j, ndbl) = idblpk (j, i) 2100 continue 2000 continue c return end c c----------------------------------------------------------------------- c subroutine dblsgl (i1, i2, i3, i4) c c This routine compares a site in IWORKS (*, 1) with the single-site c solution peaks, and returns in I1 the symmetry-operation, I2 the c origin translation, I3 the enantiomorph, and I4 the single-site c number that the site corresponds to. c c Note that the DO loops go over the space-group symmetry operations c first, then the centering operations, and finally the enantiomorphs. c c IWORKS -- c 1 : input position c 3 : symmetry-related position, with OR/EN c include 'main.cmn' c do 1100, m = 1, 16 if (m.gt.8 .and. .not.enan) go to 1100 if (m .gt. 8) then jj = m - 8 else jj = m endif if (isglor(1, jj) .lt. 0) go to 1100 c do 1000, i = 1, nsym do 1200, j = 1, 3 iworks (j, 3) = its (j, i) + isglor (j, jj) do 1300, k = 1, 3 if (m .le. 8) then iworks (j, 3) = iworks (j, 3) . + ifs (j, k, i) * (iworks (k, 1) - 1) else iworks (j, 3) = iworks (j, 3) . + ifs (j, k, i) * (1 - iworks (k, 1)) endif 1300 continue iworks (j, 3) = mod (iworks (j, 3), ncell(j)) if (iworks (j, 3) .lt. 0) iworks (j, 3) . = ncell (j) + iworks (j, 3) iworks (j, 3) = iworks (j, 3) + 1 1200 continue c do 1400, j = 1, nsglpk do 1500, k = 1, 3 if (fixxyz(k)) go to 1500 id = mod(abs (isglpk(k, j) - iworks(k, 3)), ncell(k)) if (id .gt. ntol(k)) id = ncell(k) - id if (id .gt. ntol(k)) go to 1400 1500 continue c i4 = j i3 = 1 i2 = m if (m .gt. 8) then i3 = 2 i2 = m - 8 endif i1 = i return c 1400 continue 1000 continue 1100 continue c i1 = 0 i2 = 0 i3 = 0 i4 = 0 c return end c c----------------------------------------------------------------------- c subroutine dblsho c c This routine checks all the unique self-vectors and cross-vectors c of the input sites. c c SUBROUTINES CALLED -- NUPAGE, SGLSHO, XVECTR c c IWORKS -- c 1 : single site position (i, j, k) c 2 : a single-site vector for (i, j, k) c 3 : a known single-site vector for (i, j, k) c JWORKS -- c 1 - jwork : the single-site vectors and their redundancies c include 'main.cmn' c call nupage c do 1000, i = 1, ndblsh do 1100, j = 1, 3 isglsh (j, 1) = idblsh (j, i) isglsh (j, 2) = idblsh (j+3, i) 1100 continue nsglsh = 2 call sglsho c do 1200, j = 1, 3 iworks (j, 1) = idblsh (j, i) iworks (j, 2) = idblsh (j+3, i) 1200 continue c call xvectr (ival, minx, 0) c write(lprt, 20) do 1300, j = 1, jwork if (map) then write (lprt, 30) j, (jworks(k, j), k=1, 3), . (float(jworks(k, j)-1)/float(ncell(k)), k=1, 3), . abs(jworks(4, j)) - 2000 else write (lprt, 30) j, (jworks(k, j), k=1, 3), . (float(jworks(k, j)-1)/float(ncell(k)), k=1, 3) endif 1300 continue if (map) write (lprt, 40) ival - 2000 1000 continue c 20 format (/, 5x, 'Cross vectors between the two sites ', //, . 9x, 'No.', 3x, ' I ', ' J ', ' K ', 6x, 'u', 6x, 'v', . 6x, 'w', 6x, 'Height') 30 format (10x, i2, 3x, 3i4, 3x, 3f7.3, 3x, i5) 40 format (/, 10x, 'Cross-vector search map value ', i5) c return end c c----------------------------------------------------------------------- c subroutine double c c This routine carries out a two-atom search based on the list of c single-site solutions, by applying all possible origin and c enantiomorph shifts to the single-site positions and checking the c resulting cross-vector heights. c c SUBROUTINES CALLED : SPGASU, XVECTR, SAVDBL, DBLREF, OUTDBL c c IWORKS -- c 1 : atom position #1 c 3 : atom position #2 related by OR/EN c 2 : atom position #2 by cell translation (polar axis) c c JWORKS -- c 1 - : cross vectors between the two atom positions c 1 : two-atom solution site #1 c 2 : two-atom solution site #2 c include 'main.cmn' c dimension ixyz (3) c write (log, 10) ndbl = 0 jdblct = idblct - 2000 if (jdblct .gt. 0) then jdblct = jdblct / 3 else if (jdblct .lt. 0) then jdblct = max (-500, jdblct * 2) else jdblct = -200 endif jdblct = jdblct + 2000 c do 1000, i = 1, 3 ixyz (i) = 1 if (fixxyz(i)) ixyz (i) = ncell (i) 1000 continue c c Loop 1 over the saved single site positions c do 100, i = 1, maxsgl if (isgxyz(4, i) .lt. 0) go to 500 c if (mod(i, 10) .eq. 0) write (log, 20) i do 1100, j = 1, 3 iworks (j, 1) = isgxyz (j, i) 1100 continue c c Loop 2 over the saved single site positions c do 200, j = i, maxsgl if (isgxyz(4, j) .lt. 0) go to 100 c c Loop over the possible OR/EN combinations c do 300, k = 1, 16 if (j.eq.i .and. k.eq.1 .and. .not.fixxyz(1) . .and. .not.fixxyz(2) .and. .not.fixxyz(3)) go to 300 if (k.gt.8 .and. .not.enan) go to 300 if (k .gt. 8) then kk = k - 8 else kk = k endif if (isglor(1, kk) .lt. 0) go to 300 c do 1200, m = 1, 3 if (k .gt. 8) then iworks (m, 3) = -(isgxyz(m, j) - 1) + . isglor(m, kk) + 1 else iworks (m, 3) = isgxyz(m, j) + isglor(m, kk) endif 1200 continue c c Loop over the possible polar axis directions c do 1300, L = 1, ixyz (1) iworks (1, 2) = iworks (1, 3) + L - 1 do 1400, m = 1, ixyz (2) iworks (2, 2) = iworks (2, 3) + m - 1 do 400, n = 1, ixyz (3) iworks (3, 2) = iworks (3, 3) + n - 1 c call spgasu(iasu, 1) if (iasu .eq. 1) go to 400 c call xvectr (mn, minx, jdblct) if (mn .lt. jdblct) go to 400 c do 1500, ii = 1, 3 jworks (ii, 1) = iworks (ii, 1) jworks (ii, 2) = iworks (ii, 2) 1500 continue jworks (4, 1) = isgxyz (4, i) jworks (4, 2) = isgxyz (4, j) call savdbl (mn) c 400 continue 1400 continue 1300 continue c 300 continue 200 continue 100 continue c 500 continue c call dblref call outdbl (3) c 10 format (' DOUBLE> ', 'Two-atom search based on saved ', . 'single-site positions') 20 format (' DOUBLE> ', 'Working on single-site # ', i3) c return end c c----------------------------------------------------------------------- c subroutine endgam c include 'main.cmn' c character*80 record, rekord c call nupage rewind (unit=lscrch) write (lprt, 10) 100 read (lscrch, 20, end=99) record nch = -80 rekord = record call parser (record, rekord, nch) write (lprt, 30) rekord (iw(1, 1):iw(2, nw)) go to 100 c 99 continue c write (lprt, 40) write (lprt, 50) mapdim write (lprt, 60) maxsgl, maxdbl, maxsln write (lprt, 70) maxcrs, maxsho write (lprt, 80) maxlax, maxncr, maxvec c close (unit=lscrch, status='DELETE') close (unit=lprt) c close (unit=lin) c close (unit=log) c call exit c 10 format (5x, 'Listing of input commands ---', //) 20 format (a) 30 format (10x, a) 40 format (//, 5x, 'Program storage limits ---', //) 50 format (10x, 'MAPDIM=', i7) 60 format (10x, 'MAXSGL=', i4, ', MAXDBL=', i4, ', MAXSLN=', i3) 70 format (10x, 'MAXCRS=', i4, ', MAXSHO=', i4) 80 format (10x, 'MAXLAX=', i4, ', MAXNCR=', i4, ', MAXVEC=', i4) c end c c----------------------------------------------------------------------- c subroutine euler1 (th1, th2, th3, rotmat) c c c NOTE : This is the ZY'Z" definition of Euler angles. That is, c rotation th1 is around Z, th2 is around new Y, th3 is around c new Z. The formula is from Paula Fitzgerald's MERLOT package. c c NOTE !! This matrix should be applied to positions, the transpose c should be applied for axes. c include 'main.cmn' c dimension rotmat (3, 3) c dtr = acos(-1.0) / 180.0 c c Check to see whether we are going from angles to a matrix. c if (th1 .lt. -999.0) go to 100 c tth1 = th1 * dtr tth2 = th2 * dtr tth3 = th3 * dtr c c1 = cos(tth1) s1 = sin(tth1) c2 = cos(tth2) s2 = sin(tth2) c3 = cos(tth3) s3 = sin(tth3) c rotmat (1, 1) = c1*c2*c3 - s1*s3 rotmat (2, 1) = s1*c2*c3 + c1*s3 rotmat (3, 1) = -s2*c3 rotmat (1, 2) = -c1*c2*s3 - s1*c3 rotmat (2, 2) = -s1*c2*s3 + c1*c3 rotmat (3, 2) = s2*s3 rotmat (1, 3) = c1*s2 rotmat (2, 3) = s1*s2 rotmat (3, 3) = c2 c return c 100 continue c c2 = rotmat (3, 3) if (c2 .lt. -1.0) c2 = -1.0 if (c2 .gt. 1.0) c2 = 1.0 th2 = acos (c2) s2 = sin (th2) th2 = th2 / dtr c c case 1 : th2 .ne. 0 or pi. arcsin (0.005) = 0.286 deg. c if (abs(s2) .gt. 0.005) then s1 = rotmat (2, 3) / s2 c1 = rotmat (1, 3) / s2 th1 = atan2 (s1, c1) / dtr th1 = mod (th1 + 3600.0, 360.0) s3 = rotmat (3, 2) / s2 c3 = - rotmat (3, 1) / s2 th3 = atan2 (s3, c3) / dtr th3 = mod (th3 + 3600.0, 360.0) return endif c c case 2 : th2 .eq. 0 or pi c when th2 = 0, only th1 + th3 can be determined c when th2 = pi, only th1 - th3 can be determined. c in either case, on return from the routine, th1 will carry c th1 +/- th3, and th3 will be set to 0. c s2 = 0.0 c2 = sign (1.0, c2) th2 = 90.0 - c2 * 90.0 c th1 = atan2 (c2*rotmat (2, 1), c2*rotmat (1, 1)) / dtr th1 = mod (3600.0 + th1, 360.0) th3 = 0 c return end c c----------------------------------------------------------------------- c subroutine euler2 (th1, th2, th3, rotmat) c c NOTE : This is the ZX'Z" definition of Euler angles. That is, c rotation th1 is around Z, th2 is around new X, th3 is around c new Z. See "The Molecular Replacement Method", ed. M.G. Rossmann, c p48 for details. c c NOTE !! This matrix is the TRANSPOSE of the one in Table 1 on c page 48 of the book. This is because the matrix will be applied c to positions rather than axes in the program. c include 'main.cmn' c dimension rotmat (3, 3) c dtr = acos(-1.0) / 180.0 c c Check to see whether we are going from angles to a matrix. c if (th1 .lt. -999.0) go to 100 c tth1 = th1 * dtr tth2 = th2 * dtr tth3 = th3 * dtr c c1 = cos(tth1) s1 = sin(tth1) c2 = cos(tth2) s2 = sin(tth2) c3 = cos(tth3) s3 = sin(tth3) c rotmat (1, 1) = -s1*c2*s3 + c1*c3 rotmat (2, 1) = c1*c2*s3 + s1*c3 rotmat (3, 1) = s2*s3 rotmat (1, 2) = -s1*c2*c3 - c1*s3 rotmat (2, 2) = c1*c2*c3 - s1*s3 rotmat (3, 2) = s2*c3 rotmat (1, 3) = s1*s2 rotmat (2, 3) = -c1*s2 rotmat (3, 3) = c2 c return c 100 continue c c2 = rotmat (3, 3) if (c2 .lt. -1.0) c2 = -1.0 if (c2 .gt. 1.0) c2 = 1.0 th2 = acos (c2) s2 = sin (th2) th2 = th2 / dtr c c case 1 : th2 .ne. 0 or pi. arcsin (0.005) = 0.286 deg. c if (abs(s2) .gt. 0.005) then s1 = rotmat (1, 3) / s2 c1 = - rotmat (2, 3) / s2 th1 = atan2 (s1, c1) / dtr th1 = mod (th1 + 3600.0, 360.0) s3 = rotmat (3, 1) / s2 c3 = rotmat (3, 2) / s2 th3 = atan2 (s3, c3) / dtr th3 = mod (th3 + 3600.0, 360.0) return endif c c case 2 : th2 .eq. 0 or pi c when th2 = 0, only th1 + th3 can be determined c when th2 = pi, only th1 - th3 can be determined. c in either case, on return from the routine, th1 will carry c th1 +/- th3, and th3 will be set to 0. c s2 = 0.0 c2 = sign (1.0, c2) th2 = 90.0 - c2 * 90.0 c th1 = atan2 (rotmat (2, 1), rotmat (1, 1)) / dtr th1 = mod (3600.0 + th1, 360.0) th3 = 0 c return end c c----------------------------------------------------------------------- c subroutine fillup c c This routine fills up the self- and cross-vector values in the c array ISLNHT. Also, if SOL is false, the program will dump the c array as well as all the self and cross-vectors to the output. c c SUBROUTINES CALLED : SVECTR, XVECTR, DBLSGL, SVECTR, XVECTR c c IWORKS -- c 1 : one site c 2 : another site c include 'main.cmn' c do 1000, i = 1, nsoln c do 1100, j = 1, 3 iworks (j, 1) = isoln (j, i) 1100 continue c do 1200, j = i, nsoln do 1300, k = 1, 3 iworks (k, 2) = isoln (k, j) 1300 continue c if (j .eq. i) then call svectr (ival, minhk, 0) else call xvectr (ival, minx, 0) endif c islnht (i, j) = ival islnht (j, i) = ival c 1200 continue 1000 continue c if (sol) return c call nupage write (lprt, 10) do 1400, i = 1, nsoln do 1500, j = 1, 3 iworks (j, 1) = isoln (j, i) 1500 continue call dblsgl (i1, i2, i3, i4) write (lprt, 20) i, (float(isoln(j, i) - 1) / float(ncell(j)), . j=1, 3), i4, i1, i2, i3, . (islnht (j, i) - 2000, j=1, i) 1400 continue c c Now first list all the single-site vectors for the sites c call nupage c write (lprt, 30) c do 1600, i = 1, nsoln c do 1700, j = 1, 3 iworks (j, 1) = isoln (j, i) 1700 continue c call svectr (ival, minhk, 0) c write (lprt, 40) i, . (float (iworks (j, 1) - 1)/float (ncell(j)), j=1, 3) c do 1800, j = 1, jwork if (map) then call intpol (-j, ii) ii = float(ii - 2000) / . float(jworks(4, j) - 2000) + 0.5 write (lprt, 50) j, (jworks(k, j), k=1, 3), . (float(jworks(k, j)-1)/float(ncell(k)), k=1, 3), . ii, jworks(4, j) - 2000 else write (lprt, 50) j, (jworks(k, j), k=1, 3), . (float(jworks(k, j)-1)/float(ncell(k)), k=1, 3), . jworks(4, j) endif 1800 continue if (map) write (lprt, 60) ival - 2000 1600 continue c c Now list all the cross vectors c call nupage c write (lprt, 70) c do 1900, i = 1, nsoln do 2000, j = i + 1, nsoln c do 2100, k = 1, 3 iworks (k, 1) = isoln (k, i) iworks (k, 2) = isoln (k, j) 2100 continue c call xvectr (ival, minx, 0) c write (lprt, 80) i, j do 2200, jj = 1, jwork if (map) then write (lprt, 50) jj, (jworks(k, jj), k=1, 3), . (float(jworks(k, jj)-1)/float(ncell(k)), k=1, 3), . jworks(4, jj) - 2000 else write (lprt, 50) jj, (jworks(k, jj), k=1, 3), . (float(jworks(k, jj)-1)/float(ncell(k)), k=1, 3) endif 2200 continue if (map) write (lprt, 90) ival - 2000 2000 continue 1900 continue c 10 format (5x, 'Listing of self- and cross-vector heights for the ', . 'input solution --', //) 20 format (5x, i5, 3f7.3, 2x, 4i3, 5x, 15i5) 30 format (5x, 'Listing of self-vectors for all the sites in the ', . 'input solution --', //) 40 format (//, 10x, 'Self-vectors for site #', i2, ', with ', . 'position : ', 3f7.3, //) 50 format (15x, i2, 3x, 3i4, 3x, 3f7.3, 3x, i5, 3x, i5) 60 format (/, 15x, 'Single-site solution map value ', i5) 70 format (5x, 'Listing of cross-vectors between all the sites ', . 'in the input solution --', //) 80 format (//, 10x, 'Cross-vectors between site #', i2, . ' and site #', i2, ' --', //) 90 format (/, 15x, 'Cross-vector search map value ', i5) c return end c c----------------------------------------------------------------------- c subroutine fsfinp (sum1, sum2, rmin, rmax) c c This routine reads in the Patterson map in the FSFOR format c include 'main.cmn' c dimension record (400) dimension mm(3) c c Allow 400 points across c if (nxyz(mapsec(1)) .gt. 400) then write (log, 10) nxyz(mapsec(1)), 400 call endgam endif c nxyz (1) = min (ncell (1), nxyz (1)) nxyz (2) = min (ncell (2), nxyz (2)) nxyz (3) = min (ncell (3), nxyz (3)) c c read the map in, round the values with 1000 as scale factor. c sum1 = 0.0 sum2 = 0.0 rmin = 9999.0 rmax = -999.0 c do 1000, i = 1, nxyz (mapsec(3)) mm (mapsec(3)) = i do 1100, k = 1, ncell (mapsec(2)) mm (mapsec(2)) = k read (lmap, err=98, end=99) . (record(j), j=1, nxyz (mapsec(1))) if (k .gt. nxyz(mapsec(2))) go to 1100 do 1200, j = 1, nxyz (mapsec(1)) mm (mapsec(1)) = j rmin = min (rmin, record(j)) rmax = max (rmax, record(j)) sum1 = sum1 + record(j) sum2 = sum2 + record(j) * record(j) kk = mm (1) + (mm (2) - 1) * nxyz (1) + . (mm (3) - 1) * nxyz (1) * nxyz (2) ipatmp (kk) = record (j) * 1000.0 1200 continue 1100 continue 1000 continue rewind (unit=lmap) c k = nxyz (1) * nxyz (2) * nxyz (3) sum1 = sum1 / float (k) sum2 = sum2 / float (k) - sum1 * sum1 if (sum2 .lt. 0.) sum2 = 0.0 sum2 = sqrt (sum2) c close (unit=lmap) return c c Error handling c 98 write (log, 30) call endgam c 99 write (log, 20) call endgam c 10 format (' FSFINP> *Fatal* ', 'Too many numbers per line > ', i3, . '. Current maximum is ', i3, /, ' PROGRAM ABORT') 20 format (' FSFINP> *Fatal* ', 'The Map file does not contain ', . 'enough records', /, ' PROGRAM ABORT') 30 format (' FSFINP> *Fatal* ', 'Error reading the map file', /, . ' PROGRAM ABORT') c return end c c----------------------------------------------------------------------- c subroutine gtcell c include 'main.cmn' c do 1000, i = 4, 6 if (abs(cell(i, 1)) .le. 1.0) cell(i, 1) = . acos (cell(i, 1)) / dtr 1000 continue c if (cell(1, 1) .lt. 1.0) then write (log, 10) call endgam endif c do 1200, i = 1, 3 cell (i, 2) = cell (i, 1) ** 2 cell (i + 3, 2) = cos (cell(i + 3, 1) * dtr) if (abs(cell(i+3, 2)) .lt. 0.000001) cell (i + 3, 2) = 0.0 1200 continue c cell (4, 2) = 2.0 * cell(2, 1) * cell(3, 1) * cell (4, 2) cell (5, 2) = 2.0 * cell(1, 1) * cell(3, 1) * cell (5, 2) cell (6, 2) = 2.0 * cell(1, 1) * cell(2, 1) * cell (6, 2) c 10 format (' GTCELL> *Fatal* ', 'Cell parameters must be input', /, . ' PROGRAM ABORT') c return end c c----------------------------------------------------------------------- c subroutine gthead (iopt) c include 'main.cmn' c if (maptyp .eq. 'FSFOR') then call hdrfsf (iopt) else if (maptyp .eq. 'HASSP') then c else if (maptyp(1:4) .eq. 'CCP4') then c else if (maptyp(1:4) .eq. '3D01') then c else if (maptyp .eq. 'XPLOR') then call hdrxpl (iopt) else write (log, 10) maptyp call endgam endif c 10 format (' GTHEAD> *Fatal* ', 'Unrecognized MAPType ', a, /, . ' GTHEAD> Program Abort') c return end c c----------------------------------------------------------------------- c subroutine gtmatx (iopt) c include 'main.cmn' c character*1 angle dimension rotmat (3, 3) c if (iopt .eq. 1) then th1 = rlcaxs (1, nlcaxs) th2 = rlcaxs (2, nlcaxs) th3 = rlcaxs (3, nlcaxs) angle = tlcaxs (nlcaxs) else if (iopt .eq. 2) then th1 = alcrot (1) th2 = alcrot (2) th3 = alcrot (3) angle = tlcrot (1:1) if (angle.ne.'E' .and. angle.ne.'P') then write (log, 10) angle return endif endif c if (angle .eq. 'E') then if (euler .eq. 'ZYZ') then call euler1 (th1, th2, th3, rotmat) else if (euler .eq. 'ZXZ') then call euler2 (th1, th2, th3, rotmat) else write (log, 10) 'Euler', euler if (iopt .eq. 1) nlcaxs = nlcaxs - 1 return endif else if (angle .eq. 'P') then if (polar .eq. 'XZK') then call polar1 (th1, th2, th3, rotmat) else if (polar .eq. 'XYK') then call polar2 (th1, th2, th3, rotmat) else write (log, 10) 'Polar', polar if (iopt .eq. 1) nlcaxs = nlcaxs - 1 return endif else if (angle .eq. 'V') then if (abs(th1).gt.0.01 .or. abs(th2).gt.0.01) then phi = atan2 (th2, th1) / dtr else phi = 0.0 endif psi = th1 ** 2 + th2 ** 2 + th3 ** 2 psi = acos (th3 / sqrt (psi)) / dtr th1 = phi th2 = psi th3 = 360.0 / float (mlcaxs (nlcaxs)) call polar1 (th1, th2, th3, rotmat) else if (angle .eq. 'D') then if (abs(th1).gt.0.01 .or. abs(th2).gt.0.01) then phi = atan2 (th2, th1) / dtr else phi = 0.0 endif psi = acos (th3) / dtr th1 = phi th2 = psi th3 = 360.0 / float (mlcaxs (nlcaxs)) call polar1 (th1, th2, th3, rotmat) else write (log, 20) angle nlcaxs = nlcaxs - 1 return endif c if (angle .eq. 'E') then if (iopt .eq. 2) tlcrot (2:4) = euler else if (angle .eq. 'P') then if (iopt .eq. 2) tlcrot (2:4) = polar endif c if (iopt .eq. 1) then do 1000, i = 1, 3 do 1100, j = 1, 3 rlcrot (i, j, nlcaxs) = rotmat (i, j) 1100 continue 1000 continue else if (iopt .eq. 2) then do 1200, i = 1, 3 do 1300, j = 1, 3 ematrx (i, j) = rotmat (i, j) 1300 continue 1200 continue endif c 10 format (' GTMATX> *Warning* ', 'Unknown ', a, ' Angle Type <', a, . '> -- Axis IGNORED') 20 format (' GTMATX> *Warning* ', 'Unknown Axis Specification <', a, . '>, Only E, P, V, or D supported -- IGNORED') c return end c c----------------------------------------------------------------------- c subroutine gtnabr (i, j, k, nabor) c c This routine finds the neighboring grid points for position c (i, j, k). No boundary condition is assumed, unless the grid c point is lying right at the cell edge. c c IWORKS -- c 150 : the input position c 1 - nabor : the neighbors c include 'main.cmn' c nbox = 3 nbox2 = (nbox + 1) / 2 nabor = 0 c do 1000, i1 = 1, nbox do 1100, j1 = 1, nbox do 100, k1 = 1, nbox if (i1.eq.nbox2 .and. j1.eq.nbox2 . .and. k1.eq.nbox2) go to 100 c iworks (1, 150) = i + (i1 - nbox2) * 2 iworks (2, 150) = j + (j1 - nbox2) * 2 iworks (3, 150) = k + (k1 - nbox2) * 2 c do 1200, m = 1, 3 if (iworks(m, 150) .eq. 0) then if (nxyz(m) .eq. ncell(m)) then iworks (m, 150) = ncell(m) else go to 100 endif else if (iworks(m, 150) .lt. 0) then if (nxyz(m) .eq. ncell(m)) then iworks (m, 150) = ncell(m) + iworks (m, 150) else go to 100 endif else if (iworks(m, 150) .gt. nxyz(m)) then if (nxyz(m) .eq. ncell(m)) then iworks (m, 150) = mod (iworks(m, 150), nxyz(m)) else go to 100 endif endif 1200 continue c if (nabor .eq. 149) then write (log, 10) nabor = 150 return endif c nabor = nabor + 1 do 1300, m = 1, 3 iworks (m, nabor) = iworks (m, 150) 1300 continue c 100 continue 1100 continue 1000 continue c 10 format (' GTNABR> *Warning* ', 'Only the first 150 neighbors ', . 'used in the peak search') c return end c c----------------------------------------------------------------------- c subroutine gtomdm c include 'main.cmn' c dimension omdm(3, 3), rcell(6) c ca = cos (cell(4, 1) * dtr) if (abs(ca) .lt. 0.001) ca = 0.0 cb = cos (cell(5, 1) * dtr) if (abs(cb) .lt. 0.001) cb = 0.0 cg = cos (cell(6, 1) * dtr) if (abs(cg) .lt. 0.001) cg = 0.0 sa = sin (cell(4, 1) * dtr) sb = sin (cell(5, 1) * dtr) sg = sin (cell(6, 1) * dtr) c call recell (rcell) c c BYBCX (b along Y, and b x c along X). Good for b-unique space c groups. See "The molecular replacement method", ed. M.G. Rossmann, c p48 for details. c if (ordor .eq. 'BYBCX') then co = - rcell (5) so = sqrt (1.0 - co * co) omdm (1, 1) = cell (1, 1) * sg * so omdm (1, 2) = 0.0 omdm (1, 3) = 0.0 omdm (2, 1) = cell (1, 1) * cg omdm (2, 2) = cell (2, 1) omdm (2, 3) = cell (3, 1) * ca omdm (3, 1) = cell (1, 1) * sg * co omdm (3, 2) = 0.0 omdm (3, 3) = cell (3, 1) * sa c c CZBCX (c along Z, and b x c along X). Good for c-unique space c groups. The matrix is obtained from MERLOT by P.M.D. Fitzgerald. c else if (ordor .eq. 'CZBCX') then cgs = rcell (6) sgs = sqrt (1.0 - cgs * cgs) omdm (1, 1) = cell (1, 1) * sb * sgs omdm (1, 2) = 0.0 omdm (1, 3) = 0.0 omdm (2, 1) = - cell (1, 1) * sb * cgs omdm (2, 2) = cell (2, 1) * sa omdm (2, 3) = 0.0 omdm (3, 1) = cell (1, 1) * cb omdm (3, 2) = cell (2, 1) * ca omdm (3, 3) = cell (3, 1) c c AXABZ (a along X, and a x b along Z). This is the matrix used in c the program SAMAIN in FRODO, and X-Plor, and NCODE=1 in ALMN. c else if (ordor .eq. 'AXABZ') then cas = rcell (4) sas = sqrt (1.0 - cas * cas) omdm (1, 1) = cell (1, 1) omdm (1, 2) = cell (2, 1) * cg omdm (1, 3) = cell (3, 1) * cb omdm (2, 1) = 0.0 omdm (2, 2) = cell (2, 1) * sg omdm (2, 3) = - sb * cell (3, 1) * cas omdm (3, 1) = 0.0 omdm (3, 2) = 0.0 omdm (3, 3) = cell (3, 1) * sb * sas c c Y3XYC (111 along Y, and C in XY plane). This should be used only for c rhombohedral space groups in the rhombohedral setting. Courtesy of c Dr. Murthy. c else if (ordor .eq. 'Y3XYC') then sina2 = sin (cell (4, 1) * 0.5 * dtr) omdm (1, 1) = cell (1, 1) * sina2 / sqrt (3.0) omdm (1, 2) = omdm (1, 1) omdm (1, 3) = - 2.0 * omdm (1, 1) omdm (2, 1) = cell (1, 1) . * sqrt (1.0 - 4.0 * sina2 * sina2 / 3.0) omdm (2, 2) = omdm (2, 1) omdm (2, 3) = omdm (2, 1) omdm (3, 1) = cell (1, 1) * sina2 omdm (3, 2) = - omdm (3, 1) omdm (3, 3) = 0.0 else write (log, 20) ordor call endgam endif c call matinv (omdm, alpha) c 20 format (/, ' GTOMDM> *Fatal* ', 'Unknown orthogonalization code ', . a, /, ' PROGRAM ABORT') c return end c c----------------------------------------------------------------------- c subroutine gtsgle (isame) c c This routine determines whether an input position shares the same c self-vectors with one of the sites saved in array ISGXYZ. ISAME c = 0 if there is no match, otherwise ISAME will be the index of the c single-site position. c c IWORKS -- c 1 : reserved from PATDBL c 2 : the position to be checked c 3 : space-group symmetry-related positions of 2, with OR/EN c include 'main.cmn' c do 1000, i = 1, nsym do 100, m = 1, 16 if (m.gt.8 .and. .not.enan) go to 100 if (m .gt. 8) then jj = m - 8 else jj = m endif if (isglor(1, jj) .lt. 0) go to 100 c do 1100, j = 1, 3 iworks (j, 3) = its (j, i) + isglor (j, jj) do 1200, k = 1, 3 if (m .le. 8) then iworks (j, 3) = iworks (j, 3) . + ifs (j, k, i) * (iworks (k, 2) - 1) else iworks (j, 3) = iworks (j, 3) . + ifs (j, k, i) * (1 - iworks (k, 2)) endif 1200 continue iworks (j, 3) = mod (iworks (j, 3), ncell(j)) if (iworks (j, 3) .lt. 0) iworks (j, 3) . = ncell (j) + iworks (j, 3) iworks (j, 3) = iworks (j, 3) + 1 1100 continue c do 1300, j = 1, maxsgl if (isgxyz(4, j) .lt. 0) go to 100 do 1400, k = 1, 3 if (.not.fixxyz(k) .and. . isgxyz(k, j) .ne. iworks(k, 3)) go to 1300 1400 continue isame = j return 1300 continue c 100 continue 1000 continue c isame = 0 c return end c c----------------------------------------------------------------------- c subroutine hdrfsf (iopt) c include 'main.cmn' c dimension celdim(6), itle(20), mnx(6), nuvw(3) dimension trans(3, 96) c if (iopt .eq. 1) go to 100 c open (unit=lmap, status='OLD', form='UNFORMATTED', file=mapfil, . err=99) c c The cell parameters from FSFOR will be ignored c read (lmap, err=98, end=98) itle, noset, celdim, nsym, ncent, . latt, npic c c Read in the FSFOR packed symmetry cards c do 1000, i = 1, nsym read (lmap, err=98, end=98) (trans (j, i), . (ifs (k, j, i), k=1, 2), j=1, 3) 1000 continue call putsym (trans) c c Read map grid information c read (lmap, err=98, end=98) nuvw, scale, norn, mnx, rspmx, nbit c if (norn .eq. 0) then mapsec (3) = 2 mapsec (1) = 1 mapsec (2) = 3 else if (norn .eq. 1) then mapsec (3) = 1 mapsec (1) = 2 mapsec (2) = 3 else if (norn .eq. 2) then mapsec (3) = 3 mapsec (1) = 1 mapsec (2) = 2 endif c ncell (mapsec(3)) = nuvw (3) ncell (mapsec(2)) = nuvw (2) ncell (mapsec(1)) = nuvw (1) c return c c Rewind map and skip header c 100 rewind lmap read (lmap) do 1200, i = 1, nsym read (lmap) 1200 continue read (lmap) return c 98 write (log, 20) call endgam c 99 write (log, 10) mapfil call endgam c 10 format (' HDRFSF> *Fatal* ', 'Error opening map file ', a, . /, ' HDRFSF> Program Abort') 20 format (' HDRFSF> *Fatal* ', 'Error reading from map file', /, . ' HDRFSF> ', 'Program Abort') c return end c c----------------------------------------------------------------------- c subroutine hdrxpl (iopt) c include 'main.cmn' c character*3 secdir integer mxyz(3) double precision dcell(6) c if (iopt .eq. 1) go to 100 c if (qform) then open (unit=lmap, status='OLD', form='FORMATTED', file=mapfil, . err=99) else open (unit=lmap, status='OLD', form='UNFORMATTED', . file=mapfil, err=99) endif c if (qform) then read (lmap, 40, err=98, end=98) read (lmap, 30, err=98, end=98) ntitle do 1000, i = 1, ntitle read (lmap, 40, err=98, end=98) 1000 continue else read (lmap, err=98, end=98) ntitle endif c if (qform) then read (lmap, 50, err=98, end=98) (ncell (i), mxyz (i), . mxyz (i), i=1, 3) else read (lmap, err=98, end=98) (ncell (i), mxyz (i), . mxyz (i), i=1, 3) endif c if (qform) then read (lmap, 60, err=98, end=98) (cell (i, 1), i=1, 6) else read (lmap, err=98, end=98) (dcell (i), i=1, 6) do 1200, i = 1, 6 cell (i, 1) = dcell (i) 1200 continue endif c if (qform) then read (lmap, 40, err=98, end=98) secdir else read (lmap, err=98, end=98) secdir endif c if (secdir(1:1) .eq. 'Z') mapsec (3) = 3 if (secdir(1:1) .eq. 'Y') mapsec (3) = 2 if (secdir(1:1) .eq. 'X') mapsec (3) = 1 if (secdir(2:2) .eq. 'Z') mapsec (2) = 3 if (secdir(2:2) .eq. 'Y') mapsec (2) = 2 if (secdir(2:2) .eq. 'X') mapsec (2) = 1 if (secdir(3:3) .eq. 'Z') mapsec (1) = 3 if (secdir(3:3) .eq. 'Y') mapsec (1) = 2 if (secdir(3:3) .eq. 'X') mapsec (1) = 1 c return c 100 rewind lmap c if (qform) then read (lmap, 40, err=98, end=98) read (lmap, 30, err=98, end=98) ntitle do 1100, i = 1, ntitle read (lmap, 40, err=98, end=98) 1100 continue read (lmap, 40, err=98, end=98) read (lmap, 40, err=98, end=98) read (lmap, 40, err=98, end=98) else read (lmap, err=98, end=98) read (lmap, err=98, end=98) read (lmap, err=98, end=98) read (lmap, err=98, end=98) endif return c 98 write (log, 20) call endgam c 99 write (log, 10) mapfil call endgam c 10 format (' HDRXPL> *Fatal* ', 'Error opening map file ', a, . /, ' HDRXPL> Program Abort') 20 format (' HDRXPL> *Fatal* ', 'Error reading from map file', /, . ' HDRXPL> ', 'Program Abort') 30 format (i8) 40 format (a) 50 format (9i8) 60 format (6e12.5) c return end c c----------------------------------------------------------------------- c subroutine hspinp (sum1, sum2, rmin, rmax) c c This routine reads in the Patterson map in the HASSP format c (modified for general sectioning) c include 'main.cmn' c dimension record (400) dimension mm(3) c c Allow 400 points across c if (nxyz(mapsec(1)) .gt. 400) then write (log, 10) nxyz(mapsec(1)), 400 call endgam endif c open (unit=lmap, file=mapfil, status='OLD', form='UNFORMATTED', . err=97) c c read the map in, round the values with 1000 as scale factor. c sum1 = 0.0 sum2 = 0.0 rmin = 9999.0 rmax = -999.0 c do 1000, i = 1, nxyz (mapsec(3)) mm (mapsec(3)) = i do 1100, k = 1, nxyz (mapsec(2)) mm (mapsec(2)) = k read (lmap, err=98, end=99) . (record(j), j=1, nxyz (mapsec(1))) do 1200, j = 1, nxyz (mapsec(1)) mm (mapsec(1)) = j rmin = min (rmin, record(j)) rmax = max (rmax, record(j)) sum1 = sum1 + record(j) sum2 = sum2 + record(j) * record(j) kk = mm (1) + (mm (2) - 1) * nxyz (1) + . (mm (3) - 1) * nxyz (1) * nxyz (2) ipatmp (kk) = record (j) * 1000.0 1200 continue 1100 continue 1000 continue rewind (unit=lmap) c k = nxyz (1) * nxyz (2) * nxyz (3) sum1 = sum1 / float (k) sum2 = sum2 / float (k) - sum1 * sum1 if (sum2 .lt. 0.) sum2 = 0.0 sum2 = sqrt (sum2) c close (unit=lmap) return c c Error handling c 97 write (log, 50) mapfil call endgam c 98 write (log, 30) call endgam c 99 write (log, 20) call endgam c 10 format (' HSPINP> *Fatal* ', 'Too many numbers per line > ', i3, . '. Current maximum is ', i3, /, ' PROGRAM ABORT') 20 format (' HSPINP> *Fatal* ', 'The Map file does not contain ', . 'enough records', /, ' PROGRAM ABORT') 30 format (' HSPINP> *Fatal* ', 'Error reading the map file', /, . ' PROGRAM ABORT') 50 format (' HSPINP> *Fatal* ', 'Error opening the map file -> ', a, . /, ' PROGRAM ABORT') c return end c c----------------------------------------------------------------------- c subroutine iiread (word, nch, ivalue) c include 'main.cmn' c character word*(*) character*1 digit(14) data digit/'1', '2', '3', '4', '5', '6', '7', '8', '9', '0', . '+', '-', '.', 'E'/ c idot = 0 ie = 0 do 1000, i = 1, nch if (word(i:i) .eq. '.') idot = i if (word(i:i) .eq. 'E') ie = i idif = 0 do 1100, j = 1, 14 if (word(i:i) .eq. digit(j)) idif = idif + 1 1100 continue if (idif .eq. 0) then write (log, 10) word call exit endif 1000 continue c if (idot.eq.0 .and. ie.eq.0) then ign = 1 ivalue = 0 do 1200, i = 1, nch do 1300, j = 1, 9 if (word(i:i) .eq. digit(j)) ivalue = ivalue . + j * 10 ** (nch - i) 1300 continue if (word (i:i) .eq. digit(12)) ign = -1 1200 continue ivalue = ivalue * ign else call rrread (word, nch, value) ivalue = value + sign (0.5, value) endif c 10 format (' IIREAD> *Fatal* ', 'Unrecognized character in ', . 'number > ', a, /, ' PROGRAM ABORT') c return end c c----------------------------------------------------------------------- c subroutine inispg c c Set up the space group symmetry operators for the 65 c non-centrosymmetric space groups c include 'main.cmn' c c First, the space group names c do 1000, i = 1, 230 spgnam (i) = ' ' 1000 continue c spgnam (1) = 'P1 ' c spgnam (3) = 'P2 ' spgnam (4) = 'P21 ' spgnam (5) = 'C2 ' c spgnam (16) = 'P222 ' spgnam (17) = 'P2221 ' spgnam (18) = 'P21212 ' spgnam (19) = 'P212121' spgnam (20) = 'C2221 ' spgnam (21) = 'C222 ' spgnam (22) = 'F222 ' spgnam (23) = 'I222 ' spgnam (24) = 'I212121' c spgnam (75) = 'P4 ' spgnam (76) = 'P41 ' spgnam (77) = 'P42 ' spgnam (78) = 'P43 ' spgnam (79) = 'I4 ' spgnam (80) = 'I41 ' c spgnam (89) = 'P422 ' spgnam (90) = 'P4212 ' spgnam (91) = 'P4122 ' spgnam (92) = 'P41212 ' spgnam (93) = 'P4222 ' spgnam (94) = 'P42212 ' spgnam (95) = 'P4322 ' spgnam (96) = 'P43212 ' spgnam (97) = 'I422 ' spgnam (98) = 'I4122 ' c spgnam (143) = 'P3 ' spgnam (144) = 'P31 ' spgnam (145) = 'P32 ' spgnam (146) = 'R3 ' c spgnam (149) = 'P312 ' spgnam (150) = 'P321 ' spgnam (151) = 'P3112 ' spgnam (152) = 'P3121 ' spgnam (153) = 'P3212 ' spgnam (154) = 'P3221 ' spgnam (155) = 'R32 ' c spgnam (168) = 'P6 ' spgnam (169) = 'P61 ' spgnam (170) = 'P65 ' spgnam (171) = 'P62 ' spgnam (172) = 'P64 ' spgnam (173) = 'P63 ' c spgnam (177) = 'P622 ' spgnam (178) = 'P6122 ' spgnam (179) = 'P6522 ' spgnam (180) = 'P6222 ' spgnam (181) = 'P6422 ' spgnam (182) = 'P6322 ' c spgnam (195) = 'P23 ' spgnam (196) = 'F23 ' spgnam (197) = 'I23 ' spgnam (198) = 'P213 ' spgnam (199) = 'I213 ' c spgnam (207) = 'P432 ' spgnam (208) = 'P4232 ' spgnam (209) = 'F432 ' spgnam (210) = 'F4132 ' spgnam (211) = 'I432 ' spgnam (212) = 'P4332 ' spgnam (213) = 'P4132 ' spgnam (214) = 'I4132 ' c c Next, the rotation matrices c do 1100, i = 1, 34 do 1200, j = 1, 3 do 1300, k = 1, 3 isprot (j, k, i) = 0 1300 continue 1200 continue 1100 continue c isprot (1, 1, 1) = 1 isprot (2, 2, 1) = 1 isprot (3, 3, 1) = 1 c isprot (1, 1, 2) = -1 isprot (2, 2, 2) = 1 isprot (3, 3, 2) = -1 c isprot (1, 1, 3) = 1 isprot (2, 2, 3) = -1 isprot (3, 3, 3) = -1 c isprot (1, 1, 4) = -1 isprot (2, 2, 4) = -1 isprot (3, 3, 4) = 1 c isprot (1, 2, 5) = -1 isprot (2, 1, 5) = 1 isprot (3, 3, 5) = 1 c isprot (1, 2, 6) = 1 isprot (2, 1, 6) = -1 isprot (3, 3, 6) = 1 c isprot (1, 2, 7) = 1 isprot (2, 1, 7) = 1 isprot (3, 3, 7) = -1 c isprot (1, 2, 8) = -1 isprot (2, 1, 8) = -1 isprot (3, 3, 8) = -1 c isprot (1, 2, 9) = -1 isprot (2, 1, 9) = 1 isprot (2, 2, 9) = -1 isprot (3, 3, 9) = 1 c isprot (1, 1, 10) = -1 isprot (1, 2, 10) = 1 isprot (2, 1, 10) = -1 isprot (3, 3, 10) = 1 c isprot (1, 1, 11) = -1 isprot (1, 2, 11) = 1 isprot (2, 2, 11) = 1 isprot (3, 3, 11) = -1 c isprot (1, 1, 12) = 1 isprot (2, 1, 12) = 1 isprot (2, 2, 12) = -1 isprot (3, 3, 12) = -1 c isprot (1, 1, 13) = 1 isprot (1, 2, 13) = -1 isprot (2, 2, 13) = -1 isprot (3, 3, 13) = -1 c isprot (1, 1, 14) = -1 isprot (2, 1, 14) = -1 isprot (2, 2, 14) = 1 isprot (3, 3, 14) = -1 c isprot (1, 2, 15) = 1 isprot (2, 1, 15) = -1 isprot (2, 2, 15) = 1 isprot (3, 3, 15) = 1 c isprot (1, 1, 16) = 1 isprot (1, 2, 16) = -1 isprot (2, 1, 16) = 1 isprot (3, 3, 16) = 1 c isprot (1, 3, 17) = 1 isprot (2, 1, 17) = 1 isprot (3, 2, 17) = 1 c isprot (1, 3, 18) = 1 isprot (2, 1, 18) = -1 isprot (3, 2, 18) = -1 c isprot (1, 3, 19) = -1 isprot (2, 1, 19) = -1 isprot (3, 2, 19) = 1 c isprot (1, 3, 20) = -1 isprot (2, 1, 20) = 1 isprot (3, 2, 20) = -1 c isprot (1, 2, 21) = 1 isprot (2, 3, 21) = 1 isprot (3, 1, 21) = 1 c isprot (1, 2, 22) = -1 isprot (2, 3, 22) = 1 isprot (3, 1, 22) = -1 c isprot (1, 2, 23) = 1 isprot (2, 3, 23) = -1 isprot (3, 1, 23) = -1 c isprot (1, 2, 24) = -1 isprot (2, 3, 24) = -1 isprot (3, 1, 24) = 1 c isprot (1, 2, 25) = 1 isprot (2, 1, 25) = -1 isprot (3, 3, 25) = 1 c isprot (1, 2, 26) = -1 isprot (2, 1, 26) = 1 isprot (3, 3, 26) = 1 c isprot (1, 1, 27) = 1 isprot (2, 3, 27) = 1 isprot (3, 2, 27) = -1 c isprot (1, 1, 28) = -1 isprot (2, 3, 28) = 1 isprot (3, 2, 28) = 1 c isprot (1, 1, 29) = -1 isprot (2, 3, 29) = -1 isprot (3, 2, 29) = -1 c isprot (1, 1, 30) = 1 isprot (2, 3, 30) = -1 isprot (3, 2, 30) = 1 c isprot (1, 3, 31) = 1 isprot (2, 2, 31) = 1 isprot (3, 1, 31) = -1 c isprot (1, 3, 32) = 1 isprot (2, 2, 32) = -1 isprot (3, 1, 32) = 1 c isprot (1, 3, 33) = -1 isprot (2, 2, 33) = 1 isprot (3, 1, 33) = 1 c isprot (1, 3, 34) = -1 isprot (2, 2, 34) = -1 isprot (3, 1, 34) = -1 c c The translation vectors c do 1400, i = 1, 26 do 1500, j = 1, 3 ispts (j, i) = 0 1500 continue 1400 continue c ispts (2, 2) = 6 c ispts (3, 3) = 6 c ispts (1, 4) = 6 c ispts (1, 5) = 6 ispts (2, 5) = 6 c ispts (1, 6) = 6 ispts (3, 6) = 6 c ispts (2, 7) = 6 ispts (3, 7) = 6 c ispts (3, 8) = 3 c ispts (3, 9) = 9 c ispts (1, 10) = 6 ispts (2, 10) = 6 ispts (3, 10) = 6 c ispts (2, 11) = 6 ispts (3, 11) = 3 c ispts (1, 12) = 6 ispts (3, 12) = 9 c ispts (1, 13) = 6 ispts (2, 13) = 6 ispts (3, 13) = 3 c ispts (1, 14) = 6 ispts (2, 14) = 6 ispts (3, 14) = 9 c ispts (3, 15) = 4 c ispts (3, 16) = 8 c ispts (3, 17) = 10 c ispts (3, 18) = 2 c ispts (1, 19) = 9 ispts (2, 19) = 3 ispts (3, 19) = 9 c ispts (1, 20) = 3 ispts (2, 20) = 3 ispts (3, 20) = 3 c ispts (1, 21) = 3 ispts (2, 21) = 9 ispts (3, 21) = 9 c ispts (1, 22) = 9 ispts (2, 22) = 9 ispts (3, 22) = 3 c ispts (1, 23) = 9 ispts (2, 23) = 3 ispts (3, 23) = 3 c ispts (1, 24) = 9 ispts (2, 24) = 9 ispts (3, 24) = 9 c ispts (1, 25) = 3 ispts (2, 25) = 3 ispts (3, 25) = 9 c ispts (1, 26) = 3 ispts (2, 26) = 9 ispts (3, 26) = 3 c c The space groups c do 1600, i = 1, 230 nsprot (1, i) = 1 nsplat (i) = 1 do 1700, j = 1, 48 nspts (j, i) = 1 1700 continue 1600 continue c c Triclinic c nspsym (1) = 1 c c Monoclinic c do 1800, i = 3, 5 nspsym (i) = 2 nsprot (2, i) = 2 1800 continue nspts (2, 4) = 2 nsplat (5) = 4 c c Orthorhombic c do 1900, i = 16, 24 nspsym (i) = 4 nsprot (2, i) = 4 nsprot (3, i) = 2 nsprot (4, i) = 3 1900 continue nsplat (20) = 4 nsplat (21) = 4 nsplat (22) = 5 nsplat (23) = 6 nsplat (24) = 6 c nspts (2, 17) = 3 nspts (3, 17) = 3 c nspts (3, 18) = 5 nspts (4, 18) = 5 c nspts (2, 19) = 6 nspts (3, 19) = 7 nspts (4, 19) = 5 c nspts (2, 20) = 3 nspts (3, 20) = 3 c nspts (2, 24) = 6 nspts (3, 24) = 7 nspts (4, 24) = 5 c c Tetragonal (4) c do 2000, i = 75, 80 nspsym (i) = 4 nsprot (2, i) = 4 nsprot (3, i) = 5 nsprot (4, i) = 6 2000 continue nsplat (79) = 6 nsplat (80) = 6 c nspts (2, 76) = 3 nspts (3, 76) = 8 nspts (4, 76) = 9 c nspts (3, 77) = 3 nspts (4, 77) = 3 c nspts (2, 78) = 3 nspts (3, 78) = 9 nspts (4, 78) = 8 c nspts (2, 80) = 10 nspts (3, 80) = 11 nspts (4, 80) = 12 c c Tetragonal (422) c do 2100, i = 89, 98 nspsym (i) = 8 nsprot (2, i) = 4 nsprot (3, i) = 5 nsprot (4, i) = 6 nsprot (5, i) = 2 nsprot (6, i) = 3 nsprot (7, i) = 7 nsprot (8, i) = 8 2100 continue nsplat (97) = 6 nsplat (98) = 6 c nspts (3, 90) = 5 nspts (4, 90) = 5 nspts (5, 90) = 5 nspts (6, 90) = 5 c nspts (2, 91) = 3 nspts (3, 91) = 8 nspts (4, 91) = 9 nspts (6, 91) = 3 nspts (7, 91) = 9 nspts (8, 91) = 8 c nspts (2, 92) = 3 nspts (3, 92) = 13 nspts (4, 92) = 14 nspts (5, 92) = 13 nspts (6, 92) = 14 nspts (8, 92) = 3 c nspts (3, 93) = 3 nspts (4, 93) = 3 nspts (7, 93) = 3 nspts (8, 93) = 3 c nspts (3, 94) = 10 nspts (4, 94) = 10 nspts (5, 94) = 10 nspts (6, 94) = 10 c nspts (2, 95) = 3 nspts (3, 95) = 9 nspts (4, 95) = 8 nspts (6, 95) = 3 nspts (7, 95) = 8 nspts (8, 95) = 9 c nspts (2, 96) = 3 nspts (3, 96) = 14 nspts (4, 96) = 13 nspts (5, 96) = 14 nspts (6, 96) = 13 nspts (8, 96) = 3 c nspts (2, 98) = 10 nspts (3, 98) = 11 nspts (4, 98) = 12 nspts (5, 98) = 12 nspts (6, 98) = 11 nspts (7, 98) = 10 c c Trigonal (3) c do 2200, i = 143, 146 nspsym (i) = 3 nsprot (2, i) = 9 nsprot (3, i) = 10 2200 continue nsplat (146) = 7 c nspts (2, 144) = 15 nspts (3, 144) = 16 c nspts (2, 145) = 16 nspts (3, 145) = 15 c c Trigonal (312 and 321) c do 2300, i = 149, 155 nspsym (i) = 6 nsprot (2, i) = 9 nsprot (3, i) = 10 nsprot (4, i) = 8 nsprot (5, i) = 11 nsprot (6, i) = 12 2300 continue nsprot (4, 150) = 7 nsprot (5, 150) = 13 nsprot (6, 150) = 14 nsprot (4, 152) = 7 nsprot (5, 152) = 13 nsprot (6, 152) = 14 nsprot (4, 154) = 7 nsprot (5, 154) = 13 nsprot (6, 154) = 14 nsprot (4, 155) = 7 nsprot (5, 155) = 13 nsprot (6, 155) = 14 nsplat (155) = 7 c nspts (2, 151) = 15 nspts (3, 151) = 16 nspts (4, 151) = 16 nspts (5, 151) = 15 c nspts (2, 152) = 15 nspts (3, 152) = 16 nspts (5, 152) = 16 nspts (6, 152) = 15 c nspts (2, 153) = 16 nspts (3, 153) = 15 nspts (4, 153) = 15 nspts (5, 153) = 16 c nspts (2, 154) = 16 nspts (3, 154) = 15 nspts (5, 154) = 15 nspts (6, 154) = 16 c c Hexagonal (6) c do 2400, i = 168, 173 nspsym (i) = 6 nsprot (2, i) = 9 nsprot (3, i) = 10 nsprot (4, i) = 4 nsprot (5, i) = 15 nsprot (6, i) = 16 2400 continue c nspts (2, 169) = 15 nspts (3, 169) = 16 nspts (4, 169) = 3 nspts (5, 169) = 17 nspts (6, 169) = 18 c nspts (2, 170) = 16 nspts (3, 170) = 15 nspts (4, 170) = 3 nspts (5, 170) = 18 nspts (6, 170) = 17 c nspts (2, 171) = 16 nspts (3, 171) = 15 nspts (5, 171) = 16 nspts (6, 171) = 15 c nspts (2, 172) = 15 nspts (3, 172) = 16 nspts (5, 172) = 15 nspts (6, 172) = 16 c nspts (4, 173) = 3 nspts (5, 173) = 3 nspts (6, 173) = 3 c c Hexagonal (622) c do 2500, i = 177, 182 nspsym (i) = 12 nsprot (2, i) = 9 nsprot (3, i) = 10 nsprot (4, i) = 4 nsprot (5, i) = 15 nsprot (6, i) = 16 nsprot (7, i) = 7 nsprot (8, i) = 13 nsprot (9, i) = 14 nsprot (10, i) = 8 nsprot (11, i) = 11 nsprot (12, i) = 12 2500 continue c nspts (2, 178) = 15 nspts (3, 178) = 16 nspts (4, 178) = 3 nspts (5, 178) = 17 nspts (6, 178) = 18 nspts (7, 178) = 15 nspts (9, 178) = 16 nspts (10, 178) = 17 nspts (11, 178) = 3 nspts (12, 178) = 18 c nspts (2, 179) = 16 nspts (3, 179) = 15 nspts (4, 179) = 3 nspts (5, 179) = 18 nspts (6, 179) = 17 nspts (7, 179) = 16 nspts (9, 179) = 15 nspts (10, 179) = 18 nspts (11, 179) = 3 nspts (12, 179) = 17 c nspts (2, 180) = 16 nspts (3, 180) = 15 nspts (5, 180) = 16 nspts (6, 180) = 15 nspts (7, 180) = 16 nspts (9, 180) = 15 nspts (10, 180) = 16 nspts (12, 180) = 15 c nspts (2, 181) = 15 nspts (3, 181) = 16 nspts (5, 181) = 15 nspts (6, 181) = 16 nspts (7, 181) = 15 nspts (9, 181) = 16 nspts (10, 181) = 15 nspts (12, 181) = 16 c nspts (4, 182) = 3 nspts (5, 182) = 3 nspts (6, 182) = 3 nspts (10, 182) = 3 nspts (11, 182) = 3 nspts (12, 182) = 3 c c Cubic (23) c do 2600, i = 195, 199 nspsym (i) = 12 nsprot (2, i) = 4 nsprot (3, i) = 2 nsprot (4, i) = 3 nsprot (5, i) = 17 nsprot (6, i) = 18 nsprot (7, i) = 19 nsprot (8, i) = 20 nsprot (9, i) = 21 nsprot (10, i) = 22 nsprot (11, i) = 23 nsprot (12, i) = 24 2600 continue nsplat (196) = 5 nsplat (197) = 6 nsplat (199) = 6 c nspts (2, 198) = 6 nspts (3, 198) = 7 nspts (4, 198) = 5 nspts (6, 198) = 5 nspts (7, 198) = 6 nspts (8, 198) = 7 nspts (10, 198) = 7 nspts (11, 198) = 5 nspts (12, 198) = 6 c nspts (2, 199) = 6 nspts (3, 199) = 7 nspts (4, 199) = 5 nspts (6, 199) = 5 nspts (7, 199) = 6 nspts (8, 199) = 7 nspts (10, 199) = 7 nspts (11, 199) = 5 nspts (12, 199) = 6 c c Cubic (432) c do 2700, i = 207, 214 nspsym (i) = 24 nsprot (2, i) = 4 nsprot (3, i) = 2 nsprot (4, i) = 3 nsprot (5, i) = 17 nsprot (6, i) = 18 nsprot (7, i) = 19 nsprot (8, i) = 20 nsprot (9, i) = 21 nsprot (10, i) = 22 nsprot (11, i) = 23 nsprot (12, i) = 24 nsprot (13, i) = 7 nsprot (14, i) = 8 nsprot (15, i) = 25 nsprot (16, i) = 26 nsprot (17, i) = 27 nsprot (18, i) = 28 nsprot (19, i) = 29 nsprot (20, i) = 30 nsprot (21, i) = 31 nsprot (22, i) = 32 nsprot (23, i) = 33 nsprot (24, i) = 34 2700 continue nsplat (209) = 5 nsplat (210) = 5 nsplat (211) = 6 nsplat (214) = 6 c nspts (13, 208) = 10 nspts (14, 208) = 10 nspts (15, 208) = 10 nspts (16, 208) = 10 nspts (17, 208) = 10 nspts (18, 208) = 10 nspts (19, 208) = 10 nspts (20, 208) = 10 nspts (21, 208) = 10 nspts (22, 208) = 10 nspts (23, 208) = 10 nspts (24, 208) = 10 c nspts (2, 210) = 7 nspts (3, 210) = 5 nspts (4, 210) = 6 nspts (6, 210) = 6 nspts (7, 210) = 7 nspts (8, 210) = 5 nspts (10, 210) = 5 nspts (11, 210) = 6 nspts (12, 210) = 7 nspts (13, 210) = 19 nspts (14, 210) = 20 nspts (15, 210) = 21 nspts (16, 210) = 22 nspts (17, 210) = 19 nspts (18, 210) = 22 nspts (19, 210) = 20 nspts (20, 210) = 21 nspts (21, 210) = 19 nspts (22, 210) = 21 nspts (23, 210) = 22 nspts (24, 210) = 20 c nspts (2, 212) = 6 nspts (3, 212) = 7 nspts (4, 212) = 5 nspts (6, 212) = 5 nspts (7, 212) = 6 nspts (8, 212) = 7 nspts (10, 212) = 7 nspts (11, 212) = 5 nspts (12, 212) = 6 nspts (13, 212) = 21 nspts (14, 212) = 20 nspts (15, 212) = 22 nspts (16, 212) = 19 nspts (17, 212) = 21 nspts (18, 212) = 19 nspts (19, 212) = 20 nspts (20, 212) = 22 nspts (21, 212) = 21 nspts (22, 212) = 22 nspts (23, 212) = 19 nspts (24, 212) = 20 c nspts (2, 213) = 6 nspts (3, 213) = 7 nspts (4, 213) = 5 nspts (6, 213) = 5 nspts (7, 213) = 6 nspts (8, 213) = 7 nspts (10, 213) = 7 nspts (11, 213) = 5 nspts (12, 213) = 6 nspts (13, 213) = 23 nspts (14, 213) = 24 nspts (15, 213) = 25 nspts (16, 213) = 26 nspts (17, 213) = 23 nspts (18, 213) = 26 nspts (19, 213) = 24 nspts (20, 213) = 25 nspts (21, 213) = 23 nspts (22, 213) = 25 nspts (23, 213) = 26 nspts (24, 213) = 24 c nspts (2, 214) = 6 nspts (3, 214) = 7 nspts (4, 214) = 5 nspts (6, 214) = 5 nspts (7, 214) = 6 nspts (8, 214) = 7 nspts (10, 214) = 7 nspts (11, 214) = 5 nspts (12, 214) = 6 nspts (13, 214) = 23 nspts (14, 214) = 24 nspts (15, 214) = 25 nspts (16, 214) = 26 nspts (17, 214) = 23 nspts (18, 214) = 26 nspts (19, 214) = 24 nspts (20, 214) = 25 nspts (21, 214) = 23 nspts (22, 214) = 25 nspts (23, 214) = 26 nspts (24, 214) = 24 c return end c c----------------------------------------------------------------------- c subroutine initss c include 'main.cmn' c c Commands that are supported by the program right now c data commnd/'COMM', 'CELL', 'LATT', 'SYMM', 'MAPF', . 'PRIN', 'MAPR', 'SGLC', 'SOLN', 'CROS', . 'TOLE', 'DBLC', 'MAPT', 'NCRS', 'SUMF', . 'STOP', 'PATC', 'PATX', 'TITL', 'PATP', . 'SING', 'DOUB', 'SOLU', 'TRIP', 'SGLS', . 'DBLS', 'SREG', 'MULT', 'SYMS', 'SGLF', . 'DBLF', 'ORIG', 'MAPS', . 2*' '/ data ncom/33/ c data isym/0/, icent/1/, nsoln/0/, ncross/0/, isglct/5000/, . idblct/2500/, minhk/1/, minx/1/, ipatct/1500/, ntripl/10/, . itplsg/-5000/, itpldb/0/ data ntol/2, 2, 2/, nsglsh/0/, ndblsh/0/, nsymsh/0/ data nsglft/0/, ndblft/0/ data maptyp/'HASSP'/ data ncell/120, 120, 120/ data lmap/1/, log/6/, lprt/3/, lin/5/, lscrch/2/ data mapfil/'patmap.fft'/, prtfil/'patsol.prt'/ data nptpik/0/, mptpik/900/ data title/' PROGRAM PATSOL'/ data nncr/0/, isymop/2/ data sizmol/25.0/, qlcxpd/.true./, ordor/'BYBCX'/, euler/'ZXZ'/ data polar/'XYK'/, tlcrot/'EZXZ'/, alcrot/0.0, 0.0, 0.0/ data nlcaxs/0/ data icutlo/-20000/, locsho/0/, iorign/4, 4, 4/, locsit/0/ data namspg/' '/, numspg/0/ data qform/.true./, mapsec/1, 2, 3/ c data vershn/'PATSOL - V1.2 - 8/93 '/ c call inispg c patx = .false. pat = .false. sgl = .true. dbl = .true. sol = .false. map = .false. qmulti = .false. qckall = .false. qspher = .false. qmuncr = .true. qdump = .false. qsum = .false. c title (111:132) = vershn c dtr = acos (-1.0) / 180.0 c do 1000, i = 1, 2000 iworks (1, i) = -9999 1000 continue c do 1100, i = 1, maxsgl isgxyz (4, i) = -9999 1100 continue do 1200, i = 1, 100 isglpk (4, i) = -9999 1200 continue c do 1300, i = 1, 3 fixxyz (i) = .true. do 1400, j = 1, 8 isglor (i, j) = -1 1400 continue 1300 continue enan = .true. c do 1500, i = 1, maxcrs icross (1, i) = -999 1500 continue do 1600, i = 1, maxsln isoln (4, i) = -999 1600 continue c do 1700, i = 1, 3 nxyz (i) = 0 nsxyz (i) = -1 cell (i, 1) = 0.0 cell (i + 3, 1) = 90.0 1700 continue c do 1800, i = 1, maxlax tlcaxs (i) = 'P' 1800 continue c return end c c----------------------------------------------------------------------- c subroutine inside (ind) c c This routine determines whether the input position is in the input c Patterson map region. IND=0 if it is not, IND=1 if it is. c c IWORKS -- c 2 : input grid position, in cell of 2*Icell c include 'main.cmn' c ind = 0 do 1000, i = 1, 3 j = (iworks(i, 2) - 1) * ncell (i) / ( 2 * icell) j = j + 1 if (j .le. nxyz (i)) ind = ind + 1 1000 continue c if (ind .lt. 3) then ind = 0 else ind = 1 endif c return end c c----------------------------------------------------------------------- c subroutine intpol (ii, jj) c c This routine interpolates between the supplied grid points in the c Patterson map. c c SUBROUTINE CALLED : PAKIND c c IWORKS -- c ii : input position when ii > 0 c 150 : input position c 151-158 : positions used for interpolation c c JWORKS -- c -ii : input position when ii < 0 c include 'main.cmn' c kk = 0 if (ii .gt. 0) then do 1000, i = 1, 3 iworks (i, 150) = iworks (i, ii) kk = kk + mod (iworks (i, 150), 2) 1000 continue else do 1100, i = 1, 3 iworks (i, 150) = jworks (i, -ii) kk = kk + mod (iworks (i, 150), 2) 1100 continue endif c if (kk .eq. 3) then call pakind (150, jj) jj = ipatmp (jj) return c else if (kk .eq. 2) then nabor = 2 do 1300, i = 1, 3 if (mod(iworks(i, 150), 2) .eq. 1) then iworks (i, 151) = iworks (i, 150) iworks (i, 152) = iworks (i, 150) else iworks (i, 151) = iworks (i, 150) - 1 iworks (i, 152) = iworks (i, 150) + 1 endif 1300 continue c else if (kk .eq. 1) then nabor = 0 do 1400, i = 1, 3 if (mod(iworks(i, 150), 2) .eq. 1) then iworks (i, 151) = iworks (i, 150) iworks (i, 152) = iworks (i, 150) iworks (i, 153) = iworks (i, 150) iworks (i, 154) = iworks (i, 150) else if (nabor .eq. 0) then iworks (i, 151) = iworks (i, 150) - 1 iworks (i, 152) = iworks (i, 150) - 1 iworks (i, 153) = iworks (i, 150) - 1 iworks (i, 154) = iworks (i, 150) - 1 nabor = i else iworks (i, 151) = iworks (i, 150) - 1 iworks (i, 152) = iworks (i, 150) - 1 iworks (i, 153) = iworks (i, 150) - 1 iworks (i, 154) = iworks (i, 150) - 1 iworks (nabor, 152) = iworks (nabor, 152) + 2 iworks (nabor, 153) = iworks (nabor, 153) + 2 iworks (i, 153) = iworks (i, 153) + 2 iworks (i, 154) = iworks (i, 154) + 2 endif endif 1400 continue nabor = 4 c else if (kk .eq. 0) then do 1600, j = 151, 158 do 1500, i = 1, 3 iworks (i, j) = iworks (i, 150) - 1 1500 continue 1600 continue iworks (1, 152) = iworks (1, 152) + 2 iworks (1, 153) = iworks (1, 153) + 2 iworks (2, 153) = iworks (2, 153) + 2 iworks (2, 154) = iworks (2, 154) + 2 iworks (3, 155) = iworks (3, 155) + 2 iworks (1, 156) = iworks (1, 156) + 2 iworks (3, 156) = iworks (3, 156) + 2 iworks (1, 157) = iworks (1, 157) + 2 iworks (2, 157) = iworks (2, 157) + 2 iworks (3, 157) = iworks (2, 157) + 2 iworks (2, 158) = iworks (2, 158) + 2 iworks (3, 158) = iworks (3, 158) + 2 nabor = 8 endif c kk = 0 jj = 0 do 1700, i = 1, nabor do 1800, j = 1, 3 if (iworks(j, 150+i) .gt. nxyz(j)) go to 1700 1800 continue call pakind (150 + i, j) if (ipatmp(j) .lt. 0) go to 1700 kk = kk + 1 jj = jj + ipatmp (j) 1700 continue c if (kk .le. 0) then jj = -1 else jj = jj / kk endif c return end c c----------------------------------------------------------------------- c subroutine intprt (record, rekord, nchar) c include 'main.cmn' c character*(*) record, rekord character word*4 character latice(7)*1, xyz(3)*1 data latice/'P', 'A', 'B', 'C', 'F', 'I', 'R'/ data xyz/'A', 'B', 'C'/ c if (nw .eq. 0) return c j1 = iw (1, 1) j2 = iw (2, 1) j2 = min (j1+3, j2) word = record (j1:j2) c do 1000, i = 1, ncom if (word .eq. commnd (i)) go to 100 1000 continue write (log, 20) word return c 100 go to (101, 102, 103, 104, 105, 106, 107, 108, 109, 110, . 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, . 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, . 131, 132, 133, 134, 135) . i c c COMMent c 101 return c c CELL c 102 if (nw .lt. 4) then write (log, 40) commnd (i) call endgam endif do 1100, j = 1, 3 j1 = iw (1, j + 1) j2 = iw (2, j + 1) call iiread (record(j1:j2), j2 - j1 + 1, ixyz) ncell(j) = ixyz 1100 continue return c c LATTice c 103 if (nw .lt. 2) return j1 = iw (1, 2) j2 = iw (2, 2) do 1200, j = 1, 7 if (record(j1:j1) .eq. latice(j)) then icent = j return endif 1200 continue write (log, 30) record(j1:j2) return c c SYMMetry c 104 if (nw .lt. 2) return j1 = iw (1, 2) j2 = iw (2, nw) call symtns (record (j1:j2), j2-j1+1) return c c MAPFile c 105 if (nw .lt. 2) return j1 = iw (1, 2) j2 = iw (2, 2) mapfil = rekord (j1:j2) map = .true. return c c PRINtfile c 106 if (nw .lt. 2) return j1 = iw (1, 2) j2 = iw (2, 2) prtfil = rekord (j1:j2) return c c MAPRegion c 107 if (nw .lt. 4) then write (log, 40) commnd(i) call endgam endif do 1300, j = 1, 3 j1 = iw(1, j + 1) j2 = iw(2, j + 1) call rrread (record(j1:j2), j2-j1+1, txyz) if (txyz .ge. 2.0) then nxyz (j) = txyz + 0.5 else ixyz = txyz * 1000 + sign (0.5, txyz) nxyz (j) = 1000 * ixyz endif 1300 continue return c c SGLCut c 108 if (nw .lt. 2) return j1 = iw (1, 2) j2 = iw (2, 2) call rrread (record(j1:j2), j2-j1+1, txyz) isglct = 10000.0 * txyz if (nw .lt. 3) return j1 = iw (1, 3) j2 = iw (2, 3) call iiread (record(j1:j2), j2-j1+1, ixyz) if (ixyz .gt. 1) minhk = ixyz return c c SOLutioN (SOLN) c 109 if (nw .lt. 4) then write (log, 50) commnd(i) return endif if (nsoln .ge. maxsln) then write (log, 80) maxsln return endif nsoln = nsoln + 1 do 1400, j = 1, 3 j1 = iw(1, j + 1) j2 = iw(2, j + 1) call rrread (record(j1:j2), j2-j1+1, txyz) if (abs(txyz) .ge. 2.0) then isoln (j, nsoln) = txyz + sign(0.5, txyz) else ixyz = txyz * 1000 + sign (0.5, txyz) isoln (j, nsoln) = ixyz * 1000 endif 1400 continue return c c CROSs vectors c 110 if (nw .lt. 4) then write (log, 50) commnd(i) return endif if (ncross .ge. maxcrs) then write (log, 90) maxcrs return endif ncross = ncross + 1 do 1500, j = 1, 3 j1 = iw(1, j + 1) j2 = iw(2, j + 1) call rrread (record(j1:j2), j2-j1+1, txyz) if (abs(txyz) .ge. 2.0) then icross (j, ncross) = txyz + sign(0.5, txyz) else ixyz = txyz * 1000 + sign (0.5, txyz) icross (j, ncross) = ixyz * 1000 endif 1500 continue return c c TOLErance c 111 if (nw .lt. 4) then write (log, 50) commnd(i) return endif do 1600, j = 1, 3 j1 = iw(1, j + 1) j2 = iw(2, j + 1) call iiread (record(j1:j2), j2-j1+1, ixyz) if (ixyz .gt. 0) ntol (j) = ixyz 1600 continue return c c DBLCut c 112 if (nw .lt. 2) return j1 = iw (1, 2) j2 = iw (2, 2) call rrread (record(j1:j2), j2-j1+1, txyz) if (txyz .gt. 0.0) idblct = 10000.0 * txyz if (nw .lt. 3) return j1 = iw (1, 3) j2 = iw (2, 3) call iiread (record(j1:j2), j2-j1+1, ixyz) if (ixyz .gt. 1) minx = ixyz return c c MAPTtype c 113 if (nw .lt. 2) return j1 = iw (1, 2) j2 = iw (2, 2) j2 = min (j1 + 4, j2) maptyp = record (j1:j2) if (nw.lt.3 .or. maptyp.ne.'XPLOR') return j1 = iw (1, 3) if (record(j1:j1) .eq. 'F') qform = .false. return c c NCRSymmetry c 114 if (nw .lt. 2) return call ncrint (record, rekord, nchar) return c c SUMFunction c 115 if (nw .lt. 2) return j1 = iw (1, 2) if (record(j1:j1) .eq. 'T') then qsum = .true. else if (record(j1:j1) .eq. 'F') then qsum = .false. endif return c c STOP c 116 nchar = -1 return c c PATCut c 117 if (nw .lt. 2) return j1 = iw (1, 2) j2 = iw (2, 2) call rrread (record(j1:j2), j2-j1+1, txyz) if (txyz .gt. 0.4) ipatct = 10000.0 * txyz return c c PATX c 118 if (nw .lt. 2) return j1 = iw (1, 2) if (record(j1:j1) .eq. 'T') then patx = .true. else if (record(j1:j1) .eq. 'F') then patx = .false. endif return c c TITLe c 119 if (nw .eq. 1) return j1 = iw(2, nw) - iw(1, 2) j2 = 64 - j1 / 2 title (1:j2) = ' ' title (j2+1:j2+1+j1) = rekord (iw(1, 2):iw(2, nw)) title (j2+j1+2:132) = ' ' title (111:132) = vershn return c c PATPeak-search c 120 if (nw .lt. 2) return j1 = iw (1, 2) if (record(j1:j1) .eq. 'T') then pat = .true. else if (record(j1:j1) .eq. 'F') then pat = .false. endif return c c SINGle c 121 if (nw .lt. 2) return j1 = iw (1, 2) if (record(j1:j1) .eq. 'T') then sgl = .true. else if (record(j1:j1) .eq. 'F') then sgl = .false. endif return c c DOUBle c 122 if (nw .lt. 2) return j1 = iw (1, 2) if (record(j1:j1) .eq. 'T') then dbl = .true. else if (record(j1:j1) .eq. 'F') then dbl = .false. endif return c c SOLUtion c 123 if (nw .lt. 2) return j1 = iw (1, 2) if (record(j1:j1) .eq. 'T') then sol = .true. else if (record(j1:j1) .eq. 'F') then sol = .false. endif return c c TRIPle c 124 if (nw .lt. 2) return j1 = iw (1, 2) j2 = iw (2, 2) call iiread (record(j1:j2), j2-j1+1, ixyz) if (ixyz .gt. 0) ntripl = ixyz if (nw .lt. 3) return j1 = iw (1, 3) j2 = iw (2, 3) call rrread (record(j1:j2), j2-j1+1, txyz) if (txyz .ne. 0.0) itplsg = 10000.0 * txyz if (nw .lt. 4) return j1 = iw (1, 4) j2 = iw (2, 4) call rrread (record(j1:j2), j2-j1+1, txyz) if (txyz .ne. 0.0) itpldb = 10000.0 * txyz return c c SGLShow c 125 if (nw .lt. 4) then write (log, 50) commnd(i) return endif if (nsglsh .ge. maxsho) then write (log, 210) maxsho return endif nsglsh = nsglsh + 1 do 1700, j = 1, 3 j1 = iw(1, j + 1) j2 = iw(2, j + 1) call rrread (record(j1:j2), j2-j1+1, txyz) if (abs(txyz) .ge. 2.0) then isglsh (j, nsglsh) = txyz + sign(0.5, txyz) else ixyz = txyz * 1000 + sign (0.5, txyz) isglsh (j, nsglsh) = ixyz * 1000 endif 1700 continue return c c DBLShow c 126 if (nw .lt. 7) then write (log, 50) commnd(i) return endif if (ndblsh .ge. maxsho) then write (log, 220) maxsho return endif ndblsh = ndblsh + 1 do 1800, j = 1, 6 j1 = iw(1, j + 1) j2 = iw(2, j + 1) call rrread (record(j1:j2), j2-j1+1, txyz) if (abs(txyz) .ge. 2.0) then idblsh (j, ndblsh) = txyz + sign(0.5, txyz) else ixyz = txyz * 1000 + sign (0.5, txyz) idblsh (j, ndblsh) = ixyz * 1000 endif 1800 continue return c c c SREGion c 127 if (nw .lt. 4) then write (log, 50) commnd(i) return endif do 1900, j = 1, 3 j1 = iw(1, j + 1) j2 = iw(2, j + 1) call rrread (record(j1:j2), j2-j1+1, txyz) if (txyz .ge. 2.0) then nsxyz (j) = txyz + 0.5 else ixyz = txyz * 1000 + sign (0.5, txyz) nsxyz (j) = ixyz * 1000 endif 1900 continue return c c MULTiplicity correction c 128 if (nw .lt. 2) return j1 = iw (1, 2) if (record(j1:j1) .eq. 'T') then qmulti = .true. else if (record(j1:j1) .eq. 'F') then qmulti = .false. endif return c c SYMShow c 129 if (nw .lt. 4) then write (log, 50) commnd(i) return endif if (nsymsh .ge. maxsho) then write (log, 10) maxsho return endif nsymsh = nsymsh + 1 do 2000, j = 1, 3 j1 = iw(1, j + 1) j2 = iw(2, j + 1) call rrread (record(j1:j2), j2-j1+1, txyz) if (abs(txyz) .ge. 2.0) then isymsh (j, nsymsh) = txyz + sign(0.5, txyz) else ixyz = txyz * 1000 + sign (0.5, txyz) isymsh (j, nsymsh) = ixyz * 1000 endif 2000 continue return c c SGLFit c 130 if (nw .lt. 4) then write (log, 50) commnd(i) return endif if (nsglft .ge. maxsho) then write (log, 210) maxsho return endif nsglft = nsglft + 1 do 2100, j = 1, 3 j1 = iw(1, j + 1) j2 = iw(2, j + 1) call rrread (record(j1:j2), j2-j1+1, txyz) if (abs(txyz) .ge. 2.0) then isglft (j, nsglft) = txyz + sign(0.5, txyz) else ixyz = txyz * 1000 + sign (0.5, txyz) isglft (j, nsglft) = ixyz * 1000 endif 2100 continue return c c DBLFit c 131 if (nw .lt. 7) then write (log, 50) commnd(i) return endif if (ndblft .ge. maxsho) then write (log, 220) maxsho return endif ndblft = ndblft + 1 do 2200, j = 1, 6 j1 = iw(1, j + 1) j2 = iw(2, j + 1) call rrread (record(j1:j2), j2-j1+1, txyz) if (abs(txyz) .ge. 2.0) then idblft (j, ndblft) = txyz + sign(0.5, txyz) else ixyz = txyz * 1000 + sign (0.5, txyz) idblft (j, ndblft) = ixyz * 1000 endif 2200 continue return c c ORIGin c 132 if (nw .lt. 4) return do 2300, j = 1, 3 j1 = iw (1, j + 1) j2 = iw (2, j + 1) call iiread (record(j1:j2), j2-j1+1, ixyz) if (ixyz .gt. 0) iorign (i) = ixyz 2300 continue return c c MAPSection c 133 if (nw .lt. 2) return j1 = iw (1, 2) j2 = iw (2, 2) if (j2 .lt. j1+2) return do 2400, j = 1, 3 if (record(j1:j1) .eq. xyz(j)) then mapsec (3) = j go to 300 endif 2400 continue write (log, 60) record (j1:j1) call endgam 300 continue c j1 = j1 + 1 do 2500, j = 1, 3 if (record(j1:j1) .eq. xyz(j)) then mapsec (1) = j go to 400 endif 2500 continue write (log, 60) record (j1:j1) call endgam 400 continue c j1 = j1 + 1 do 2600, j = 1, 3 if (record(j1:j1) .eq. xyz(j)) then mapsec (2) = j go to 500 endif 2600 continue write (log, 60) record (j1:j1) call endgam 500 continue c if (mapsec(1).eq.mapsec(2) .or. mapsec(2).eq.mapsec(3) .or. . mapsec(1).eq.mapsec(3)) then write (log, 70) mapsec call endgam endif c 134 continue 135 continue c 10 format (' INTPRT> *Warning* ', 'Number of positions ', . 'reached the limit of ', . i4, /, ' INTPRT> ', 'This input position is ignored') 20 format (' INTPRT> *Warning* ', 'Unrecognized command > ', a, . '. IGNORED') 30 format (' INTPRT> *Warning* ', . 'Unknown Lattice Type > ', a, '

ASSUMED') 40 format (' INTPRT> *Fatal* ', 'Not enough arguments for ', . 'command > ', a, /, ' PROGRAM ABORT') 50 format (' INTPRT> *Warning* ', 'Not enough arguments for ', . 'command > ', a, '. IGNORED') 60 format (' INTPRT> *Fatal* ', 'Unrecognized sectioning index ', . a, /, ' PROGRAM ABORT') 70 format (' INTPRT> *Fatal* ', 'Illegal sectioning information ', . 3i4, /, ' PROGRAM ABORT') 80 format (' INTPRT> *Warning* ', 'Number of solutions reached ', . 'the limit of ', i4, /, ' INTPRT> ', . 'This input solution is ignored') 90 format (' INTPRT> *Warning* ', 'Number of cross-vectors ', . 'reached the limit of ', i4, /, ' INTPRT> ', . 'This input vector is ignored') 210 format (' INTPRT> *Warning* ', 'Number of single-sites ', . 'reached the limit of ', . i4, /, ' INTPRT> ', 'This input site is ignored') 220 format (' INTPRT> *Warning* ', 'Number of double-sites ', . 'reached the limit of ', . i4, /, ' INTPRT> ', 'These input sites are ignored') c return end c c----------------------------------------------------------------------- c subroutine locsym c include 'main.cmn' c dimension rlcang (3, maxncr) dimension a(3, 3), b(3, 3), c(3, 3) c c First, check out each of the input rotation axes - c if (nlcaxs .eq. 0) return c call nupage write (lprt, 10) write (lprt, 20) do 1000, i = 1, nlcaxs call ckaxis (i, th1, th2, th3) write (lprt, 30) i, th1, th2, th3, mlcaxs(i), . ((rlcrot(j, k, i), k=1, 3), j=1, 3) 1000 continue c c Now generate all other equivalent positions c c First, expand each rotation axis c nncr = nlcaxs if (.not. qlcxpd) go to 200 do 1100, i = 1, nlcaxs do 1200, j = 1, 3 do 1300, k = 1, 3 a (j, k) = rlcrot (j, k, i) b (j, k) = rlcrot (j, k, i) 1300 continue 1200 continue do 1400, j = 1, mlcaxs (i) - 1 if (j .gt. 1) then do 1500, m = 1, 3 do 1600, n = 1, 3 b (m, n) = c (m, n) 1600 continue 1500 continue endif call matmpy (a, b, c) c do 3200, k = 1, nncr do 3300, m = 1, 3 do 3400, n = 1, 3 b (m, n) = rlcrot (m, n, k) 3400 continue 3300 continue isame = 1 call smmatx (b, c, isame) if (isame .eq. 1) go to 1400 3200 continue c nncr = nncr + 1 if (nncr .gt. maxncr) then write (log, 40) maxncr call endgam endif do 1700, m = 1, 3 do 1800, n = 1, 3 rlcrot (m, n, nncr) = c (m, n) 1800 continue 1700 continue 1400 continue 1100 continue c c Second, generate the cross-rotation matrices between the c different symmetry operators. Repeat till we get no more. c icycle = 0 100 continue icycle = icycle + 1 mlcsym = nncr do 2000, i = 1, nncr do 2100, m = 1, 3 do 2200, n = 1, 3 a (m, n) = rlcrot (m, n, i) 2200 continue 2100 continue do 2300, j = 1, nncr if (i .eq. j) go to 2300 do 2400, m = 1, 3 do 2500, n = 1, 3 b (m, n) = rlcrot (m, n, j) 2500 continue 2400 continue c call matmpy (a, b, c) do 2600, k = 1, mlcsym do 2700, m = 1, 3 do 2800, n = 1, 3 b (m, n) = rlcrot (m, n, k) 2800 continue 2700 continue isame = 1 call smmatx (b, c, isame) if (isame .eq. 1) go to 2300 2600 continue c mlcsym = mlcsym + 1 if (mlcsym .gt. maxncr) then write (log, 40) maxncr c nncr = mlcsym c go to 200 call endgam endif do 2900, m = 1, 3 do 3000, n = 1, 3 rlcrot (m, n, mlcsym) = c (m, n) 3000 continue 2900 continue 2300 continue 2000 continue c if (mlcsym .ne. nncr) then nncr = mlcsym if (icycle .gt. 5) then write (log, 80) call endgam endif go to 100 endif 200 continue c c Find the identity symmetry operator, and put it as the first one c do 3500, i = 1, nncr do 3600, j = 1, 3 do 3700, k = 1, 3 if (j.ne.k .and. abs(rlcrot(j, k, i)).gt.0.001) . go to 3500 if (j.eq.k .and. abs(rlcrot(j, k, i)-1.0).gt.0.001) . go to 3500 3700 continue 3600 continue c do 3800, j = 1, 3 do 3900, k = 1, 3 rlcrot (j, k, i) = rlcrot (j, k, 1) rlcrot (j, k, 1) = 0.0 3900 continue rlcrot (j, j, 1) = 1.0 3800 continue go to 300 c 3500 continue c write (log, 90) nncr = nncr + 1 if (nncr .gt. maxncr) then write (log, 40) maxncr call endgam endif do 4000, j = 1, 3 do 4100, k = 1, 3 rlcrot (j, k, nncr) = rlcrot (j, k, 1) rlcrot (j, k, 1) = 0.0 4100 continue rlcrot (j, j, 1) = 1.0 4000 continue c 300 continue c c Now save all the rotation angles c do 4200, i = 1, nncr call ckaxis (i, th1, th2, th3) rlcang (1, i) = th1 rlcang (2, i) = th2 rlcang (3, i) = th3 4200 continue c c Now sort the angles, th3, th2, th1 order, ascending c dangle = 1.0 do 4300, i = 1, nncr do 4400, j = i + 1, nncr 400 continue th1 = rlcang (3, i) - rlcang (3, j) if (th1 .lt. 0.0) th1 = min (-dangle, th1) if (th1 .gt. 0.0) th1 = max (dangle, th1) if (th1 .lt. -dangle) go to 4400 th2 = rlcang (2, i) - rlcang (2, j) if (th2 .lt. 0.0) th2 = min (-dangle, th2) if (th2 .gt. 0.0) th2 = max (dangle, th2) if (th1.le.dangle .and. th2.lt.-dangle) go to 4400 th3 = rlcang (1, i) - rlcang (1, j) if (th3 .lt. 0.0) th3 = min (-dangle, th3) if (th3 .gt. 0.0) th3 = max (dangle, th3) if (th1.le.dangle .and. th2.le.dangle .and. . th3.lt.-dangle) go to 4400 do 4500, m = 1, 3 th1 = rlcang (m, i) rlcang (m, i) = rlcang (m, j) rlcang (m, j) = th1 do 4600, n = 1, 3 th1 = rlcrot (m, n, i) rlcrot (m, n, i) = rlcrot (m, n, j) rlcrot (m, n, j) = th1 4600 continue 4500 continue go to 400 4400 continue 4300 continue c c Now dump out all the symmetry operators c if (nncr+nlcaxs .gt. 50) call nupage write (lprt, 50) nncr, polar write (log, 110) nncr write (lprt, 60) do 3100, i = 1, nncr write (lprt, 70) i, (rlcang (k, i), k=1, 3), . ((rlcrot(j, k, i), k=1, 3), j=1, 3) 3100 continue c 10 format (5x, 'Listing of the Input Rotation Axes --', //) 20 format (7x, 'No.', 3x, 'Phi', 5x, 'Psi', 4x, 'Kappa', 1x, . 'Fold', 5x, 'R o t a t i o n M a t r i x', /) 30 format (7x, i3, 3f8.2, i4, 3(3x, 3f8.4)) 40 format (/, ' LOCSYM> *Fatal* ', 'Number of operators exceeds ', . 'the maximum [', i3, ']', /, ' PROGRAM ABORT') 50 format (//, 5x, 'Listing of the ', i3, ' Local Symmetry ', . 'Operators (Polar Angle Convention : ', a, ') --', //) 60 format (7x, 'No.', ' Phi ', ' Psi ', ' Kappa ', 10x, . 'R o t a t i o n M a t r i x', /) 70 format (7x, i3, 3f8.2, 4x, 3(3x, 3f8.4)) 80 format (/, ' LOCSYM> *Fatal* ', 'Too many iterations to ', . 'generate the local symmetry elements', /, . ' PROGRAM ABORT') 90 format (' LOCSYM> *Warning* ', 'No identity operation found in ', . 'the local symmetry', /, . ' LOCSYM> ', 'One will be inserted as the first ', . 'operation') 110 format (' LOCSYM> ', 'Number of local symmetry operators ', i3) c return end c c----------------------------------------------------------------------- c subroutine m3dinp c return end c c----------------------------------------------------------------------- c subroutine mapinp c c This is the driver routine for inputting Patterson map values. c Maps of different format can be treated by inserting a special c call here. c c The input map will be scaled such that the rms value is 100 c (the map average will be set to zero). Then the map will be c truncated between -1999 and 2000 (-20 to 20 rms). A constant c of 2000 is added to each grid point. c include 'main.cmn' c character*1 xyz(3) data xyz/'a', 'b', 'c'/ c i = nnxyz (1) * nnxyz (2) * nnxyz (3) if (i .gt. mapdim) then write (log, 20) i, mapdim call endgam endif c call gthead (1) c write (log, 80) c if (maptyp .eq. 'FSFOR') then call fsfinp (avg, sd, rmin, rmax) else if (maptyp .eq. 'HASSP') then call hspinp (avg, sd, rmin, rmax) else if (maptyp(1:4) .eq. 'CCP4') then call ccpinp (avg, sd, rmin, rmax) else if (maptyp(1:4) .eq. '3D01') then call m3dinp (avg, sd, rmin, rmax) else if (maptyp .eq. 'XPLOR') then call xplinp (avg, sd, rmin, rmax) else write (log, 10) maptyp call endgam endif c call nupage write (lprt, 30) mapfil write (lprt, 40) maptyp write (lprt, 90) xyz(mapsec(3)), xyz(mapsec(1)), xyz(mapsec(2)) c do 1000, i = 1, 3 write (lprt, 50) nxyz(i), ncell(i), xyz(i) 1000 continue write (lprt, 60) avg, sd, rmax, rmin c c Now scales the map make rms = 100.0 c if (sd .le. 0.0) then write (log, 70) call endgam endif c sd = 100.0 / sd do 1100, i = 1, nxyz(1) * nxyz(2) * nxyz(3) rmin = float (ipatmp (i)) / 1000.0 rmin = rmin * sd if (rmin .gt. 2000.0) rmin = 2000.0 if (rmin .le. -2000.0) rmin = -1999.0 ipatmp (i) = rmin + sign (0.5, rmin) ipatmp (i) = ipatmp (i) + 2000 1100 continue c 10 format (' MAPINP> *Fatal* ', 'Unknown map type -> ', a, /, . ' PROGRAM ABORT') 20 format (' MAPINP> *Fatal* ', 'Number of grid points [', i8, . '] exceeds the maximum [', i8, '].', /, ' PROGRAM ABORT') 30 format (5x, 'Input Patterson map file name -- ', a) 40 format (5x, 'Input Patterson map file format -- ', a, /) 50 format (5x, 'The map covers ', i3, ' grid points out of ', i3, . ' total along ', a) 60 format (/, 5x, 'The Patterson map has an average value of ', . f6.1, ' with a standard deviation of ', f6.1, /, . 5x, 'The maximum and minimum values of the map ', 2f9.1, . //, 5x, 'The average value will be subtracted from all ', . 'grid points.', /, 5x, 'Then the map values will be ', . 'scaled so that the standard deviation is 100.') 70 format (' MAPINP> *Fatal* ', 'Map standard deviation is 0', /, . ' PROGRAM ABORT') 80 format (' MAPINP> ', 'Reading in the Patterson map values') 90 format (5x, 'Input Patterson map has ', a, ' sections, ', a, . ' across, and ', a, ' down') c return end c c----------------------------------------------------------------------- c subroutine matinv (a, b) c include 'main.cmn' c dimension a(3, 3), b(3, 3), c(3, 3) c c b = a ** -1 c c To avoid over or underflow, divide all matrix elements by the c maximum element. c rmax = -1.0 do 2000, i = 1, 3 do 2100, j = 1, 3 rmax = max (rmax, abs(a (i, j))) 2100 continue 2000 continue c if (rmax .le. 0.0) then write (log, 20) call endgam endif c do 2200, i = 1, 3 do 2300, j = 1, 3 c (i, j) = a (i, j) / rmax 2300 continue 2200 continue c det = c(1, 1) * c(2, 2) * c(3, 3) + c(2, 1) * c(3, 2) * c(1, 3) . + c(1, 2) * c(2, 3) * c(3, 1) - c(1, 3) * c(2, 2) * c(3, 1) . - c(1, 2) * c(2, 1) * c(3, 3) - c(3, 2) * c(2, 3) * c(1, 1) if (abs(det) .lt. 1.0e-3) then write (log, 10) call endgam endif c b (1, 1) = c(2, 2) * c(3, 3) - c(3, 2) * c(2, 3) b (1, 2) = c(2, 1) * c(3, 3) - c(3, 1) * c(2, 3) b (1, 3) = c(2, 1) * c(3, 2) - c(3, 1) * c(2, 2) b (2, 1) = c(1, 2) * c(3, 3) - c(3, 2) * c(1, 3) b (2, 2) = c(1, 1) * c(3, 3) - c(3, 1) * c(1, 3) b (2, 3) = c(1, 1) * c(3, 2) - c(3, 1) * c(1, 2) b (3, 1) = c(1, 2) * c(2, 3) - c(2, 2) * c(1, 3) b (3, 2) = c(1, 1) * c(2, 3) - c(2, 1) * c(1, 3) b (3, 3) = c(1, 1) * c(2, 2) - c(2, 1) * c(1, 2) c b(1, 2) = - b(1, 2) b(2, 1) = - b(2, 1) b(3, 2) = - b(3, 2) b(2, 3) = - b(2, 3) c do 1000, i = 1, 3 b (i, i) = b (i, i) / det do 1100, j = i + 1, 3 tt = b (i, j) / det b (i, j) = b (j, i) / det b (j, i) = tt 1100 continue 1000 continue c do 1200, i = 1, 3 do 1300, j = 1, 3 b (i, j) = b (i, j) / rmax 1300 continue 1200 continue c 10 format (/, ' MATINV> *Fatal* ', 'Determinant of the matrix is 0 ', . /, ' PROGRAM ABORT') 20 format (/, ' MATINV> *Fatal* ', 'Maximum element of the matrix ', . 'is 0 ', /, ' PROGRAM ABORT') c return end c c----------------------------------------------------------------------- c subroutine matmpy (a, b, c) c dimension a(3, 3), b(3, 3), c(3, 3) c c C = A x B c do 1000, i = 1, 3 do 1100, j = 1, 3 c (i, j) = 0.0 do 1200, k = 1, 3 c (i, j) = c (i, j) + a (i, k) * b (k, j) 1200 continue 1100 continue 1000 continue c return end c c----------------------------------------------------------------------- c subroutine ncrdbl c c This routine carries out a cross-vector search with local symmetry. c c SUBROUTINES CALLED : CKLMTS, NCRMIN, SAVSGL, NCRPIK, NUPAGE c c IWORKS -- c 1 : search position c include 'main.cmn' c dimension linval(40) c if (.not. map) return c call cklmts (1) c jj = 0 c do 1000, i = 1, nlimts (1) if (qdump) then call nupage write (lprt, 50) i write (lprt, 60) (k, k=1, nlimts (3)) endif do 1100, j = 1, nlimts (2) do 1200, k = 1, nlimts (3) c iworks (1, 1) = i iworks (2, 1) = j iworks (3, 1) = k c call ncrmin (ival) if (qdump) linval (k) = max (ival - 2000, -99) c if (ival .lt. idblct) go to 1200 call savsgl (ival, jj) c 1200 continue if (qdump) write (lprt, 70) j, (linval(k), k=1, nlimts(3)) 1100 continue 1000 continue c c Peak search c call ncrpik c if (nsglpk .le. 0) then write (log, 10) return endif c call nupage write (lprt, 40) idblct - 2000 write (lprt, 20) nsglpk, min (nsglpk, 20) do 1300, n = 1, min (nsglpk, 20) i = isglpk (1, n) j = isglpk (2, n) k = isglpk (3, n) x = slimts (1, 1) + float (i - 1) * slimts (3, 1) y = slimts (1, 2) + float (j - 1) * slimts (3, 2) z = slimts (1, 3) + float (k - 1) * slimts (3, 3) write (lprt, 30) n, i, j, k, x, y, z, isglpk (4, n) - 2000 1300 continue c 10 format (' NCRDBL> ', 'No positions found in cross-vector search.') 20 format (5x, 'Number of positions found in the cross-vector ', . 'search ', i3, //, 5x, 'A list of the top ', i2, . ' positions follows --', //, 7x, 'No.', 2x, ' I ', . ' J ', ' K ', 8x, 'x', 9x, 'y', 9x, 'z', 6x, 'Height', . //) 30 format (7x, i3, 2x, 3i4, 3x, 3f10.4, 3x, i4) 40 format (5x, 'Cross-vector peak height cut-off ', i4, /) 50 format (5x, 'Dumping the cross-vector search map -- ', //, . 5x, 'Section I = ', i3, ' with K across and J down') 60 format (//, 9x, 40i3) 70 format (4x, i3, 2x, 40i3) c return end c c----------------------------------------------------------------------- c subroutine ncrint (record, rekord, nchar) c include 'main.cmn' c character*(*) record, rekord character word*4 c if (nw .lt. 3) return c j1 = iw (1, 2) j2 = iw (2, 2) j2 = min (j1 + 3, j2) word = record (j1:j2) c if (word .eq. 'CHEC') then j1 = iw (1, 3) if (record(j1:j1) .eq. 'T') qckall = .true. if (record(j1:j1) .eq. 'F') qckall = .false. else if (word .eq. 'CELL') then if (nw .lt. 5) then write (log, 30) 'CELL' return endif do 1300, j = 1, min (6, nw - 2) j1 = iw (1, j + 2) j2 = iw (2, j + 2) call rrread (record(j1:j2), j2-j1+1, txyz) cell (j, 1) = txyz 1300 continue else if (word .eq. 'CUTL') then j1 = iw (1, 3) j2 = iw (2, 3) call rrread (record(j1:j2), j2-j1+1, txyz) icutlo = txyz * 10000.0 else if (word .eq. 'DUMP') then j1 = iw (1, 3) if (record(j1:j1) .eq. 'T') qdump = .true. if (record(j1:j1) .eq. 'F') qdump = .false. else if (word .eq. 'EULE') then j1 = iw (1, 3) j2 = iw (2, 3) euler = record (j1:j2) else if (word .eq. 'LOCE') then j1 = iw (1, 3) if (record(j1:j1) .eq. 'T') qlcxpd = .true. if (record(j1:j1) .eq. 'F') qlcxpd = .false. else if (word .eq. 'LOCR') then if (nw .lt. 5) then write (log, 30) 'LOCRotate' return endif do 1500, j = 1, 3 j1 = iw (1, j + 2) j2 = iw (2, j + 2) call rrread (record(j1:j2), j2-j1+1, txyz) alcrot (j) = txyz 1500 continue if (nw .ge. 6) then j1 = iw (1, 6) tlcrot (1:1) = record (j1:j1) endif call gtmatx (2) else if (word .eq. 'LOCS') then if (nw .lt. 5) then write (log, 30) 'LOCSymmetry' call endgam endif nlcaxs = nlcaxs + 1 if (nlcaxs .gt. maxlax) then write (log, 40) maxlax call endgam endif do 1600, j = 1, 3 j1 = iw (1, j + 2) j2 = iw (2, j + 2) call rrread (record(j1:j2), j2 - j1 + 1, txyz) rlcaxs (j, nlcaxs) = txyz 1600 continue if (nw .eq. 5) go to 200 j1 = iw (1, 6) j2 = iw (2, 6) call rrread (record(j1:j2), j2 - j1 + 1, txyz) if (txyz .gt. 0.0) then if (txyz .lt. 10.0) then mlcaxs (nlcaxs) = txyz else mlcaxs (nlcaxs) = 360.0 / mod(txyz, 360.0) + 0.5 endif endif if (nw .eq. 6) go to 200 j1 = iw (1, 7) tlcaxs (nlcaxs) = record (j1:j1) 200 continue call gtmatx (1) else if (word .eq. 'MULT') then j1 = iw (1, 3) if (record(j1:j1) .eq. 'T') qmuncr = .true. if (record(j1:j1) .eq. 'F') qmuncr = .false. else if (word .eq. 'ORTH') then j1 = iw (1, 3) j2 = iw (2, 3) ordor = record (j1:j2) else if (word .eq. 'POLA') then j1 = iw (1, 3) j2 = iw (2, 3) polar = record (j1:j2) else if (word .eq. 'PXYZ') then if (nw .lt. 5) return do 1100, i = 1, 3 j1 = iw (1, i + 2) j2 = iw (2, i + 2) call rrread (record(j1:j2), j2-j1+1, txyz) pinput (i) = txyz 1100 continue else if (word .eq. 'SHOW') then if (nw .lt. 5) return if (locsho .ge. maxsho) then write (log, 10) maxsho return endif locsho = locsho + 1 do 1200, i = 1, 3 j1 = iw (1, i + 2) j2 = iw (2, i + 2) call rrread (record(j1:j2), j2-j1+1, txyz) shoncr (i, locsho) = txyz 1200 continue else if (word .eq. 'SITE') then if (nw .lt. 5) return if (locsit .ge. maxsho) then write (log, 10) maxsho return endif locsit = locsit + 1 do 1400, i = 1, 3 j1 = iw (1, i + 2) j2 = iw (2, i + 2) call rrread (record(j1:j2), j2-j1+1, txyz) sitncr (i, locsit) = txyz 1400 continue else if (word .eq. 'SIZE') then j1 = iw (1, 3) j2 = iw (2, 3) call rrread (record(j1:j2), j2-j1+1, txyz) if (txyz .gt. 0.0) sizmol = txyz else if (word .eq. 'SLIM') then if (nw .lt. 6) return j1 = iw (1, 3) j2 = iw (2, 3) call iiread (record(j1:j2), j2-j1+1, ixyz) if (ixyz .lt. 1) return if (ixyz .gt. 3) return do 1000, i = 1, 3 j1 = iw (1, i + 3) j2 = iw (2, i + 3) call rrread (record(j1:j2), j2-j1+1, txyz) slimts (i, ixyz) = txyz 1000 continue else if (word .eq. 'SPHE') then j1 = iw (1, 3) if (record(j1:j1) .eq. 'T') qspher = .true. if (record(j1:j1) .eq. 'F') qspher = .false. else if (word .eq. 'SYMM') then j1 = iw (1, 3) j2 = iw (2, 3) call iiread (record(j1:j2), j2-j1+1, ixyz) if (ixyz .gt. 1) isymop = ixyz else write (log, 20) word endif c 10 format (' NCRINT> *Warning* ', 'Number of sites for SHOW ', . 'exceeds the maximum of ', i2, ' -- IGNORED') 20 format (' NCRINT> *Warning* ', 'Unrecognized option for ', . 'NCRSymmetry --> ', a) 30 format (' NCRINT> *Warning* ', 'Not enough arguments for ', . 'option -> ', a, '. IGNORED') 40 format (' NCRINT> *Fatal* ', 'Number of local symmetry axes ', . 'exceeds the maximum of ', i2, /, ' PROGRAM ABORT') c return end