The length of the prefix string to use.
C
C
C
C
C close_ Close the creation CIF. MUST be used if pfile_ is used.
C [subroutine call]
C
C
C____________________________________________________________________________
C
C
C
C....The CIF tool box also provides variables for data access control:
C
C
C alias_ Logical variable: if left .true. then all calls to
C CIFtbx functions may use aliases of data item names.
C The preferred synonym from the dictionary will be
C subsituted internally, provided aliased data names were
C supplied by an input dictionary (via dict_). The
C default is .true., but alias_ may be set to .false.
C in an application.
C
C aliaso_ Logical variable: if set .true. then cif output
C routines will convert aliases to the names to preferred
C synonyms from the dictionary. The default is .false., but
C aliaso_ may be set to .true. in an application. The
C setting of aliaso_ is independent of the setting of
C alias_.
C
C align_ Logical variable signals alignment of loop_ lists during
C the creation of a CIF. The default is .true.
C
C bloc_ Character*(NUMCHAR) variable: the current block name.
C
C dictype_ Character*(NUMCHAR) variable: the precise data type code
C (see test_)
C
C diccat_ Character*(NUMCHAR) variable: the category (see test_)
C
C dicname_ Character*(NUMCHAR) variable: the root alias (see test_)
C
C esdlim_ Integer variable: Specifies the upper limit of esd's
C produced by pnumb_, and, implicitly, the lower limit.
C The default value is 19, which limits esd's to the range
C 2-19. Typical values of esdlim_ might be 9 (limiting
C esd's to the range 1-9), 19, or 29 (limiting esd's
C to the range 3-29)
C
C file_ Character*(MAXBUF) variable: the filename of the current file.
C
C line_ Integer variable: Specifies the input/output line limit
C for processing a CIF. The default value is 80 characters.
C This may be set by the program. The max value is MAXBUF
C which has a default value of 200.
C
C list_ Integer variable: the loop block number (see test_).
C
C long_ Integer variable: the length of the data string in strg_.
C
C longf_ Integer variable: the length of the filename in file_.
C
C loop_ Logical variable signals if another loop packet is present.
C
C pposdec_ Integer variable giving the position of the decimal point
C for the next number to be written.
C
C pposend_ Integer variable giving the ending column of the next
C number or quoted character value to be written. Used to
C pad with zeros or blanks.
C
C pposnam_ Integer variable giving the starting column of the next
C name or comment or data block to be written.
C
C pposval_ Integer variable giving the starting column of the next
C data value to be written by pchar_, pnumb_ or pnumd_.
C Also used to set the position of the initial "loop_"
C in a ploop_ call or to set the position of a terminal "save_"
C for a save frame in a pdata_ call for which saveo_ is .true.
C
C posdec_ Integer variable giving the position of the decimal point
C for the last number read.
C
C posend_ Integer variable giving the ending column of the last
C data value read, not including a terminal quote.
C
C posnam_ Integer variable giving the starting column of the last
C name or comment or data block read.
C
C posval_ Integer variable giving the starting column of the last
C data value read. Also reports the column of the
C terminal "save_" of a save frame.
C
C pquote_ Character variable giving the quotation symbol to be
C used for the next string written.
C
C precn_ Integer variable: Reports the record number of the last
C line written to the output cif. Set to zero by init_. Also
C set to zero by pfile_ and close_ if the output cif file name
C was not blank.
C
C ptabx_ Logical variable signals tab character expansion to blanks
C during the creation of a CIF. The default is .true.
C
C quote_ Character variable giving the quotation symbol found
C delimiting the last string read.
C
C recn_ Integer variable: Reports the record number of the last
C line read from the direct access copy of the input cif.
C
C save_ Logical variable signals that the current data block
C is actually a save-frame (.true. for a save-frame).
C
C saveo_ Logical variable signals that the output data block from
C pdata_ is actually a save-frame (.true. for a save-frame).
C
C strg_ Character*(MAXBUF) variable: the current data item.
C
C tabl_ Logical variable signals tab-stop alignment of output
C during the creation of a CIF. The default is .true.
C
C tabx_ Logical variable signals tab character expansion to blanks
C during the reading of a CIF. The default is .true.
C
C text_ Logical variable signals if another text line is present.
C
C type_ Character*4 variable: the data type code (see test_).
C
C
C
C_____________________________________________________________________________
C
C
C >>>>>> Set the device numbers.
C
function init_(devcif,devout,devdir,deverr)
C
logical init_
include 'ciftbx.sys'
integer devcif,devout,devdir,deverr
integer ii,kdig
real ytest
double precision ztest
C
init_=.true.
cifdev=devcif
outdev=devout
dirdev=devdir
errdev=deverr
recn_=0
precn_=0
C
C recompute decimal single precision precision
C This is found by computing the smallest power of
C 10 which, when added to 1, produces a change
C and then backing off by 1
C
decprc = .1
do ii = 1,6
ytest = 1.+decprc/10.
if (ytest.eq.1.) go to 100
decprc = decprc/10.
enddo
100 continue
decprc=decprc*10.
C
C recompute decimal double precision precision
C
kdig = 1
dpprc = .1D0
do ii = 1,15
ztest = 1.D0+dpprc/10.
if (ztest.eq.1.D0) go to 200
dpprc = dpprc/10.D0
kdig = kdig+1
enddo
200 continue
dpprc=dpprc*10.D0
write(ndpfmt,'(5h(d30.,i2,1h))') kdig-1
C
C recompute decimal single precision minimum power of ten
C
decmin = .1
do ii = 1,37
ytest = decmin/10.
if (ytest.eq.0.) go to 300
decmin = decmin/10.
enddo
300 continue
C
C recompute decimal double precision minimum power of 10
C and its log base 10 (minexp)
C
dpmin = .1D0
minexp = -1
do ii = 1,307
ztest = dpmin/10.
if (ztest.eq.0.D0) go to 400
dpmin = dpmin/10.D0
minexp = minexp-1
enddo
400 continue
call clearfp
return
end
C
C
C
C
C
C >>>>>> Read a CIF dictionary and prepare for checks
C
function dict_(fname,checks)
C
logical dict_
logical ocif_
logical data_
logical char_
integer lastnb
include 'ciftbx.sys'
character locase*(MAXBUF)
character fname*(*),checks*(*)
character temp*24,codes(4)*5,name*(MAXBUF),bxname*(NUMCHAR)
character bcname*(NUMCHAR),biname*(NUMCHAR),bname*(NUMCHAR)
character baname*(NUMCHAR),ganame*(NUMCHAR),btname*(NUMCHAR)
character batag*(NUMCHAR)
integer lbcname,lbaname,lbtname,lbname
integer kdict,kadict,ifind,jfind,iafind
integer i,j,nmatch,mycat,ksmatch,ii
C
C Control flags for matching categories, names and types
C
C icloop is the loop number of the block for the
C current category
C ictype is the type of the current category
C 0 - none found yet
C 1 - _item.category.id
C 2 - _category
C 3 - _category.id
C inloop is the loop number of the block for the
C current name
C intype is the type of the current name
C 0 - none found yet
C 1 - _item.name
C 2 - _name
C ialoop is the loop number of the block for the
C current alias
C iatype is the type for the current alias
C 0 - none found yet
C 1 - _item_aliases.alias_name
C itloop is the loop number of the block for the
C current type
C ittype is the type of the current type
C 0 - none found yet
C 1 - _item_type.code
C 2 - _type
C
integer icloop,ictype,inloop,intype,ialoop,iatype,
* itloop,ittype
C
character*4 map_type(12),map_to(12),mapped
character*(NUMCHAR) dt(2),ct(3),nt(2),at(1),tt(2)
data map_type
* /'floa','int ','yyyy','symo','ucha','ucod','name','idna',
* 'any ','code','line','ulin'/
data map_to
* /'numb','numb','char','char','char','char','char','char',
* 'char','char','char','char'/
data dt
* /'_dictionary.title ',
* '_dictionary_name '/
data ct
* /'_item.category_id ',
* '_category ',
* '_category.id '/
data nt
* /'_item.name ',
* '_name '/
data at
* /'_item_aliases.alias_name '/
data tt
* /'_item_type.code ',
* '_type '/
C
data codes /'valid','dtype','reset','close'/
C
C....... Are the codes OK
C
temp=checks
i=0
120 i=i+1
if(i.ge.24) goto 190
if(temp(i:i).eq.' ') goto 120
do 150 j=1,4
if(temp(i:i+4).eq.codes(j)) goto 170
150 continue
dict_=.false.
goto 500
170 i=i+4
if(j.eq.1) vcheck='yes'
if(j.eq.2) tcheck='yes'
if(j.eq.3) then
vcheck = 'no '
tcheck = 'no '
goto 170
endif
if(j.eq.4) then
vcheck = 'no '
tcheck = 'no '
ndcname = 0
ndict = 0
if(nname.gt.0) then
do 180 i = 1,nname
dtype(i)=' '
dxtyp(i)=' '
cindex(i)=0
ddict(i)=0
180 continue
endif
dict_=.true.
goto 500
endif
goto 120
C
C if no category names have been loaded, clean up
C the hash table for dictionary category names
C
190 if(ndcname.eq.0) then
call hash_init(dcname,dcchain,NUMDICT,ndcname,dchash,
* NUMHASH)
endif
C
C if no dictionary names have been loaded, clean up
C the hash table for dictionary names
C
if(ndict.eq.0) then
call hash_init(dicnam,dicchain,NUMDICT,ndict,dichash,
* NUMHASH)
endif
C
C....... Open and store the dictionary
C
dict_=.true.
if(fname.eq.' ') goto 500
if(nname.gt.0) call err(' Dict_ must precede ocif_')
dict_=ocif_(fname)
if(.not.dict_) goto 500
dictfl='yes'
C
C....... Loop over data blocks; extract _name's, _type etc.
C
200 if(.not.data_(' ')) goto 400
if(bloc_(1:1).eq.'_') then
bname=locase(bloc_)
else
bname='_'//locase(bloc_)
endif
lbname=lastnb(bname)
C
C see if this is a dictionary defining block
C
do i = 1,2
if(char_(dt(i),name)) goto 200
enddo
C
Cdbg WRITE(6,*) ndict,bloc_
C
C Analyze loop structure for categories, names and types
C
C
C initalize loop info
C
icloop = -1
inloop = -1
ialoop = -1
itloop = -1
ictype = 0
intype = 0
iatype = 0
ittype = 0
bcname = ' '
lbcname = 1
baname = ' '
batag = ' '
lbaname = 1
btname = ' '
lbtname = 1
biname=bloc_
mycat=0
loop_=.false.
loopnl=0
nmatch=0
ksmatch=0
C
C Process categories
C
do i = 1,3
if(char_(ct(i),name)) then
if(ictype.ne.0)
* call warn(' Multiple DDL 1 and 2 category definitions ')
ictype = i
if(loop_) icloop = loopnl
bcname=locase(name(1:long_))
lbcname=long_
call hash_store(bcname,
* dcname,dcchain,NUMDICT,ndcname,dchash,NUMHASH,mycat)
if(mycat.eq.0) then
call err(' Dictionary category names > NUMDICT ')
endif
C
C if this is not a loop of categories, we expect a match
C against the block name
C
if(.not.loop_) then
if(ictype.eq.1) then
if(bname(1:lbcname+2).ne.
* '_'//bcname(1:lbcname)//'.') then
call warn(' Category id does not match block name')
endif
else
if(ictype.eq.2) then
if(bcname.ne.'dictionary_definition') then
if(bname(1:lbcname+2).ne.
* '_'//bcname(1:lbcname)//'_') then
if(bname(1:lbcname+2).ne.
* '_'//bcname(1:lbcname)//' ') then
call warn(' Category id does not match block name')
endif
endif
endif
endif
endif
endif
endif
loop_ = .false.
loopnl = 0
enddo
C
C Process names
do i = 1,2
if(char_(nt(i),name)) then
if(intype.ne.0)
* call warn(' Multiple DDL 1 and 2 name definitions ')
intype = i
bxname=locase(name(1:long_))
if(loop_) inloop = loopnl
endif
loop_ = .false.
loopnl=0
enddo
if(intype.eq.0.and.ictype.ne.3)
* call warn (' No name defined in block')
loop_ = .false.
if(char_(at(1),name)) then
iatype=1
baname = locase(name(1:long_))
batag = name(1:long_)
lbaname = long_
if(loop_) ialoop = loopnl
endif
loop_ = .false.
loopnl=0
if(ictype.ne.3) then
do i=1,2
if(char_(tt(i),name)) then
if(ittype.ne.0)
* call warn(' Multiple DDL 1 and 2 type definitions ')
ittype = i
btname = locase(name(1:long_))
if(loop_) itloop = loopnl
endif
loop_ = .false.
loopnl=0
enddo
endif
C
C Now test for consistent combinations
C
if(inloop.ne.-1) then
if(icloop.ne.-1.and.icloop.ne.inloop)
* call warn(
* ' Categories and names in different loops')
if(iatype.ne.0.and.ialoop.ne.inloop) then
if(ialoop.eq.-1) then
if(bxname.ne.bname)
* call warn(
* ' One alias, looped names, linking to first')
else
call warn(
* ' Aliases and names in different loops '
* //' only using first alias ')
endif
endif
if(itloop.ne.-1.and.itloop.ne.inloop)
* call warn(
* ' Types and names in different loops')
else
if(icloop.ne.-1)
* call warn(
* ' Multiple categories for one name')
if(itloop.ne.-1)
* call warn(
* ' Multiple types for one name')
endif
C
C This is the main loop
C
if(intype.eq.0) go to 200
250 if(.not.char_(nt(intype),name)) goto 200
kdict=ndict+1
call hash_store(locase(name(1:long_)),dicnam,dicchain,
* NUMDICT,ndict,dichash,NUMHASH,ifind)
if(ifind.eq.0) call err(' Cifdic names > NUMDICT')
if(ifind.eq.kdict)dictag(ifind)=name(1:long_)
if(dicnam(ifind).eq.bname) nmatch=ifind
if(dicnam(ifind)(1:lbname).eq.bname) ksmatch=ifind
Cdbg if(dicnam(ifind).ne.bname)
Cdbg * call warn (' Name mismatch: '//dicnam(ifind)//bname)
if(inloop.ge.0)then
C
C We are in a loop of names. If it is the same loop as
C for categories, we need to extract the matching category
C
if(inloop.eq.icloop) then
mycat=0
if(char_(ct(ictype),name)) then
bcname=locase(name(1:long_))
lbcname=long_
call hash_store(bcname,
* dcname,dcchain,NUMDICT,ndcname,dchash,NUMHASH,mycat)
if(mycat.eq.0) then
call err(' Dictionary category names > NUMDICT ')
endif
endif
endif
C
C If it is the same loop as for types, we need to extract
C the matching type
C
if(inloop.eq.itloop) then
btname=' '
if(char_(ct(ittype),name)) then
btname=locase(name(1:long_))
lbtname=long_
endif
endif
C
C If it is the same loop as for aliases, we need to extract
C the matching alias
C
if(inloop.eq.ialoop) then
baname=' '
batag=' '
if(char_(at(1),name)) then
baname = locase(name(1:long_))
batag = name(1:long_)
lbaname = long_
endif
endif
endif
C
C now we have a name stored in dicnam at location ifind
C the index of the category in mycat, the type in btname,
C the alias in baname
C
C First verify match between the name and category, if
C we have one, or extract from the block name
C
if (mycat.eq.0) then
if (dcindex(ifind).eq.0) then
if (dicnam(ifind).eq.bloc_) then
call excat(dicnam(ifind),bcname,lbcname)
Cdbg call warn(' Extracting category name from block name '
Cdbg * //bloc_(1:max(1,lastnb(bloc_))))
if(bcname(1:1).ne.' ') then
ictype = 1
call hash_store(bcname,
* dcname,dcchain,NUMDICT,ndcname,dchash,NUMHASH,mycat)
if(mycat.eq.0) then
call err(' Dictionary category names > NUMDICT ')
endif
else
call warn(' No category defined in block '
* //bloc_(1:max(1,lastnb(bloc_)))//' and name '
* //dicnam(ifind)(1:max(1,lastnb(dicnam(ifind))))
* //' does not match')
endif
endif
endif
else
if (bcname(1:lbcname).ne.'dictionary_definition') then
if (dicnam(ifind)(1:lbcname+1).ne.'_'//bcname(1:lbcname)
* .or.( dicnam(ifind)(lbcname+2:lbcname+2).ne.'_' .and.
* dicnam(ifind)(lbcname+2:lbcname+2).ne.'.' .and.
* dicnam(ifind)(lbcname+2:lbcname+2).ne.' ' )) then
call warn(' Item name '//
* dicnam(ifind)(1:max(1,lastnb(dicnam(ifind))))//' '//
* ' does not match category name '//bcname(1:lbcname))
endif
endif
endif
C
C We will need the type in what follows. cifdic.m96 defines
C some higher level types. We map them to primitive types
C
mapped = btname(1:4)
do i = 1,12
if (btname(1:4).eq.map_type(i)) mapped = map_to(i)
enddo
if (mapped.ne.'char' .and.
* mapped.ne.'text' .and.
* mapped.ne.' ' .and.
* mapped.ne.'null' .and.
* mapped.ne.'numb' ) then
if (tcheck .eq. 'yes') call warn (' Item type '//
* btname(1:max(1,lastnb(btname)))//' not recognized')
endif
C
C There are two cases to consider, one if the name is new to
C the dictionary, the other, if it is not
C
if(ifind.eq.kdict) then
aroot(ifind)=0
alias(ifind)=0
dcindex(ifind)=mycat
dictyp(ifind)=mapped
dicxtyp(ifind)=btname
else
if(dcindex(ifind).ne.mycat) then
if(dcindex(ifind).eq.0) then
jfind=ifind
if (aroot(ifind).ne.0) jfind=ifind
255 continue
dcindex(jfind)=mycat
jfind=alias(jfind)
if(jfind.ne.0) goto 255
else
if(mycat.ne.0.and.
* (vcheck.eq.'yes'.or.tcheck.eq.'yes'))
* call warn(' Attempt to redefine category for item')
endif
endif
if(dictyp(ifind).ne.mapped .or.
* dicxtyp(ifind).ne.btname) then
if(dictyp(ifind).eq.' ') then
jfind=ifind
if (aroot(ifind).ne.0) jfind=ifind
256 continue
dictyp(jfind)=mapped
dicxtyp(jfind)=btname
jfind=alias(jfind)
if(jfind.ne.0) go to 256
else
if(mapped.ne.' '.and.tcheck.eq.'yes')
* call warn(' Attempt to redefine type for item')
endif
endif
endif
C
C now deal with alias, if any.
C
if(baname.ne.' ') then
kadict=ndict+1
call hash_store(baname(1:lbaname),dicnam,dicchain,
* NUMDICT,ndict,dichash,NUMHASH,iafind)
if(iafind.eq.0) call err(' Cifdic names > NUMDICT')
if(iafind.eq.kadict) then
dictag(iafind) =batag
aroot(iafind) =aroot(ifind)
if(aroot(iafind).eq.0) aroot(iafind)=ifind
alias(iafind) =0
alias(ifind) =iafind
dcindex(iafind) =dcindex(ifind)
dictyp(iafind) =dictyp(ifind)
dicxtyp(iafind) =dicxtyp(ifind)
else
if(aroot(iafind).ne.0) then
if(aroot(iafind).eq.ifind .or.
* aroot(iafind).eq.aroot(ifind)) then
call warn(' Duplicate definition of same alias')
else
call warn(' Conflicting definition of alias')
endif
else
if((dcindex(iafind).eq.0.or.
* dcindex(iafind).eq.dcindex(ifind)).and.
* (dictyp(iafind).eq.' '.or.
* (dictyp(iafind).eq.dictyp(ifind) .and.
* dicxtyp(iafind).eq.dicxtyp(ifind)))) then
dcindex(iafind) =dcindex(ifind)
dictyp(iafind) =dictyp(ifind)
dicxtyp(iafind) =dicxtyp(ifind)
endif
aroot(iafind) =aroot(ifind)
if(aroot(iafind).eq.0) aroot(iafind)=ifind
alias(ifind) =iafind
endif
endif
endif
if(inloop.ge.0) then
baname = ' '
batag = ' '
endif
C
if(inloop.ge.0.and.loop_) go to 250
if(nmatch.eq.0) then
if (ksmatch.eq.0.or.inloop.lt.0) then
call warn(' No name in the block matches the block name')
endif
endif
C
C check for aliases
C we execute this loop only in the case of unlooped name
C with looped alias
C
if(inloop.lt.0.and.ialoop.ge.0) then
loop_=.false.
loopnl=0
ganame=baname
260 if(.not.char_(at(iatype),name)) goto 200
baname=locase(name(1:long_))
batag=name(1:long_)
lbaname=long_
if(baname.eq.ganame) then
if(loop_) go to 260
go to 200
endif
if(baname.ne.' ') then
kadict=ndict+1
call hash_store(baname(1:lbaname),dicnam,dicchain,
* NUMDICT,ndict,dichash,NUMHASH,iafind)
if(iafind.eq.0) call err(' Cifdic names > NUMDICT')
if(iafind.eq.kadict) then
dictag(iafind) =batag
aroot(iafind) =aroot(ifind)
if(aroot(iafind).eq.0) aroot(iafind)=ifind
alias(iafind) =0
alias(ifind) =iafind
dcindex(iafind) =dcindex(ifind)
dictyp(iafind) =dictyp(ifind)
dicxtyp(iafind) =dicxtyp(ifind)
ifind=iafind
else
if(aroot(iafind).ne.0) then
if(aroot(iafind).eq.ifind .or.
* aroot(iafind).eq.aroot(ifind)) then
call warn(' Duplicate definition of same alias')
else
call warn(' Conflicting definition of alias')
endif
else
if((dcindex(iafind).eq.0.or.
* dcindex(iafind).eq.dcindex(ifind)).and.
* (dictyp(iafind).eq.' '.or.
* (dictyp(iafind).eq.dictyp(ifind) .and.
* dicxtyp(iafind).eq.dicxtyp(ifind)))) then
dcindex(iafind) =dcindex(ifind)
dictyp(iafind) =dictyp(ifind)
dicxtyp(iafind) =dicxtyp(ifind)
ifind=iafind
endif
aroot(iafind) =aroot(ifind)
if(aroot(iafind).eq.0) aroot(iafind)=ifind
alias(ifind) =iafind
endif
endif
endif
if(loop_) go to 260
endif
go to 200
C
400 bloc_=' '
if (ndcname.ne.0) then
do ii = 1,ndict
if (aroot(ii).eq.0.and.dcindex(ii).eq.0)
* call warn(' No category specified for name '//
* dicnam(ii)(1:max(1,lastnb(dicnam(ii)))))
enddo
endif
do ii = 1,ndict
if (dicxtyp(ii).eq.' ') then
dicxtyp(ii) = 'null'
dictyp(ii) = 'null'
if (tcheck.eq.'yes')
* call warn(' No type specified for name '//
* dicnam(ii)(1:max(1,lastnb(dicnam(ii)))))
endif
enddo
close(dirdev)
nrecd=0
dictfl='no '
500 continue
if(tcheck.eq.'yes') vcheck='yes'
Cdbg WRITE(6,'(i5,3x,a,2x,a)') (i,dicnam(i),dictyp(i),i=1,ndict)
return
end
C
C
C
C
C
C >>>>>> Find position of last non_blank in a string
C
function lastnb(str)
C
integer lastnb
include 'ciftbx.sys'
character*(*) str
integer lenn,ii
lenn = len(str)
do 100 ii=lenn,1,-1
if(str(ii:ii).eq.' ') goto 100
if(str(ii:ii).ne.tab) goto 120
100 continue
ii=1
120 lastnb = ii
return
end
C
C
C
C
C
C >>>>>> Extract the item.category_id from a save frame name
C
subroutine excat(sfname,bcname,lbcname)
C
character*(*) sfname,bcname
integer lbcname,ii,ic,lastnb,lenn
C
C Note that this logic works only for item.category_id
C not for category.id
C
lenn = lastnb(sfname)
bcname = ' '
lbcname = 1
if (lenn.eq.0.or.sfname(1:1).ne.'_') return
do ii = 1,lenn-2
ic = 1+lenn-ii
if (sfname(ic:ic).eq.'.') then
bcname = sfname(2:ic-1)
lbcname = ic-2
return
endif
enddo
return
end
C
C
C
C
C
C >>>>>> Open a CIF and copy its contents into a direct access file.
C
function ocif_(fname)
C
logical ocif_
integer lastnb
include 'ciftbx.sys'
logical test
character fname*(*)
integer case,i
C
save_=.false.
jchar=MAXBUF
lastch=0
if(line_.gt.MAXBUF) call err(' Input line_ value > MAXBUF')
if(nrecd.ne.0) close(dirdev)
nrecd=0
lrecd=0
case=ichar('a')-ichar('A')
tab=char(05)
if(case.lt.0) goto 100
tab=char(09)
bloc_=' '
C
C....... Make sure the CIF is available to open
C
100 file_=fname
do 120 i=1,MAXBUF
if(file_(i:i).eq.' ') goto 140
120 continue
140 longf_=i-1
if (longf_.gt.0) then
inquire(file=file_(1:longf_),exist=test)
ocif_=test
if(.not.ocif_) goto 200
else
file_ = ' '
longf_ = 1
ocif_ = .true.
endif
C
C....... Open up the CIF and a direct access formatted file as scratch
C
if (file_(1:1).ne.' ')
* open(unit=cifdev,file=fname,status='OLD',access='SEQUENTIAL',
* form='FORMATTED')
open(unit=dirdev,status='SCRATCH',access='DIRECT',
* form='FORMATTED',recl=MAXBUF)
C
C....... Copy the CIF to the direct access file
C
160 read(cifdev,'(a,a)',end=180) buffer
nrecd=nrecd+1
irecd=nrecd
if (lastnb(buffer(1:MAXBUF)).gt.line_)
* call warn(' Input line length exceeds line_')
write(dirdev,'(a)',rec=nrecd) buffer
Cdbg WRITE(6,'(i5,1x,a)') nrecd,buffer(1:70)
goto 160
C
180 lrecd=0
jrecd=0
jrect=-1
irecd=0
recn_=0
if (file_(1:1).ne.' ') close(cifdev)
200 return
end
C
C
C
C
C
C >>>>>> Close off direct access file of the current CIF
C and reset all data name tables and pointers
C
subroutine purge_
C
include 'ciftbx.sys'
C
if(nrecd.ne.0) close(dirdev)
recn_=0
save_=.false.
jchar=MAXBUF
lastch=0
nrecd=0
lrecd=0
irecd=0
nname=0
nhash=0
iname=0
loopct=0
loopnl=0
loop_=.false.
text_=.false.
return
end
C
C
C
C
C
C >>>>>> Store the data names and pointers for the requested data block
C
function data_(name)
C
logical data_
integer lastnb
include 'ciftbx.sys'
character name*(*),flag*4,temp*(NUMCHAR),ltype*4
character ctemp*(NUMCHAR)
character locase*(MAXBUF),isbuf*(MAXBUF)
integer ndata,idata,nitem,npakt,i,ii,j,k,kchar,krecd
integer fcatnum,lctemp,isrecd,isjchr,islast
integer pnname,itpos,ipp,ipj
C
jchar=MAXBUF
nname=0
ndata=0
nhash=0
nitem=0
idata=0
iname=0
loopct=0
loopnl=0
ltype=' '
posnam_=0
posval_=0
posdec_=0
posend_=0
data_=.false.
loop_=.false.
text_=.false.
do ii = 1,MAXBOOK
ibkmrk(1,ii)=-1
enddo
irecd=lrecd
lrecd=nrecd
if(name(1:1).ne.' ') irecd=0
call hash_init(dname,dchain,NUMBLOCK,nname,dhash,
* NUMHASH)
call hash_init(cname,cchain,NUMBLOCK,ncname,chash,
* NUMHASH)
isrecd=irecd
isjchr=jchar
islast=lastch
isbuf=' '
if(lastch.gt.0)isbuf(1:lastch)=buffer(1:lastch)
C
C....... Find the requested data block in the file
C
100 call getstr
isjchr=jchar
if(irecd.ne.isrecd) then
isrecd=irecd
islast=lastch
isbuf=' '
if(lastch.gt.0)isbuf(1:lastch)=buffer(1:lastch)
endif
if(type_.eq.'fini') goto 500
if(type_.ne.'text') goto 120
110 call getlin(flag)
if(buffer(1:1).ne.';') goto 110
jchar=2
goto 100
120 continue
if(type_.eq.'save') then
if(long_.lt.6) then
if(.not.save_)
* call err(' Save frame terminator found out of context ')
save_=.false.
goto 100
else
if(save_)
* call err(' Prior save frame not terminated ')
save_=.true.
if(name.eq.' ') goto 150
if(strg_(6:long_).ne.name) goto 100
goto 150
endif
endif
if(type_.ne.'data') goto 100
if(name.eq.' ') goto 150
if(strg_(6:long_).ne.name) goto 100
150 data_=.true.
bloc_=strg_(6:long_)
itpos=jchar-long_
if(tabx_) then
itpos=0
do ipp=1,jchar-long_
itpos=itpos+1
if(buffer(ipp:ipp).eq.tab) itpos=((itpos+7)/8)*8
enddo
endif
posnam_=itpos
C
C....... Get the next token and identify
C
200 call getstr
Cdbg if(dictfl.eq.'no ')
Cdbg * WRITE(6,*) ltype,type_,loop_,nitem,ndata,idata,iname,nname
C
if(ltype.ne.'name') goto 201
if(type_.eq.'numb') goto 203
if(type_.eq.'char') goto 203
if(type_.eq.'text') goto 203
if(type_.eq.'null') goto 203
if(type_.eq.'name'.and.loop_) goto 204
call err(' Illegal tag/value construction')
201 if(ltype.ne.'valu') goto 204
if(type_.eq.'numb') goto 202
if(type_.eq.'char') goto 202
if(type_.eq.'text') goto 202
if(type_.eq.'null') goto 202
goto 204
202 if(nitem.gt.0) goto 205
call err(' Illegal tag/value construction')
203 ltype='valu'
goto 205
204 ltype=type_
C
205 if(type_.eq.'name') goto 206
if(type_.eq.'loop') goto 210
if(type_.eq.'data') goto 210
if(type_.eq.'save') goto 210
if(type_.ne.'fini') goto 220
206 if(loop_) goto 270
210 if(nitem.eq.0) goto 215
C
C....... End of loop detected; save pointers
C
npakt=idata/nitem
if(npakt*nitem.ne.idata) call err(' Item miscount in loop')
loopni(loopct)=nitem
loopnp(loopct)=npakt
nitem=0
idata=0
215 if(type_.eq.'name') goto 270
if(type_.eq.'data') goto 300
if(type_.eq.'save') goto 300
if(type_.eq.'fini') goto 300
C
C....... Loop_ line detected; incr loop block counter
C
loop_=.true.
loopct=loopct+1
if(loopct.gt.NUMLOOP) call err(' Number of loop_s > NUMLOOP')
loorec(loopct)=irecd
loopos(loopct)=jchar-long_
if(quote_.ne.' ') loopos(loopct)=jchar-long_-1
itpos=0
do ipp=1,loopos(loopct)
itpos=itpos+1
if(buffer(ipp:ipp).eq.tab) itpos=((itpos+7)/8)*8
enddo
loopox(loopct)=itpos
goto 200
C
C....... This is the data item; store char position and length
C
220 if(loop_ .and. nitem.eq.0)
* call err(' Illegal tag/value construction')
loop_=.false.
C
i=nname
if(nitem.gt.0) i=i-nitem+mod(idata,nitem)+1
if(i.lt.1) call err(' Illegal tag/value construction')
if(dtype(i).ne.'test') goto 223
if(dictfl.eq.'yes') goto 223
if(tcheck.eq.'no ') goto 223
C>>>> if(long_.eq.1.and.strg_(1:1).eq.'?') goto 223
C>>>> if(long_.eq.1.and.strg_(1:1).eq.'.') goto 223
if(type_.eq.'null') goto 223
if(type_.eq.'numb') goto 223
call warn( ' Numb type violated '//dname(i))
223 if(nitem.le.0) goto 224
idata=idata+1
if(dtype(i).eq.'null') dtype(i)=type_
if(dtype(i).eq.'numb' .and.
* (type_.eq.'char'.or.type_.eq.'text')) dtype(i)='char'
224 if(nname.eq.ndata) goto 230
ndata=ndata+1
if(iloop(ndata).gt.1) goto 225
krecd=irecd
kchar=jchar-long_-1
if(quote_.ne.' ')kchar=kchar-1
225 continue
if(dtype(ndata).eq.' ') dtype(ndata)=type_
drecd(ndata)=krecd
dchar(ndata)=kchar
if(nloop(ndata).gt.0) goto 230
nloop(ndata)=0
iloop(ndata)=long_
C
C....... Skip text lines if present
C
230 if(type_.ne.'text') goto 200
if(nloop(ndata).eq.0) dchar(ndata)=0
if(nloop(ndata).eq.0) iloop(ndata)=long_
250 call getlin(flag)
if(buffer(1:1).eq.';') then
jchar=2
goto 200
endif
if(flag.eq.'fini') call err(' Unexpected end of data')
goto 250
C
C....... This is a data name; store name and loop parameters
C
270 temp=locase(strg_(1:long_))
k=0
if(dictfl.ne.'yes' .and. ndict.gt.0) then
call hash_find(temp,
* dicnam,dicchain,NUMDICT,ndict,dichash,NUMHASH,k)
if(k.ne.0) then
if(alias_ .and. aroot(k).ne.0) temp=dicnam(aroot(k))
endif
endif
pnname=nname
call hash_store(temp,
* dname,dchain,NUMBLOCK,nname,dhash,
* NUMHASH,j)
if(j.eq.pnname+1) then
dtag(j)=strg_(1:long_)
if(k.ne.0) dtag(j)=dictag(k)
trecd(j)=irecd
tchar(j)=jchar-long_
if(quote_.ne.' ') tchar(j)=jchar-long_-1
itpos=0
do ipp=1,tchar(j)
itpos=itpos+1
if(buffer(ipp:ipp).eq.tab) itpos=((itpos+7)/8)*8
enddo
xchar(j)=itpos
endif
if(j.eq.0)
* call err(' Number of data names > NUMBLOCK')
if(k.ne.0)temp=dicnam(k)
if(j.ne.pnname+1) then
call warn(' Duplicate data item '//
* temp(1:max(1,lastnb(temp))))
goto 200
endif
dtype(nname)=' '
dxtyp(nname)=' '
cindex(nname)=0
ddict(nname)=0
ctemp='(none)'
lctemp=6
C
if(dictfl.eq.'yes' .or. vcheck.eq.'no ') goto 290
j=k
if(j.ne.0) then
ddict(nname)=j
cindex(nname)=dcindex(j)
dxtyp(nname)=dicxtyp(j)
dtype(nname)=dictyp(j)
if(vcheck.eq.'no ') goto 280
if(dictyp(j).eq.'numb') then
dtype(nname)='test'
endif
if(cindex(nname).ne.0) then
ctemp=dcname(cindex(nname))
lctemp=lastnb(ctemp)
goto 290
endif
goto 280
endif
call warn(' Data name '//
* temp(1:max(1,lastnb(temp)))
* //' not in dictionary!')
280 call excat(temp,ctemp,lctemp)
if (ctemp.eq.' '.or.'_'//ctemp.eq.temp) then
ctemp = '(none)'
lctemp= 6
if (ndcname.ne.0.and.vcheck.eq.'yes')
* call warn(' No category defined for '
* //temp)
else
call hash_find(ctemp,
* dcname,dcchain,NUMDICT,ndcname,dchash,NUMHASH,j)
if(j.ne.0) then
cindex(nname) = j
else
ipj=ncname
call hash_store(ctemp(1:lctemp),
* cname,cchain,NUMBLOCK,ncname,chash,NUMHASH,j)
if (j.eq.0)
* call err(' Number of categories > NUMBLOCK ')
cindex(nname) = -j
if (ndcname.gt.0.and.j.eq.ipj+1.and.vcheck.eq.'yes')
* call warn(' Category '//
* ctemp(1:lctemp)//' first implicitly defined in cif ')
endif
endif
C
290 lloop(nname)=0
nloop(nname)=0
iloop(nname)=0
if (nitem.eq.0) fcatnum=cindex(nname)
if(.not.loop_) goto 200
nitem=nitem+1
if(nitem.gt.NUMITEM)
* call err(' Items per loop packet > NUMITEM')
nloop(nname)=loopct
iloop(nname)=nitem
if (fcatnum.ne.cindex(nname)) then
temp = '(none)'
if (fcatnum.gt.0) temp=dcname(fcatnum)
if (fcatnum.lt.0) temp=cname(-fcatnum)
if (ctemp(1:lctemp).ne.temp(1:lastnb(temp)))
* call warn (' Heterogeneous categories in loop '//
* ctemp(1:lastnb(ctemp))//' vs '//
* temp(1:lastnb(temp)))
fcatnum=cindex(nname)
endif
goto 200
300 continue
C
C....... Are names checked against dictionary?
C
if(dictfl.eq.'yes') goto 500
if(vcheck.eq.'no '.or.ndict.eq.0) goto 500
do 350 i=1,nname
if(dtype(i).eq.'test') dtype(i)='numb'
350 continue
C
C....... End of data block; tidy up loop storage
C
500 lrecd=irecd-1
if(type_.eq.'save'.and.long_.lt.6) then
itpos=jchar-long_
if(tabx_) then
itpos=0
do ipp=1,jchar-long_
itpos=itpos+1
if(buffer(ipp:ipp).eq.tab) itpos=((itpos+7)/8)*8
enddo
endif
posval_=itpos
endif
irecd=isrecd
jchar=isjchr
lastch=islast
recn_=irecd
buffer=' '
if(lastch.gt.0)buffer=isbuf(1:lastch)
jrecd=irecd
loop_=.false.
loopct=0
if(ndata.ne.nname) call err(' Syntax construction error')
C
Cdbg WRITE(6,'(a)')
Cdbg * ' data name type recd char loop leng'
Cdbg WRITE(6,'(a,1x,a,4i5)') (dname(i),dtype(i),drecd(i),dchar(i),
Cdbg * nloop(i),iloop(i),i=1,nname)
Cdbg WRITE(6,'(3i5)') (i,loopni(i),loopnp(i),i=1,loopct)
C
return
end
C
C
C
C
C
C
C >>>>>> Get the attributes of data item associated with data name
C
function test_(temp)
C
logical test_
include 'ciftbx.sys'
character temp*(*),name*(NUMCHAR)
character otestf*3
character locase*(MAXBUF)
C
otestf=testfl
testfl='yes'
name=locase(temp)
test_=.true.
if(otestf.eq.'no ') goto 100
if(name.eq.nametb) goto 200
100 call getitm(name)
200 list_ =loopnl
if(type_.eq.'null') test_=.false.
return
end
C
C
C
C
C
C >>>>>> Set or Reference a bookmark
C
function bkmrk_(mark)
C
logical bkmrk_
include 'ciftbx.sys'
C
integer mark,ii,nitem
character*4 flag
bkmrk_=.true.
if(mark.eq.0) then
do ii=1,MAXBOOK
if(ibkmrk(1,ii).lt.0) goto 100
enddo
bkmrk_=.false.
call warn(' More than MAXBOOK bookmarks requested')
return
100 mark=ii
ibkmrk(1,ii)=iname
ibkmrk(2,ii)=irecd
ibkmrk(3,ii)=jchar
if(iname.gt.0) then
ibkmrk(2,ii) = trecd(iname)
ibkmrk(3,ii) = tchar(iname)
endif
ibkmrk(4,ii)=0
if(iname.gt.0) then
if(nloop(iname).ne.0.and.
* loopnl.eq.nloop(iname).and.loopct.ne.0) then
nitem=loopni(nloop(iname))
ibkmrk(2,ii)=looprd(1)
ibkmrk(3,ii)=max(0,loopch(1)-1)
ibkmrk(4,ii)=loopct
endif
endif
else
if(ibkmrk(1,mark).lt.0) then
bkmrk_=.false.
return
endif
iname=ibkmrk(1,mark)
irecd=ibkmrk(2,mark)
loopct=ibkmrk(4,mark)
loop_=.false.
text_=.false.
loopnl=-1
testfl='no '
if(iname.gt.0) then
if(nloop(iname).ne.0.and.loopct.ne.0) then
nitem=loopni(nloop(iname))
looprd(nitem+1)=ibkmrk(2,mark)
loopch(nitem+1)=ibkmrk(3,mark)
do ii = 1,nitem
lloop(ii+iname-iloop(iname))=loopct-1
enddo
loopct=loopct-1
if(lloop(iname).gt.0) then
loop_=.true.
loopnl=nloop(iname)
endif
endif
endif
jchar=MAXBUF
if(irecd.gt.0) then
irecd=irecd-1
call getlin(flag)
jchar=ibkmrk(3,mark)
endif
ibkmrk(1,mark)=-1
mark=0
endif
return
end
C
C
C
C
C
C
C >>>>>> Find the location of the requested item in the CIF
C The argument "name" may be a data item name, blank
C for the next such item. The argument "type" may be
C blank for unrestricted acceptance of any non-comment
C string (use cmnt_ to see comments), including loop headers,
C "name" to accept only the name itself and "valu"
C to accept only the value, or "head" to position to the
C head of the CIF. Except when the "head" is requested,
C the position is left after the data item provided.
C
function find_(name,type,strg)
C
logical find_
include 'ciftbx.sys'
character name*(*),type*(*),strg*(*),flag*4
character jjbuf*(MAXBUF)
integer jjchar,jjrecd,jjlast,jjlrec,jjjrec
C
find_ = .false.
strg = ' '
long_ = 0
jjchar = jchar
jjrecd = lrecd
jjlast = lastch
jjlrec = lrecd
jjjrec = jrecd
jjbuf = ' '
if(lastch.gt.0) jjbuf(1:lastch)=buffer(1:lastch)
if(type.eq.'head') then
lrecd = nrecd
irecd=0
jchar=MAXBUF+1
call getlin(flag)
if(flag.eq.'fini') goto 300
find_=.true.
lrecd=jjlrec
return
endif
if(name.ne.' ') then
testfl='no '
call getitm(name)
if(iname.eq.0) goto 300
if(type.eq.'valu') then
list_=loopnl
strg=strg_(1:long_)
find_=.true.
return
endif
if(type.eq.'name'.or.loopnl.eq.0) then
irecd=trecd(iname)-1
call getlin(flag)
jchar=tchar(iname)
posnam_=jchar+1
call getstr
strg=strg_(1:long_)
recn_=irecd
find_=.true.
return
endif
if(type.eq.' ') then
irecd=loorec(loopnl)-1
call getlin(flag)
jchar=loopos(loopnl)
call getstr
posval_=loopos(loopnl)
if(tabx_) posval_=loopox(loopnl)
strg=strg_(1:long_)
recn_=irecd
find_=.true.
return
endif
call err(' Call to find_ with invalid arguments')
endif
if(name.eq.' ') then
200 call getstr
if(type_.eq.'fini') goto 300
if(type.ne.' '.and.
* (type_.eq.'data'.or.type_.eq.'save')) goto 300
if(type.eq.'name'.and.type_.ne.'name') goto 200
if(type.eq.'valu'.and.
* type_.ne.'numb'.and.type_.ne.'text'
* .and.type_.ne.'char'.and.type_.ne.'null') goto 200
find_=.true.
strg=strg_(1:long_)
if(type_.eq.'name') then
posnam_=jchar-long_
else
posval_=jchar-long_
if(quote_.ne.' ') posval_=posval_-1
endif
recn_=irecd
return
endif
C
C Search failed, restore pointers
C
300 irecd = jjrecd
lastch = jjlast
lrecd = jjlrec
jchar = jjchar
buffer = ' '
if(lastch.gt.0)buffer(1:lastch)=jjbuf(1:lastch)
jrecd = jjjrec
if(jrecd.ne.irecd) jrecd=-1
recn_ = irecd
C
return
end
C
C
C
C
C
C
C >>>>>> Get the next data name in the data block
C
function name_(temp)
C
logical name_
include 'ciftbx.sys'
character temp*(*)
C
name_=.false.
temp=' '
iname=iname+1
if(iname.gt.nname) goto 100
name_=.true.
temp=dtag(iname)
if(ddict(iname).ne.0) temp=dictag(ddict(iname))
100 return
end
C
C
C
C
C
C
C >>>>>> Extract a number data item and its standard deviation
C This version return single precision numbers
C
function numb_(temp,numb,sdev)
C
logical numb_
include 'ciftbx.sys'
character temp*(*),name*(NUMCHAR)
character locase*(MAXBUF)
real numb,sdev
C
name=locase(temp)
if(testfl.eq.'no ') goto 100
if(name.eq.nametb) goto 150
C
100 call getitm(name)
C
150 numb_=.false.
if(type_.ne.'numb') goto 200
numb_=.true.
numb =numbtb
if(sdevtb.ge.0.0) sdev=sdevtb
C
200 testfl='no '
return
end
C
C
C
C
C
C
C >>>>>> Extract a number data item and its standard deviation
C This version returns double precision numbers
C
function numd_(temp,numb,sdev)
C
logical numd_
include 'ciftbx.sys'
character temp*(*),name*(NUMCHAR)
character locase*(MAXBUF)
double precision numb,sdev
C
name=locase(temp)
if(testfl.eq.'no ') goto 100
if(name.eq.nametb) goto 150
C
100 call getitm(name)
C
150 numd_=.false.
if(type_.ne.'numb') goto 200
numd_=.true.
numb =numbtb
if(sdevtb.ge.0.0) sdev=sdevtb
C
200 testfl='no '
return
end
C
C
C
C
C
C
C >>>>>> Extract a character data item.
C
function char_(temp,strg)
C
logical char_
include 'ciftbx.sys'
character temp*(*),name*(NUMCHAR)
character strg*(*),flag*4
character locase*(MAXBUF)
integer icpos,itpos,ixpos,ixtpos,ipp,iepos,ispos
C
name=locase(temp)
if(testfl.eq.'yes') goto 100
if(.not.text_) goto 120
if(name.ne.nametb) goto 120
char_=.false.
text_=.false.
strg=' '
long_=0
call getlin(flag)
if(flag.eq.'fini') goto 200
if(buffer(1:1).eq.';') then
jchar=2
goto 200
endif
quote_=' '
jchar=lastch+1
long_=lastch
strg_(1:long_)=buffer(1:long_)
goto 150
C
100 if(name.eq.nametb) goto 150
C
120 call getitm(name)
if(type_.eq.'null') then
char_=.false.
text_=.false.
strg_=' '
long_=0
goto 200
endif
C
150 char_=.true.
text_=.false.
if(tabx_) then
call detab
icpos=jchar-long_
if(quote_.ne.' ') icpos=icpos-1
iepos=icpos+long_-1
itpos=0
do ipp=1,icpos
itpos=itpos+1
if(buffer(ipp:ipp).eq.tab) itpos=((itpos+7)/8)*8
enddo
ispos=itpos
160 ixpos=index(buffer(icpos:iepos),tab)
ixtpos=itpos+ixpos-1
if(ixpos.gt.0.and.ixtpos.le.MAXBUF) then
ixtpos=((ixtpos+7)/8)*8
icpos=icpos+ixpos
itpos=ixtpos+1
if(icpos.le.iepos) goto 160
else
strg =
* bufntb(ispos:min(MAXBUF,itpos+iepos-icpos))
long_=min(MAXBUF,itpos+iepos-icpos)-ispos+1
if(ispos.eq.1.and.strg(1:1).eq.';')
* strg(1:1) = ' '
endif
else
strg=strg_(1:long_)
endif
if(type_.eq.'char') goto 200
char_=.false.
if(type_.ne.'text') goto 200
char_=.true.
call getlin(flag)
jchar=MAXBUF+1
if(flag.eq.'fini') goto 200
if(buffer(1:1).eq.';')then
jchar=2
goto 200
endif
irecd=irecd-1
text_=.true.
C
200 testfl='no '
return
end
C
C
C
C
C
C
C >>>>>> Extract a comment field.
C
function cmnt_(strg)
C
logical cmnt_
integer lastnb
include 'ciftbx.sys'
character strg*(*),flag*4,c*1,
* jjbuf*(MAXBUF)
integer jjchar,jjrecd,jjlast,jjlrec,jjjrec
integer ipp,itpos,ixpos
C
jjchar = jchar
jjrecd = irecd
jjlast = lastch
jjlrec = lrecd
jjjrec = jrecd
jjbuf=' '
if(lastch.gt.0)jjbuf(1:lastch)=buffer(1:lastch)
lrecd = nrecd
if(bloc_.eq.' ') then
if(irecd.eq.0) jchar=MAXBUF
endif
strg=' '
long_=0
cmnt_=.false.
goto 105
100 jchar=jchar+1
105 if(jchar.le.lastch) goto 140
C
C....... Read a new line
C
110 call getlin(flag)
if(flag.eq.'fini') then
strg='fini'
jchar=MAXBUF+1
long_=4
cmnt_=.false.
return
endif
jchar=1
strg=char(0)
long_=1
posnam_=0
goto 220
140 if(lastch.eq.1.and.buffer(1:1).eq.' ') go to 200
C
C....... Process this character in the line
C
150 c=buffer(jchar:jchar)
if(c.eq.' ') goto 100
if(c.eq.tab.and.(.not.tabx_)) goto 190
if(c.eq.tab) goto 100
if(c.eq.'#') goto 200
goto 300
C
C For a tab, when not expanding to blanks, accept
C that single character as a comment
C
190 long_=1
strg=tab
posnam_=jchar
jchar=jchar+1
goto 220
C
C....... Accept the remainder of the line as a comment
C
200 long_=lastch-jchar
itpos=jchar
if(tabx_) then
itpos=0
do ipp=1,jchar
itpos=itpos+1
if(buffer(ipp:ipp).eq.tab) itpos=((itpos+7)/8)*8
enddo
endif
210 posnam_=itpos
if(long_.gt.0) then
if(tabx_) then
call detab
ixpos= lastnb(bufntb)
strg = bufntb(itpos+1:ixpos)
else
strg = buffer(jchar+1:lastch)
endif
endif
if(long_.le.0) then
strg=' '
long_=1
endif
jchar=MAXBUF+1
220 lrecd=jjlrec
cmnt_=.true.
return
C
C....... Found a non-comment field, restore pointers
C
300 irecd = jjrecd
lastch = jjlast
lrecd = jjlrec
jchar = jjchar
buffer=' '
if(lastch.gt.0)buffer(1:lastch)=jjbuf(1:lastch)
jrecd=jjjrec
if(jrecd.ne.irecd) jrecd=-1
recn_=irecd
return
end
C
C
C
C
C
C >>>>> Convert name string to lower case
C
function locase(name)
C
include 'ciftbx.sys'
character locase*(MAXBUF)
character temp*(MAXBUF),name*(*)
character low*26,cap*26,c*1
integer i,j
data cap /'ABCDEFGHIJKLMNOPQRSTUVWXYZ'/
data low /'abcdefghijklmnopqrstuvwxyz'/
C
temp=name
do 100 i=1,MAXBUF
c=temp(i:i)
if(c.eq.' ') goto 200
if(c.eq.tab) goto 200
j=index(cap,c)
if(j.ne.0) temp(i:i)=low(j:j)
100 continue
200 locase=temp
return
end
C
C
C
C
C
C >>>>>> Get the data item associated with the tag.
C
subroutine getitm(name)
C
include 'ciftbx.sys'
SAVE
character name*(*)
character flag*4
integer iitem,nitem,npakt
integer kchar,loopi,i,j,itpos,ipp
C
C....... Find requested dataname in hash list
C
nametb=name
posnam_=0
posval_=0
posdec_=0
posend_=0
quote_=' '
if(name(1:1).eq.'_') goto 100
type_='null'
dictype_='null'
diccat_='(none)'
dicname_=name
tagname_=' '
strg_=' '
long_=1
goto 1000
100 call hash_find(nametb,
* dname,dchain,NUMBLOCK,nname,dhash,NUMHASH,
* iname)
if(iname.gt.0) goto 180
if(dictfl.ne.'yes') then
call hash_find(nametb,
* dicnam,dicchain,NUMDICT,ndict,dichash,NUMHASH,j)
if(j.ne.0) then
dictype_=dicxtyp(j)
if(dcindex(j).ne.0) diccat_=dcname(dcindex(j))
dicname_=nametb
if(aroot(j).ne.0) then
dicname_=dictag(aroot(j))
call hash_find(dicnam(aroot(j)),
* dname,dchain,NUMBLOCK,nname,dhash,NUMHASH,
* iname)
if(iname.gt.0) goto 180
endif
type_='null'
tagname_=' '
strg_=' '
long_=1
go to 1000
endif
endif
160 continue
type_='null'
dictype_='null'
diccat_='(none)'
dicname_=name
long_=1
goto 1000
C
C
180 tagname_=dtag(iname)
if(ddict(iname).ne.0) tagname_=dictag(ddict(iname))
posnam_=tchar(iname)
if(tabx_)posnam_=xchar(iname)
if(nloop(iname).le.0) goto 500
C
C....... Process loop packet if first item request
C
if(nloop(iname).ne.loopnl) goto 200
if(lloop(iname).lt.loopct) goto 300
if(loop_) goto 230
200 loop_=.true.
loopct=0
loopnl=nloop(iname)
nitem=loopni(loopnl)
npakt=loopnp(loopnl)
irecd=drecd(iname)-1
call getlin(flag)
jchar=max(0,dchar(iname)-1)
Cdbg if(jchar.lt.0) write(6,'(7H dchar ,i5)') jchar
do 220 i=1,nitem
220 lloop(i+iname-iloop(iname))=0
goto 240
C
C....... Read a packet of loop items
C
230 nitem=loopni(loopnl)
npakt=loopnp(loopnl)
irecd=looprd(nitem+1)-1
call getlin(flag)
jchar=loopch(nitem+1)
Cdbg if(jchar.lt.0) write(6,'(7H loopch,i5)') jchar
240 iitem=0
250 iitem=iitem+1
if(iitem.le.nitem) goto 255
loopch(iitem)=jchar
looprd(iitem)=irecd
goto 270
255 call getstr
loopch(iitem)=jchar-long_
if(quote_.ne.' ')loopch(iitem)=jchar-long_-1
loopln(iitem)=long_
looprd(iitem)=irecd
if(buffer(1:1).ne.';') goto 250
260 call getlin(flag)
if(buffer(1:1).ne.';') goto 260
jchar=2
goto 250
270 loopct=loopct+1
if(loopct.lt.npakt) goto 300
loop_=.false.
C
C....... Point to the loop data item
C
300 lloop(iname)=lloop(iname)+1
loopi=iloop(iname)
irecd=looprd(loopi)-1
call getlin(flag)
long_=loopln(loopi)
kchar=loopch(loopi)
goto 550
C
C....... Point to the non-loop data item
C
500 irecd=drecd(iname)-1
call getlin(flag)
kchar=dchar(iname)+1
long_=iloop(iname)
loop_=.false.
loopct=0
loopnl=0
C
C....... Place data item into variable string and make number
C
550 type_=dtype(iname)
dictype_=dxtyp(iname)
diccat_='(none)'
if(cindex(iname).gt.0) diccat_=dcname(cindex(iname))
if(cindex(iname).lt.0) diccat_=cname(-cindex(iname))
if(diccat_.eq.' ') diccat_='(none)'
dicname_=dtag(iname)
if(ddict(iname).ne.0) then
if (aroot(ddict(iname)).ne.0) then
dicname_=dictag(aroot(ddict(iname)))
endif
endif
strg_(1:long_)=buffer(kchar:kchar+long_-1)
itpos=kchar
if(tabx_) then
itpos=0
do ipp=1,kchar
itpos=itpos+1
if(buffer(ipp:ipp).eq.tab) itpos=((itpos+7)/8)*8
enddo
endif
posval_=itpos
posend_=itpos+long_-1
jchar=kchar+long_
if(jchar.le.MAXBUF) then
if(buffer(jchar:jchar).ne.' ' .and.
* buffer(jchar:jchar).ne.tab) jchar=jchar+1
endif
quote_=' '
if(kchar.gt.1) then
if(buffer(kchar-1:kchar-1).ne.' ' .and.
* buffer(kchar-1:kchar-1).ne.tab) then
quote_=buffer(kchar-1:kchar-1)
endif
endif
if(type_.eq.'char' .and. kchar.eq.1 .and.
* buffer(1:1).eq.';') type_='text'
if(type_.eq.'text') then
if(buffer(1:1).eq.';') then
strg_(1:1)=' '
else
type_='char'
endif
endif
if(type_.eq.'numb') then
call ctonum
if(posdec_.gt.0) posdec_=posval_+posdec_-1
endif
if(quote_.ne.' ') goto 1000
if(long_.eq.1.and.strg_(1:1).eq.'?') type_='null'
if(long_.eq.1.and.strg_(1:1).eq.'.') type_='null'
C
1000 return
end
C
C
C
C
C
C
C
C >>>>>> Read the next string from the file
C
subroutine getstr
C
include 'ciftbx.sys'
integer i,j,jj(11),im
logical quoted
character c*1,num*21,flag*4
data num/'0123456789+-.()EDQedq'/
C
quoted=.false.
quote_=' '
if(irecd.gt.0.and.
* jchar.eq.1.and.lastch.gt.0) goto 140
100 jchar=jchar+1
if(jchar.le.lastch) goto 150
C
C....... Read a new line
C
110 call getlin(flag)
type_='fini'
dictype_=type_
diccat_='(none)'
dicname_=' '
Cdbg write(6,'(/5i5,a)')
Cdbg * irecd,jrecd,lrecd,nrecd,lastch, buffer(1:lastch)
if(flag.eq.'fini') goto 500
C
C....... Test if the new line is the start of a text sequence
C
140 if(buffer(1:1).ne.';') goto 150
type_='text'
jchar=lastch+1
long_=lastch
strg_(1:long_)=buffer(1:long_)
strg_(1:1)=' '
goto 500
C
C....... Process this character in the line
C
150 c=buffer(jchar:jchar)
if(c.eq.' ') goto 100
if(c.eq.tab) goto 100
if(c.eq.'#') goto 110
if(c.eq.'''') goto 300
if(c.eq.'"') goto 300
if(c.ne.'_') goto 200
type_='name'
goto 210
C
C....... Span blank delimited token; test if a number or a character
C
200 type_='numb'
im=0
do 205 i=1,11
205 jj(i)=0
210 do 250 i=jchar,lastch
if(buffer(i:i).eq.' ') goto 400
if(buffer(i:i).eq.tab) goto 400
if(type_.ne.'numb') goto 250
j=index(num,buffer(i:i))
if(j.eq.0) type_='char'
if(j.le.10) then
im=im+1
goto 250
endif
if(j.gt.13.and.im.eq.0) type_='char'
jj(j-10)=jj(j-10)+1
250 continue
i=lastch+1
if(type_.ne.'numb') goto 400
do 270 j=1,5
if((jj(j).gt.1.and.j.gt.2) .or.
* jj(j).gt.2) type_='char'
270 continue
goto 400
C
C....... Span quote delimited token; assume character
C
300 type_='char'
quoted=.true.
jchar=jchar+1
do 320 i=jchar,lastch
if(buffer(i:i).ne.c) goto 320
if(i+1.ge.lastch) goto 400
if(buffer(i+1:i+1).eq.' ') goto 400
if(buffer(i+1:i+1).eq.tab) goto 400
320 continue
Cdbg write(6,'(a,4i5,a)')
Cdbg * '**** ',irecd,lastch,i,jchar,buffer(jchar:i)
call warn(' Quoted string not closed')
C
C....... Store the string for the getter
C
400 long_=i-jchar
strg_(1:long_)=buffer(jchar:i-1)
jchar=i
quote_=' '
if(quoted) then
quote_=buffer(jchar:jchar)
jchar =jchar+1
endif
Cdbg write(6,'(5x,8i5,5x,a)')
Cdbg * irecd,jrecd,lrecd,nrecd,lastch,i,jchar,long_,strg_(1:long_)
if(type_.ne.'char'.or.quoted) goto 500
if(strg_(1:5).eq.'data_') type_='data'
if(strg_(1:5).eq.'loop_') type_='loop'
if(long_.eq.1.and.strg_(1:1).eq.'?') type_='null'
if(long_.eq.1.and.strg_(1:1).eq.'.') type_='null'
if(strg_(1:5).eq.'save_') type_='save'
C
500 return
end
C
C
C
C
C
C
C >>>>>> Convert a character string into a number and its esd
C
C Q
C D+
C E-
C + +
C number string -xxxx.xxxx-xxx(x)
C component count CCNT 11111222223333444
C (with at least 1 digit in the mantissa)
C
subroutine ctonum
C
integer lastnb
include 'ciftbx.sys'
character test*22,c*1
integer*4 m,nchar
integer*4 ccnt,expn,msin,esin,ndec,ids,nmd
double precision numb,sdev,ntemp,mant
data test /'0123456789+.-()EDQedq '/
C
numbtb=0.D0
sdevtb=-1.D0
numb=1.D0
sdev=0.D0
ccnt=0
mant=0.D0
expn=0.
msin=+1
esin=+1
ndec=0
ids=0
nmd=0
type_='char'
posdec_=0
if(long_.eq.1.and.
* index('0123456789',strg_(1:1)).eq.0) goto 500
C
C....... Loop over the string and identify components
C
C The scan works in phases
C ccnt = 0 processing looking for first digit
C ccnt = 1 processing before decimal point
C ccnt = 2 processing after decimal point
C ccnt = 3 processing exponent
C ccnt = 4 processing standard deviation
C
do 400 nchar=1,long_
C
c=strg_(nchar:nchar)
m=index(test,c)
if(m.eq.0) goto 500
if(m.gt.10) goto 300
C
C....... Process the digits
C
if(ccnt.eq.0) ccnt=1
if(ccnt.eq.2) ndec=ndec+1
if(ccnt.gt.2) goto 220
ntemp=m-1
mant=mant*10.D0+ntemp
nmd=nmd+1
if(ccnt.eq.1.and.mant.ne.0.D0) ids=ids+1
goto 400
220 if(ccnt.gt.3) goto 240
expn=expn*10+m-1
goto 400
240 ntemp=m-1
sdev=sdev*10.D0+ntemp
sdevtb=1.D0
goto 400
C
C....... Process the characters . + - ( ) E D Q
C
300 if(c.ne.'.') goto 320
if(ccnt.gt.1) goto 500
posdec_=nchar
ccnt=2
goto 400
C
320 if(nmd.eq.0.and.m.gt.13) goto 500
if(c.ne.'(') goto 340
if(posdec_.eq.0) posdec_=nchar
ccnt=4
goto 400
C
340 if(posdec_.eq.0.and.ccnt.gt.0) posdec_=nchar
if(c.eq.' ') goto 400
if(m.gt.13) m = 11
if(ccnt.eq.3) goto 500
if(ccnt.gt.0) goto 360
ccnt=1
msin=12-m
goto 400
360 ccnt=3
esin=12-m
C
400 continue
if(posdec_.eq.0) posdec_=lastnb(strg_(1:long_))+1
C
C....... String parsed; construct the numbers
C
expn=expn*esin-ndec
if(expn+ids.gt.-minexp) then
call warn(' Exponent overflow in numeric input')
expn=-minexp-ids
endif
if(expn.lt.minexp) then
call warn(' Exponent underflow in numeric input')
expn=minexp
endif
if(expn.lt.0) numb=1./10.D0**abs(expn)
if(expn.gt.0) numb=10.D0**expn
if(sdevtb.gt.0.0) sdevtb=numb*sdev
ntemp=msin
numbtb=numb*mant*ntemp
type_='numb'
C
500 return
end
C
C
C
C
C
C
C >>>>>> Read a new line from the direct access file
C
subroutine getlin(flag)
C
integer lastnb
include 'ciftbx.sys'
character flag*4
C
irecd=irecd+1
jchar=1
if(irecd.eq.jrecd) goto 200
if(irecd.le.lrecd) goto 100
buffer=' '
lastch=0
jchar=MAXBUF+1
jrecd=-1
flag='fini'
goto 200
100 read(dirdev,'(a)',rec=irecd) buffer
recn_=irecd
lastch=max(1,lastnb(buffer))
jrecd=irecd
flag=' '
200 return
end
C
C
C
C
C
C
C >>>>>> Detab buffer into bufntb
C
subroutine detab
C
include 'ciftbx.sys'
integer icpos,itpos,ixpos,ixtpos
if(jrecd.eq.jrect) return
icpos=1
itpos=1
bufntb=' '
if(lastch.gt.0) then
100 ixpos=index(buffer(icpos:lastch),tab)
ixtpos=ixpos+itpos-1
if(ixpos.gt.0.and.ixtpos.le.MAXBUF) then
ixtpos=((ixtpos+7)/8)*8
if(ixpos.gt.1) then
bufntb(itpos:ixtpos)=
* buffer(icpos:ixpos+icpos-2)
else
bufntb(itpos:ixtpos)=' '
endif
itpos=ixtpos+1
icpos=ixpos+icpos
goto 100
else
bufntb(itpos:max(MAXBUF,itpos+lastch-icpos))=
* buffer(icpos:lastch)
endif
endif
jrect=jrecd
return
end
C
C
C
C
C
C
C >>>>>> Write error message and exit.
C
subroutine err(mess)
character*(*) mess
call cifmsg('error',mess)
stop
end
C
C
C
C
C
C
C >>>>>> Write warning message and continue.
C
subroutine warn(mess)
character*(*) mess
call cifmsg('warning',mess)
return
end
C
C
C
C
C
C
C >>>>>> Write a message to the error device
C
subroutine cifmsg(flag,mess)
C
integer lastnb
include 'ciftbx.sys'
character*(*) flag
character*(*) mess
character*(MAXBUF) tline
character*5 btype
integer ll,ls,ltry,ii,i
C
btype = 'data_'
if(save_) btype = 'save_'
tline= ' ciftbx '//flag//': '
* //file_(1:longf_)//' '//btype
* //bloc_(1:max(1,lastnb(bloc_)))//' line:'
ll = max(1,lastnb(tline))
write(errdev,'(a,i7)')tline(1:ll),irecd
ll=len(mess)
ls=1
100 if(ll-ls.le.79) then
write(errdev,'(1X,a)') mess(ls:ll)
return
else
ltry = min(ll,ls+79)
do ii = ls+1,ltry
i = ltry-ii+ls+1
if(mess(i:i).eq.' ') then
write(errdev,'(1X,a)') mess(ls:i-1)
ls=i+1
if(ls.le.ll) go to 100
return
endif
enddo
write(errdev,'(1X,a)') mess(ls:ltry)
ls=ltry+1
if(ls.le.ll) go to 100
return
endif
end
C
C
C
C
C >>>>>> Create a named file.
C
function pfile_(fname)
C
logical pfile_
include 'ciftbx.sys'
logical test
integer i
character fname*(*)
C
C....... Test if a file by this name is already open.
C
if(pfilef.eq.'yes') call close_
pfilef='no '
file_=fname
do 120 i=1,MAXBUF
if(file_(i:i).eq.' ') goto 140
120 continue
140 if (i.gt.1) then
inquire(file=file_(1:i-1),exist=test)
pfile_=.false.
longf_ = i-1
if(test) goto 200
else
file_ = ' '
pfile_ = .true.
longf_ = 1
endif
C
C....... Open up a new CIF
C
if (file_(1:1) .ne. ' ') then
open(unit=outdev,file=fname,status='NEW',access='SEQUENTIAL',
* form='FORMATTED')
precn_=0
endif
pfile_=.true.
pfilef='yes'
nbloc=0
pchar=1+lprefx
pcharl=0
obuf=prefx
obuf(pchar:MAXBUF)=' '
200 return
end
C
C
C
C
C
C >>>>>> Store a data block command in the CIF
C Call with blank name to close current block only
C
function pdata_(name)
C
logical pdata_
include 'ciftbx.sys'
character name*(*),temp*(MAXBUF)
character dbloc(100)*(NUMCHAR)
integer i
C
pdata_=.true.
if(ploopn.ne.0) call eoloop
if(ptextf.eq.'yes') call eotext
if(psaveo) then
pchar=-1
if(pposval_.ne.0) then
pchar=lprefx+1
call putstr(' ')
pchar=lprefx+pposval_
pposval_=0
endif
call putstr('save_')
psaveo=.false.
endif
C
C....... Check for duplicate data name
C
temp=name
if(temp.eq.' ') goto 200
if(saveo_) goto 130
pdata_=.false.
do 120 i=1,nbloc
if(temp.eq.dbloc(i)) goto 200
120 continue
C
C....... Save block name and put data_ statement
C
nbloc=nbloc+1
if(nbloc.le.100) dbloc(nbloc)=temp
130 pchar=-1
temp='data_'//name
if(saveo_) temp='save_'//name
psaveo=saveo_
if(pposnam_.gt.0) then
pchar=lprefx+1
call putstr(' ')
pchar=lprefx+pposnam_
pposnam_=0
endif
call putstr(temp)
pchar=lprefx
pdata_=.true.
C
200 return
end
C
C
C
C
C
C
C >>>>>> Put a number into the CIF, perhaps with an esd appended
C
function pnumb_(name,numb,sdev)
C
logical pnumb_
include 'ciftbx.sys'
logical flag,tflag
character name*(*),temp*(NUMCHAR)
real numb,sdev
double precision dnumb,dsdev,dprec
C
pnumb_=.true.
flag =.true.
tflag =.true.
temp=name
if(ptextf.eq.'yes') call eotext
C
if(name(1:1).eq.' ') goto 120
if(vcheck.eq.'no ') goto 100
call dcheck(temp,'numb',flag,tflag)
if (aliaso_.and.xdchk.ne.0) then
if (aroot(xdchk).ne.0)
* temp=dictag(aroot(xdchk))
endif
pnumb_=flag
100 if(ploopn.ne.0) call eoloop
pchar=-1
if(pposnam_.ne.0)pchar=pposnam_+lprefx
call putstr(temp)
C
120 if(ploopf.eq.'yes') ploopc=0
ploopf='no '
dprec=decprc
dnumb=numb
dsdev=sdev
call putnum(dnumb,dsdev,dprec)
if(.not.flag) then
if(.not.tabl_) pchar=lprefx+57
call putstr('#< not in dictionary')
endif
if(.not.tflag) then
if(.not.tabl_) pchar=lprefx+57
call putstr('#< not correct type')
endif
C
150 pposnam_=0
pposval_=0
pposdec_=0
pposend_=0
return
end
C
C
C
C
C
C
C >>>>>> Put a double precision number into the CIF, perhaps
C with an esd appended
C
function pnumd_(name,numb,sdev)
C
logical pnumd_
include 'ciftbx.sys'
logical flag,tflag
character name*(*),temp*(NUMCHAR)
double precision numb,sdev
C
pnumd_=.true.
flag =.true.
tflag =.true.
temp=name
if(ptextf.eq.'yes') call eotext
C
if(name(1:1).eq.' ') goto 120
if(vcheck.eq.'no ') goto 100
call dcheck(temp,'numb',flag,tflag)
if (aliaso_.and.xdchk.ne.0) then
if (aroot(xdchk).ne.0)
* temp=dictag(aroot(xdchk))
endif
pnumd_=flag
100 if(ploopn.ne.0) call eoloop
pchar=-1
if(pposnam_.ne.0)pchar=pposnam_+lprefx
call putstr(temp)
C
120 if(ploopf.eq.'yes') ploopc=0
ploopf='no '
call putnum(numb,sdev,dpprc)
if(.not.flag) then
if(.not.tabl_) pchar=lprefx+57
call putstr('#< not in dictionary')
endif
if(.not.tflag) then
if(.not.tabl_) pchar=lprefx+57
call putstr('#< not correct type')
endif
C
150 pposnam_=0
pposval_=0
pposdec_=0
pposend_=0
return
end
C
C
C
C
C
C
C >>>>>> Put a character string into the CIF.
C
function pchar_(name,string)
C
logical pchar_
include 'ciftbx.sys'
logical flag,tflag
character name*(*),temp*(NUMCHAR),string*(*)
character line*(MAXBUF),strg*(MAXBUF)
integer i,j
C
pchar_=.true.
flag =.true.
tflag =.true.
temp =name
if(ptextf.eq.'yes') call eotext
C
if(name(1:1).eq.' ') goto 110
if(vcheck.eq.'no ') goto 100
call dcheck(temp,'char',flag,tflag)
if (aliaso_.and.xdchk.ne.0) then
if (aroot(xdchk).ne.0)
* temp=dictag(aroot(xdchk))
endif
pchar_=flag
100 if(ploopn.ne.0) call eoloop
pchar=-1
if(pposnam_.gt.0) pchar=posnam_+lprefx
call putstr(temp)
C
110 if(ploopf.eq.'yes') ploopc=0
ploopf='no '
line=string
do 120 i=MAXBUF,2,-1
if(line(i:i).ne.' ') goto 130
120 continue
130 if(pposval_.ne.0.and.pposend_.ge.pposval_)
* i=max(i,pposend_-pposval_+1)
if(pquote_.ne.' ') goto 150
do 140 j=i,1,-1
if(line(j:j).eq.' ') goto 150
140 continue
if((line(1:1).eq.'_'
* .or. line(i:i).eq.'_'
* .or. line(1:1).eq.''''
* .or. line(1:1).eq.'"'
* .or. line(1:1).eq.';')
* .and.line(1:i).ne.'''.'''
* .and.line(1:i).ne.'''?'''
* .and.line(1:i).ne.'"."'
* .and.line(1:i).ne.'"?"') goto 150
strg=line(1:i)
goto 200
150 if(pquote_.eq.';') goto 190
if(pquote_.eq.'''') goto 165
if(pquote_.eq.'"') goto 185
do 160 j=1,i
if(line(j:j).eq.'''') goto 170
160 continue
165 strg=''''//line(1:i)//''''
i=i+2
pquote_=''''
goto 200
170 do 180 j=1,i
if(line(j:j).eq.'"') goto 190
180 continue
185 strg='"'//line(1:i)//'"'
i=i+2
pquote_='"'
goto 200
190 pchar=-1
strg='; '//line(1:i)
i=i+2
ptextf='yes'
call putstr(strg(1:i))
pchar=-1
ptextf='no '
call putstr(';')
pchar=lprefx
call putstr(' ')
call warn(' Converted pchar_ output to text for: '
* //strg(3:i))
goto 210
C
200 if(pposval_.ne.0) then
pchar=pposval_+lprefx
if(pquote_.ne.' ') pchar=pchar-1
endif
call putstr(strg(1:i))
210 if(.not.flag) then
if(.not.tabl_) pchar=lprefx+57
call putstr('#< not in dictionary')
endif
if((.not.tflag).and.line(1:i).ne.'.'.and.
* line(1:i).ne.'?'.and.pquote_.eq.' ') then
if(.not.tabl_) pchar=lprefx+57
call putstr('#< not correct type')
endif
250 pposval_=0
pposdec_=0
pposnam_=0
pposend_=0
pquote_=' '
return
end
C
C
C
C
C
C >>>>>> Put a comment in the output CIF
C
function pcmnt_(string)
C
logical pcmnt_
include 'ciftbx.sys'
character string*(*), temp*(MAXBUF)
C
if(ptextf.eq.'yes') call eotext
if(pposnam_.ne.0) pchar=pposnam_+lprefx
if(string.eq.' '.or.
* (string.eq.char(0)) .or.
* (string.eq.tab.and.(.not.ptabx_))) then
if(string.eq.' ') pchar=-1
call putstr(string)
if(string.eq.' ') call putstr(char(0))
else
temp='#'//string
call putstr(temp)
call putstr(char(0))
endif
pcmnt_=.true.
pposnam_=0
if(string.ne.tab)pchar=lprefx+1
return
end
C
C
C
C
C
C
C
C >>>>>> Put a text sequence into the CIF.
C
function ptext_(name,string)
C
logical ptext_
integer lastnb
include 'ciftbx.sys'
logical flag,tflag
integer ll
character name*(*),temp*(NUMCHAR),string*(*),store*(NUMCHAR)
character temp2*(MAXBUF)
data store/' '/
C
ptext_=.true.
flag =.true.
tflag =.true.
ll=lastnb(string)
if(ploopf.eq.'yes') ploopc=0
ploopf='no '
temp=name
if(ptextf.eq.'no ') goto 100
if(temp.eq.store) goto 150
call eotext
C
100 if(name(1:1).ne.' ') goto 110
if(ptextf.eq.'yes') goto 150
goto 130
C
110 if(ploopn.ne.0) call eoloop
if(vcheck.eq.'no ') goto 120
call dcheck(name,'char',flag,tflag)
if (aliaso_.and.xdchk.ne.0) then
if (aroot(xdchk).ne.0)
* temp=dictag(aroot(xdchk))
endif
ptext_=flag
120 pchar=-1
if(pposnam_.ne.0) pchar=pposnam_+lprefx
call putstr(temp)
if(.not.flag) then
if(.not.tabl_) pchar=lprefx+57
call putstr('#< not in dictionary')
endif
if(.not.tflag) then
if(.not.tabl_) pchar=lprefx+57
call putstr('#< not correct type')
endif
130 ptextf='yes'
store=temp
if(string(1:1).eq.' '.and.ll.gt.1) then
pchar=-1
temp2=';'//string(2:ll)
call putstr(temp2)
pchar=-1
return
endif
pchar=-1
call putstr(';')
pchar=-1
if(string.eq.' ') return
150 pchar=-1
call putstr(string(1:max(1,ll)))
pchar=-1
pposnam_=0
pposval_=0
pposdec_=0
pposend_=0
return
end
C
C
C
C
C
C
C >>>>>> Put a loop_ data name into the CIF.
C
function ploop_(name)
C
logical ploop_
include 'ciftbx.sys'
logical flag,tflag
character name*(*),temp*(NUMCHAR)
C
ploop_=.true.
flag =.true.
if(ptextf.eq.'yes') call eotext
if(ploopf.eq.'no ') call eoloop
temp=' '
if(name(1:1).eq.' ') goto 100
C
if(tabl_.and.pposnam_.eq.0) then
temp=' '//name
else
temp=name
endif
if(vcheck.eq.'no ') goto 100
call dcheck(name,' ',flag,tflag)
if (aliaso_.and.xdchk.ne.0) then
if (aroot(xdchk).ne.0) then
if(tabl_.and.pposnam_.eq.0) then
temp=' '//dictag(aroot(xdchk))
else
temp=dictag(aroot(xdchk))
endif
endif
endif
ploop_=flag
100 if(ploopn.ne.0) goto 120
ploopf='yes'
pchar=-1
if(pposval_.ne.0) then
pchar=lprefx+1
call putstr(' ')
pchar=pposval_+lprefx
else
if(pposnam_.ne.0) then
pchar=lprefx+1
call putstr(' ')
pchar=pposnam_+lprefx+1
endif
endif
call putstr('loop_')
pchar=-1
if(name(1:1).eq.' ') then
ploopn=-1
return
endif
120 if(pposnam_.ne.0) pchar=pposnam_+lprefx
call putstr(temp)
if(flag) goto 130
if(.not.tabl_) pchar=lprefx+57
call putstr('#< not in dictionary')
130 pchar=lprefx+1
ploopn=max(ploopn,0)+1
C
150 return
end
C
C
C
C
C
C >>>>>> Create or clear a prefix string
C Any change in the length of the prefix string flushes
C pending text, if any, loops and partial output lines
C
function prefx_(strg,lstrg)
C
logical prefx_
include 'ciftbx.sys'
character strg*(*)
integer lstrg,mxline
C
mxline=MAXBUF
if(line_.gt.0) mxline=min(line_,MAXBUF)
if(lstrg.ne.lprefx.and.pcharl.gt.0) then
pchar=-1
call putstr(' ')
endif
if (lstrg.le.0) then
prefx=' '
if(pchar.ge.lprefx+1)pchar=pchar-lprefx
lprefx=0
else
if(lstrg.gt.mxline) then
call warn(' Prefix string truncated')
endif
prefx=strg
if(pchar.ge.lprefx+1)pchar=pchar-lprefx+lstrg
obuf(1:min(mxline,lstrg))=prefx
lprefx=lstrg
if(mxline-lprefx.lt.NUMCHAR) then
call warn(' Output prefix may force line overflow')
endif
endif
prefx_=.true.
return
end
C
C
C
C
C
C
C >>>>>> Close the CIF
C
subroutine close_
C
include 'ciftbx.sys'
C
if(ptextf.eq.'yes') call eotext
if(ploopn.ne.0) call eoloop
if(pcharl.ge.lprefx+1) then
pchar=-1
call putstr(' ')
endif
if (file_(1:1) .ne. ' ') then
close(outdev)
precn_=0
endif
return
end
C
C
C
C
C
C >>>>>> Put the string into the output CIF buffer
C
subroutine putstr(string)
C
integer lastnb
include 'ciftbx.sys'
SAVE
character string*(*),temp*(MAXBUF),bfill*(MAXBUF)
character temp2*(MAXBUF)
integer i,ii,mxline,ioffst,ifree,icpos,itpos
integer ixpos,ixtpos,it,im,kbin,kpass
logical pflush,waslop
data waslop /.false./
C
bfill = ' '
mxline=MAXBUF
if(line_.gt.0) mxline=min(line_,MAXBUF)
temp=string
temp2=temp
pflush=.false.
if(pchar.lt.0) pflush=.true.
C
do 100 i=MAXBUF,1,-1
if(temp(i:i).eq.' ') goto 100
if(ptabx_.and.temp(i:i).eq.tab) goto 100
goto 110
100 continue
i=0
it=i
C
C....... Organise the output of loop_ items
C
110 if(i.eq.0) goto 130
if(i.eq.1.and.string.eq.tab) goto 130
if(i.eq.1.and.string.eq.char(0)) then
pcharl=MAXBUF
goto 200
endif
if(temp(1:1).eq.'#') goto 130
if(ploopf.eq.'yes') goto 130
if(ptextf.eq.'yes') goto 130
if(ploopn.le.0) goto 130
ploopc=ploopc+1
if(align_.or.tabl_) then
if(ploopc.gt.ploopn) then
if(pcharl.gt.lprefx) pflush=.true.
ploopc=1
if(pchar.gt.0) pchar=lprefx+1
endif
if(pchar.lt.0) goto 130
if(tabl_) then
kbin=(mxline-lprefx)/8
if(ploopn.lt.kbin) then
if(kbin/(ploopn+1).gt.1) then
pchar=9+lprefx+
* (ploopc-1)*8*(kbin/(ploopn+1))
else
pchar=1+lprefx+
* (ploopc-1)*8*(kbin/ploopn)
endif
else
if(ploopc.le.kbin) then
pchar=1+lprefx+(ploopc-1)*8
else
kpass=(ploopc-kbin-1)/(kbin-1)+1
pchar=2*kpass+1+lprefx+
* mod(ploopc-kbin-1,kbin-1)*8
endif
endif
else
if(ptabx_) then
icpos=1
itpos=1
120 ixpos=index(temp(icpos:i),tab)
ixtpos=(pchar+itpos-1+ixpos)
ixtpos=((ixtpos+7)/8)*8
if(ixpos.gt.0) then
if(ixpos.gt.1) then
temp2(itpos:ixtpos-pchar+1)=temp(icpos:ixpos-1)
else
temp2(itpos:ixtpos-pchar+1)=' '
endif
icpos=ixpos+1
itpos=ixtpos+2-pchar
if(icpos.le.i) goto 120
it=itpos-1
else
temp2(itpos:itpos+i-icpos)=temp(icpos:i)
it=itpos+i-icpos
endif
endif
if((pchar+i).gt.mxline+1.or.
* (ptabx_.and.pchar+it.gt.mxline+1)) then
if(pcharl.gt.lprefx)pflush=.true.
pchar=lprefx+1
endif
endif
else
if(ploopc.le.ploopn) goto 130
ploopc=1
endif
C
C....... Is the buffer full and needs flushing?
C
130 if(i.eq.1.and.string.eq.tab) then
if(pcharl.gt.lprefx) then
if(obuf(pcharl:pcharl).eq.' ') pcharl=pcharl-1
endif
endif
if(pchar.le.pcharl.and.pcharl.gt.lprefx) pflush=.true.
pchar=max(lprefx+1,pchar)
if((ploopf.eq.'yes'.or.ploopn.le.0).and.tabl_)
* pchar=((pchar-lprefx+6)/8)*8+1+lprefx
if(ptabx_) then
icpos=1
itpos=1
135 ixpos=index(temp(icpos:i),tab)
ixtpos=(pchar+itpos-1+ixpos)
ixtpos=((ixtpos+7)/8)*8
if(ixpos.gt.0) then
if(ixpos.gt.1) then
temp2(itpos:ixtpos-pchar+1)=temp(icpos:ixpos-1)
else
temp2(itpos:ixtpos-pchar+1)=' '
endif
icpos=ixpos+1
itpos=ixtpos+2-pchar
if(icpos.le.i) goto 135
it=itpos-1
else
temp2(itpos:itpos+i-icpos)=temp(icpos:i)
it=itpos+i-icpos
endif
endif
if((pchar+i).gt.mxline+1.or.
* (ptabx_.and.pchar+it.gt.mxline+1)) then
pflush=.true.
pchar=mxline+1-i
pchar=max(lprefx+1,pchar)
endif
if(.not.pflush) goto 150
140 if(pcharl.gt.lprefx) then
if(waslop.or.(.not.tabl_)) goto 145
ioffst=0
pcharl=max(lastnb(obuf(1:pcharl)),lprefx+1)
ifree=mxline-pcharl
if(ifree.gt.0) then
im=numtab+2
if(numtab.gt.0.and.numtab.le.MAXTAB) then
if(obuf(itabp(numtab):itabp(numtab)).eq.'#')
* im=im-1
endif
if(ifree.ge.16.and.im.lt.4.and.
* (obuf(1+lprefx:1+lprefx).ne.'#'
* .and.obuf(1+lprefx:1+lprefx).ne.';'
* .and.obuf(1+lprefx:1+lprefx).ne.'_'
* .and.obuf(1+lprefx:1+lprefx).ne.' '
* .and.obuf(1+lprefx:5+lprefx).ne.'data_'
* .and.obuf(1+lprefx:5+lprefx).ne.'save_'
* .and.obuf(1+lprefx:5).ne.'loop_')) then
temp(1+lprefx:pcharl)=obuf(1+lprefx:pcharl)
obuf(1+lprefx:pcharl+8)=
* bfill(1:8)//temp(1+lprefx:pcharl)
ioffst = 8
ifree=ifree-8
pcharl=pcharl+8
endif
do ii=1,min(MAXTAB,numtab)
icpos=itabp(ii)+ioffst
if(icpos.gt.pcharl) goto 145
if(im.lt.4) then
itpos=(max(icpos-lprefx,
* ii*(mxline-lprefx)/im)+6)/8
itpos=itpos*8+1+lprefx
else
itpos=(max(icpos-lprefx,
* ii*(mxline-lprefx)/im)+4)/6
itpos=itpos*6+1+lprefx
endif
if((obuf(icpos:icpos).eq.''''.or.
* obuf(icpos:icpos).eq.'"').and.
* itpos.gt.icpos) itpos=itpos-1
if(itpos-icpos.gt.ifree) itpos=icpos+ifree
if(itpos.gt.icpos) then
temp(1:pcharl-icpos+1)=
* obuf(icpos:pcharl)
if(i.lt.numtab) then
ixpos=itabp(ii+1)+ioffst
if(ixpos.gt.icpos+itpos-icpos+1) then
if(obuf(ixpos-(itpos-icpos+1):ixpos-1).eq.
* bfill(1:itpos-icpos+1)) then
temp(ixpos-itpos+1:pcharl-itpos+1)=
* obuf(ixpos:pcharl)
pcharl=pcharl-(itpos-icpos)
endif
endif
endif
obuf(icpos:pcharl+itpos-icpos)=
* bfill(1:itpos-icpos)//temp(1:pcharl-icpos+1)
ifree=ifree-(itpos-icpos)
ioffst=ioffst+itpos-icpos
pcharl=pcharl+itpos-icpos
endif
if(ifree.le.0) goto 145
enddo
endif
145 pcharl=max(1,lastnb(obuf))
write(outdev,'(a)') obuf(1:pcharl)
else
if(precn_.gt.0) then
if(lprefx.gt.0) then
write(outdev,'(a)') obuf(1:lprefx)
else
write(outdev,'(a)')
endif
else
precn_=precn_-1
endif
endif
waslop=.false.
precn_=precn_+1
do ii = 1,MAXTAB
itabp(ii)=0
enddo
numtab=0
if(lprefx.gt.0) then
obuf=prefx(1:lprefx)
else
obuf=' '
endif
C
C....... Load the next item into the buffer
C
150 pcharl=pchar+i
if(ptabx_) pcharl=pchar+it
waslop= ploopf.eq.'no '.and.ploopn.gt.0.and.align_
if(i.eq.0) then
if(pcharl.eq.lprefx+1.and.
* obuf(lprefx+1:lprefx+1).eq.' ') pcharl=pcharl-1
pchar=pcharl+1
goto 200
endif
if(ptabx_) then
obuf(pchar:pcharl)=temp2(1:it)
else
if(string.eq.tab) pcharl=pcharl-1
obuf(pchar:pcharl)=string(1:i)
endif
if(pchar.gt.1+lprefx) then
numtab=numtab+1
if(numtab.le.MAXTAB) itabp(numtab)=pchar
endif
pchar=pcharl+1
if(pchar.gt.mxline+2) then
call warn(' Output CIF line longer than line_')
endif
C
200 return
end
C
C
C
C
C
C >>>>>> Convert the number and esd to string nnnn(m), limited
C by relative precision prec
C
subroutine putnum(numb,sdev,prec)
C
include 'ciftbx.sys'
character string*30,temp*30,c*1,sfmt*8
double precision numb,sdev,prec,xxnumb,xsdev,slog
integer i,iexp,ifp,ii,jj,j,jlnz,jn,kexp,m,ixsdev,islog
integer kdecp,ibexp
C
kdecp=0
if (sdev.gt.abs(numb)*prec) then
if (esdlim_.ne.esdcac) then
C
C determine the number of digits set by esdlim_
C
if (esdlim_.lt.9 .or.esdlim_.gt.99999) then
call warn(' Invalid value of esdlim_ reset to 19')
esdlim_ = 19
endif
C
C determine the number of esd digits
C
esddig = 1.+alog10(float(esdlim_))
esdcac = esdlim_
endif
C
C determine kexp, the power of 10 necessary
C to present sdev as an integer in the range
C (esdlim_/10,esdlim_]
C
slog = dlog10(sdev)
islog = slog+1000.
islog = islog-1000
kexp = -islog+esddig
C
C Adjust exponent kexp, so that sdev*10**kexp
C is in the interval (esdlim_/10,esdlim_]
C
20 if (kexp.lt.minexp) then
call warn(' Underflow of esd')
ixsdev = 0
go to 30
endif
if (kexp.gt.-minexp) then
call warn(' Overflow of esd')
ixsdev = 99999
go to 30
endif
xsdev = sdev*10.D0**kexp
ixsdev = xsdev+.5
if (ixsdev.gt.esdlim_) then
kexp = kexp -1
go to 20
endif
if (ixsdev.lt.(esdlim_+5)/10) then
kexp = kexp+1
go to 20
endif
C
C We need to present the number to the same scaling
C at first, but will adjust to avoid Ennn notation
C if possible
C
30 xxnumb = dabs(numb)*10.d0**kexp+.5
if(xxnumb*prec .gt.1.D0) then
call warn(' ESD less than precision of machine')
ixsdev=0
endif
if(numb.lt.0.d0) xxnumb = -xxnumb
write(string,ndpfmt)xxnumb
C
C Extract the power of 10
C
iexp = 0
ibexp = 0
do ii = 0,4
i = 30-ii
c = string(i:i)
m = index('0123456789',c)
if (m.gt.0) then
iexp = iexp+(m-1)*10**(ii-ibexp)
else
if (c.eq.' ') then
ibexp = ibexp+1
else
if (c.eq.'-') iexp=-iexp
goto 40
endif
endif
enddo
call err(' Internal error in putnum')
C
C Scan the rest of the string shifting the
C decimal point to get an integer
C
40 ifp = 0
j=1
do ii = 1,i-1
c = string(ii:ii)
if (c.ne.' ')then
m=index('0123456789+-',c)
if(m.ne.0) then
temp(j:j)=c
if(j.gt.1.or.c.ne.'0')j=j+1
if(j.eq.3.and.temp(1:2).eq.'-0')j=j-1
if(ifp.ne.0)then
iexp=iexp-1
if(iexp.le.0) goto 50
endif
else
if(c.eq.'.') then
ifp=1
if(iexp.le.0) goto 50
endif
endif
endif
enddo
C
C The string from 1 to j-1 has an integer
C If iexp < 0, we present a 0. If iexp > 0
C we pad with zeros
C
50 if(j.eq.1.or.iexp.lt.0) then
temp(1:1)='0'
j=2
iexp = 0
endif
if (iexp.gt.0) then
do ii = 1,iexp
temp(j:j)='0'
j=j+1
enddo
iexp=0
endif
string=temp(1:j-1)
C
C We have the number for which the presentation
C would be nnnnnE-kexp. If kexp is gt 0, we can
C decrease it and introduce a decimal point
C
jj=0
if(index('0123456789',temp(1:1)).eq.0) jj=1
if(kexp.gt.0.and.kexp.lt.j-jj+8) then
if(kexp.lt.j-1) then
string=temp(1:j-1-kexp)//'.'//
* temp(j-kexp:j-1)
kexp = 0
j=j+1
else
if(jj.ne.0)string(1:1)=temp(1:1)
string(1+jj:1+jj)='.'
do ii=1,kexp-(j-1-jj)
string(1+jj+ii:1+jj+ii)='0'
enddo
string(2+jj+(kexp-(j-1-jj)):30)=
* temp(1+jj:j-1)
j=j+1+kexp-(j-1-jj)
kexp=0
endif
endif
kdecp=index(string(1:j-1),'.')
if(kdecp.eq.0) kdecp=j
if(kexp.ne.0) then
write(temp(1:5),'(i5)') -kexp
string(j:j)='E'
j=j+1
do ii=1,5
c=temp(ii:ii)
if(c.ne.' ') then
string(j:j)=c
j=j+1
endif
enddo
endif
C
C if there is a standard deviation
C append it in parentheses
C
if(ixsdev.ne.0) then
write(temp(1:5),'(i5)') ixsdev
string(j:j)='('
j=j+1
do ii=1,5
c=temp(ii:ii)
if(c.ne.' ') then
string(j:j)=c
j=j+1
endif
enddo
string(j:j)=')'
j=j+1
endif
else
C
C There is no standard deviation, just write numb
C But limit to the digits implied by prec
C
slog = dlog10(min(.1D0,max(prec,dpprc)))
islog = slog+1000.5
islog = islog-1000
kexp = -islog
write(sfmt,'(5h(D30.,i2,1h))') kexp
write(temp,sfmt)numb
C
C Now have the number in the form
C [sign][0].nnnnnnnnDeee
C which, while sufficient, is not neat
C we reformat for the case 0<=eee<=kexp
C
C
C Extract the power of 10
C
iexp = 0
ibexp = 0
do ii = 0,4
i = 30-ii
c = temp(i:i)
m = index('0123456789',c)
if (m.gt.0) then
iexp = iexp+(m-1)*10**(ii-ibexp)
else
if (c.eq.' ') then
ibexp = ibexp+1
else
if (c.eq.'-') iexp=-iexp
goto 140
endif
endif
enddo
call err(' Internal error in putnum')
C
C Scan the rest of the string shifting the
C decimal point to get a number with exponent 0,
C if possible
C
140 ifp = 0
j=1
do ii = 1,i-1
jn=ii
c = temp(ii:ii)
if (c.ne.' ')then
m=index('0123456789+-',c)
if(m.ne.0) then
string(j:j)=c
if(j.gt.1.or.c.ne.'0')j=j+1
if(j.eq.3.and.string(1:2).eq.'-0')j=j-1
if(ifp.ne.0)then
iexp=iexp-1
if(iexp.le.0) goto 150
endif
else
if(c.eq.'.') then
ifp = -1
if(iexp.le.0) goto 150
endif
endif
endif
enddo
150 string(j:j)='.'
ifp = j
j = j+1
jlnz = j-1
155 do ii = jn+1,i-1
c = temp(ii:ii)
if (c.ne.' ')then
m=index('0123456789',c)
if(m.ne.0) then
string(j:j)=c
j=j+1
if(m.ne.1)jlnz=j
if(m.eq.1.and.ifp.ge.1.and.
* pposdec_.ne.0.and.pposend_.ne.0) then
if(j-1-ifp-min(iexp,0).le.pposend_-pposdec_)
* jlnz=j
endif
else
goto 160
endif
endif
enddo
160 j=jlnz
if(j.eq.1) then
string(1:1)='0'
j=2
endif
if(iexp.lt.0.and.iexp.gt.-7.and.ifp.lt.j-1.and.
* ifp.ne.0.and.j-ifp-iexp.le.kexp) then
temp(1:ifp)=string(1:ifp)
do ii = 1,-iexp
temp(ifp+ii:ifp+ii) = '0'
enddo
temp(ifp-iexp+1:j-iexp-1) = string(ifp+1:j-1)
j = j-iexp
iexp=0
string(1:j-1) = temp(1:j-1)
endif
kdecp=index(string(1:j-1),'.')
if(kdecp.eq.0) kdecp=j
if(iexp.ne.0) then
write(temp(1:5),'(i5)')iexp
string(j:j)='E'
j=j+1
do ii=1,5
c=temp(ii:ii)
if(c.ne.' ') then
string(j:j)=c
j=j+1
endif
enddo
endif
endif
C
if(j.lt.1) then
string(1:1)='0'
j=2
endif
if(kdecp.lt.1)kdecp=j
if(pposdec_.ne.0) then
pchar=lprefx+pposdec_-kdecp+1
else
if(pposval_.ne.0)pchar=lprefx+pposval_
endif
call putstr(string(1:j-1))
return
end
C
C
C
C
C
C >>>>>> Check dictionary for data name validation
C
subroutine dcheck(name,type,flag,tflag)
C
include 'ciftbx.sys'
logical flag,tflag
character name*(*),temp*(NUMCHAR),
* locase*(MAXBUF),type*4
C
flag=.true.
tflag=.true.
temp=locase(name)
call hash_find(temp,
* dicnam,dicchain,NUMDICT,ndict,dichash,NUMHASH,xdchk)
if(xdchk.eq.0) goto 150
if(tcheck.eq.'no ') goto 200
if(type.eq.dictyp(xdchk)) goto 200
if(type.eq.' ') goto 200
if(dictyp(xdchk).eq.'text' .and. type.eq.'char') goto 200
tflag=.false.
goto 200
150 flag=.false.
200 continue
return
end
C
C
C
C
C
C >>>>>> End of text string
C
subroutine eotext
C
include 'ciftbx.sys'
C
if(ptextf.ne.'yes') then
call warn(' Out-of-sequence call to end text block')
return
endif
ptextf='no '
pchar=-1
call putstr(';')
call putstr(char(0))
return
end
C
C
C
C
C
C >>>>>> End of loop detected; check integrity and tidy up pointers
C
subroutine eoloop
C
include 'ciftbx.sys'
integer i
C
if(ploopn.eq.0) goto 200
if(ploopn.eq.-1) then
call putstr('_DUMMY')
ploopn=1
ploopc=0
call warn(
* ' Missing: missing loop_ name set as _DUMMY')
endif
if(ploopn.eq.ploopc) goto 200
do 150 i=ploopc+1,ploopn
150 call putstr('DUMMY')
call warn(
* ' Missing: missing loop_ items set as DUMMY')
C
200 ploopc=0
ploopn=0
return
end
C
C
C
C
C
C
C >>>>>> Set common default values
C
block data
C
include 'ciftbx.sys'
data cifdev /1/
data outdev /2/
data dirdev /3/
data errdev /6/
data loopct /0/
data nhash /0/
data ndict /0/
data nname /0/
data nbloc /0/
data ploopn /0/
data ploopc /0/
data ploopf /'no '/
data ptextf /'no '/
data pfilef /'no '/
data testfl /'no '/
data vcheck /'no '/
data tcheck /'no '/
data align_ /.true./
data tabl_ /.true./
data tabx_ /.true./
data ptabx_ /.true./
data text_ /.false./
data loop_ /.false./
data ndcname /0/
data ncname /0/
data save_ /.false./
data saveo_ /.false./
data alias_ /.true./
data aliaso_ /.false./
data dchash /NUMHASH*0/
data dichash /NUMHASH*0/
data dhash /NUMHASH*0/
data dcchain /NUMDICT*0/
data aroot /NUMDICT*0/
data cindex /NUMBLOCK*0/
data line_ /80/
data lastch /0/
data dictype_ /' '/
data dicname_ /' '/
data diccat_ /' '/
data tagname_ /' '/
data prefx /' '/
data lprefx /0/
data esdlim_ /19/
data esdcac /19/
data esddig /2/
data esdfmt /'(e12.2)'/
data edpfmt /'(d12.2)'/
data ndpfmt /'(d30.14)'/
data decprc /1.e-6/
data dpprc /1.d-14/
data decmin /1.e-37/
data dpmin /1.d-307/
data minexp /-307/
data itabp /MAXTAB*0/
data jrect /-1/
data numtab /0/
data recn_ /0/
data precn_ /0/
data posnam_ /0/
data posval_ /0/
data posdec_ /0/
data posend_ /0/
data pposnam_ /0/
data pposval_ /0/
data pposdec_ /0/
data pposend_ /0/
data quote_ /' '/
data pquote_ /' '/
data ibkmrk /MAXBOOK*-1,MAXBOOK*-1,
* MAXBOOK*-1,MAXBOOK*-1/
end
C
C
C change the following include to include 'clearfp_sun.f'
C for use on a SUN
C
include 'clearfp.f'
"CUT_HERE_CUT_HERE_CUT_HERE"
#___________ THE FILE 'ciftbx.f' ENDS TWO LINES BEFORE THIS ___________
if ( 125819 != `wc -c <'ciftbx.f'` ) then
echo "cshar: 'ciftbx.f' unpacked with wrong size"
endif
# end of 'ciftbx.f'
endif
if ( -e 'ciftbx.sys' && ${1} != "-c" ) then
echo "cshar: Will not clobber existing file 'ciftbx.sys'"
sed "s/^X//" >'/dev/null' <<"CUT_HERE_CUT_HERE_CUT_HERE"
else
echo "cshar: Extracting 'ciftbx.sys' (12073 characters)"
#___________ THE FILE 'ciftbx.sys' BEGINS TWO LINES AFTER THIS __________
sed "s/^X//" >'ciftbx.sys' <<"CUT_HERE_CUT_HERE_CUT_HERE"
C
C
C >>>>>> Common declararations 'ciftbx.sys'
C
C For inclusion in 'ciftbx.f'
C
C Include all user common definitions, but not functions
include 'ciftbx.cmv'
C
C Define parameters controlling the sizes of things
C
integer NUMDICT,NUMHASH,NUMBLOCK,NUMLOOP,NUMITEM,NUMIP1,MAXTAB,
* MAXBOOK
C
C Number of entries in dictionary tables
PARAMETER (NUMDICT=2500)
C Number of hash table entries (a modest prime)
PARAMETER (NUMHASH=53)
C Number of entries in data block tables
PARAMETER (NUMBLOCK=500)
C Number of loops in a data block
PARAMETER (NUMLOOP=50)
C Number of items in a loop
PARAMETER (NUMITEM=50)
C Maximum number of tabs in output cif line
PARAMETER (MAXTAB=10)
C Maximum number of simultaneous bookmarks
PARAMETER (MAXBOOK=1000)
C
C Derived value for loop arrays
PARAMETER (NUMIP1=NUMITEM+1)
C
C Flag if test_ last called (yes/no )
character testfl*3
C Name of current data item
character nametb*(NUMCHAR)
C Tab character for this machine
character tab*1
C Character buffer for reading lines
character buffer*(MAXBUF)
C Character buffer for tab-expanded lines
character bufntb*(MAXBUF)
C Dictionary validation check
character vcheck*3
C Dictionary data type check
character tcheck*3
C Dictionary flag
character dictfl*3
C Dictionary names
character dicnam(NUMDICT)*(NUMCHAR)
C Dictionary names with upper/lower case
character dictag(NUMDICT)*(NUMCHAR)
C Dictionary data types
character dictyp(NUMDICT)*4
C Data names in data block
character dname(NUMBLOCK)*(NUMCHAR)
C Data names in data block with u/l case
character dtag(NUMBLOCK)*(NUMCHAR)
C Data type of data item
character dtype(NUMBLOCK)*4
C Flag signalling output CIF open
character pfilef*3
C Flag signalling loop_ being loaded
character ploopf*3
C Flag signalling text being loaded
character ptextf*3
C Output string buffer
character obuf*(MAXBUF)
C Output prefix string
character prefx*(MAXBUF)
C Dictionary category names
character dcname(NUMDICT)*(NUMCHAR)
C Category names in use
character cname(NUMBLOCK)*(NUMCHAR)
C Dictionary extended types
character dicxtyp(NUMDICT)*(NUMCHAR)
C Datablock extended types
character dxtyp(NUMBLOCK)*(NUMCHAR)
C Format for single precision esd's
character*13 esdfmt
C Format for double precision esd's
character*13 edpfmt
C Format for writing double precision numbers
character*8 ndpfmt
C Record number containing data item
integer drecd(NUMBLOCK)
C Character position of item in record
integer dchar(NUMBLOCK)
C Record number containing tag of item
integer trecd(NUMBLOCK)
C Character position of tag of item in record
integer tchar(NUMBLOCK)
C Character position of tag of item in record
C tab-expanded
integer xchar(NUMBLOCK)
C Loop block number (0 for non-loop)
integer nloop(NUMBLOCK)
C Item count in loop packet
integer iloop(NUMBLOCK)
C Loop line counter; initially zero
integer lloop(NUMBLOCK)
C Number of items per packet in each loop
integer loopni(NUMLOOP)
C Number of packets per loop
integer loopnp(NUMLOOP)
C Record number of loop header
integer loorec(NUMLOOP)
C Character position of loop header
integer loopos(NUMLOOP)
C Character position of loop header
C tab-expanded
integer loopox(NUMLOOP)
C Number of last non-blank chars in input
integer lastch
C Number of items in current loop packet
integer ploopc
C Number of items in output loop packet
integer ploopn
C Number of current loop block
integer loopnl
C Count of packets in current loop
integer loopct
C Number of data names in hash table
integer nhash
C Number of data names in data block
integer nname
C Current number of data name in block
integer iname
C Number of dictionary names
integer ndict
C Number of records in CIF
integer nrecd
C Record number of requested line
integer irecd
C Record number of current line
integer jrecd
C Record number of tab-expanded line
integer jrect
C Last record number of current block
integer lrecd
C Character pointer of current input line
integer jchar
C Character pointer of output CIF line
C next location to store
integer pchar
C Character pointer of output CIF line
C last location stored
integer pcharl
C Number of data block names stored
integer nbloc
C Device number of input CIF
integer cifdev
C Device number of direct access file
integer dirdev
C Device number of error message file
integer errdev
C Device number of output CIF
integer outdev
C Cached copy of esdlim_
integer esdcac
C Cached esd digits
integer esddig
C Number of Dictionary Category names stored
integer ndcname
C Bookmark data
integer ibkmrk(4,MAXBOOK)
C Indices of Dictionary Category
C names for items
integer dcindex(NUMDICT)
C Hash Table for Dic. Category names
integer dchash(NUMHASH)
C Chain pointers for hash searches
C of Dictionary Category names
integer dcchain(NUMDICT)
C Hash table for Dictionary name
integer dichash(NUMHASH)
C Chain pointers for hash search of
C Dictionary names
integer dicchain(NUMDICT)
C Block hash table
integer dhash(NUMHASH)
C Category hash table
integer chash(NUMHASH)
C Chain pointers for hash search of
C block names
integer dchain(NUMBLOCK)
C Chain pointers for hash search of
C categories
integer cchain(NUMBLOCK)
C Number of Category names stored
integer ncname
C Character starting position of loop items
integer loopch(NUMIP1)
C String length of loop items
integer loopln(NUMIP1)
C Record number of loop items
integer looprd(NUMIP1)
C Indices of Category
C names for items
integer cindex(NUMBLOCK)
C Dictionary pointer from data block
integer ddict(NUMBLOCK)
C Alias links, 0 for no alias or index
C of item name in dicnam which
C is the root of its definitions
integer aroot(NUMDICT)
C Alias links, 0 for no next alias or
C index of next alias in dicnam
integer alias(NUMDICT)
C dcheck result from dicnam
integer xdchk
C minimum exponent for DP calculations
integer minexp
C length of output prefix string
integer lprefx
C Potential tabs in output line
integer itabp(MAXTAB)
C Count of potential tabs
integer numtab
C Returned number
double precision numbtb
C Returned standard deviation
double precision sdevtb
C Decimal double precision precision
double precision dpprc
C Decimal double precision minimum power of 10
double precision dpmin
C Decimal single precision precision
real decprc
C Decimal single precision minimum power of 10
real decmin
C Internal save of saveo_
logical psaveo
C
common/tbxc/ buffer,dname,dtype,tab,dicnam,dictag,dictyp,
* dictfl,dtag,nametb,testfl,vcheck,tcheck,pfilef,
* ploopf,ptextf,dcname,cname,dicxtyp,dxtyp,esdfmt,
* edpfmt,ndpfmt,prefx,obuf,bufntb
C
common/tbxi/ nrecd,drecd,trecd,irecd,lrecd,dchar,tchar,xchar,
* ndict,outdev,nloop,iloop,lloop,loopct,loopni,
* loopnp,loopnl,nname,nhash,cifdev,dirdev,errdev,
* jchar,pchar,iname,ploopn,ploopc,nbloc,jrecd,
* lastch,esdcac,esddig,ndcname,dcindex,ncname,
* cindex,aroot,alias,dchash,dcchain,dichash,
* dicchain,dhash,dchain,ddict,xdchk,loopch,loopln,
* looprd,minexp,lprefx,itabp,numtab,pcharl,loopos,
* loopox,loorec,jrect,ibkmrk,chash,cchain
C
common/tbxdp/numbtb,sdevtb,dpprc,dpmin
C
common/tbxr/ decprc,decmin
C
common/tbxl/ psaveo
C
"CUT_HERE_CUT_HERE_CUT_HERE"
#___________ THE FILE 'ciftbx.sys' ENDS TWO LINES BEFORE THIS ___________
if ( 12073 != `wc -c <'ciftbx.sys'` ) then
echo "cshar: 'ciftbx.sys' unpacked with wrong size"
endif
# end of 'ciftbx.sys'
endif
if ( -e 'clearfp.f' && ${1} != "-c" ) then
echo "cshar: Will not clobber existing file 'clearfp.f'"
sed "s/^X//" >'/dev/null' <<"CUT_HERE_CUT_HERE_CUT_HERE"
else
echo "cshar: Extracting 'clearfp.f' (383 characters)"
#___________ THE FILE 'clearfp.f' BEGINS TWO LINES AFTER THIS __________
sed "s/^X//" >'clearfp.f' <<"CUT_HERE_CUT_HERE_CUT_HERE"
subroutine clearfp
C
C This is a null version of a
C subroutine to clear IEEE floating point exceptions
C for inexact and underflow under SUN OS 4 f77
C For most other systems, no action is needed.
C
C character*1 out
C ii = ieee_flags('clear','exception','underflow',out)
C ii = ieee_flags('clear','execption','inexact',out)
return
end
"CUT_HERE_CUT_HERE_CUT_HERE"
#___________ THE FILE 'clearfp.f' ENDS TWO LINES BEFORE THIS ___________
if ( 383 != `wc -c <'clearfp.f'` ) then
echo "cshar: 'clearfp.f' unpacked with wrong size"
endif
# end of 'clearfp.f'
endif
if ( -e 'clearfp_sun.f' && ${1} != "-c" ) then
echo "cshar: Will not clobber existing file 'clearfp_sun.f'"
sed "s/^X//" >'/dev/null' <<"CUT_HERE_CUT_HERE_CUT_HERE"
else
echo "cshar: Extracting 'clearfp_sun.f' (298 characters)"
#___________ THE FILE 'clearfp_sun.f' BEGINS TWO LINES AFTER THIS __________
sed "s/^X//" >'clearfp_sun.f' <<"CUT_HERE_CUT_HERE_CUT_HERE"
subroutine clearfp
C
C subroutine to clear IEEE floating point exceptions
C for inexact and underflow under SUN OS 4 f77
C
character*1 out
ii = ieee_flags('clear','exception','underflow',out)
ii = ieee_flags('clear','execption','inexact',out)
return
end
"CUT_HERE_CUT_HERE_CUT_HERE"
#___________ THE FILE 'clearfp_sun.f' ENDS TWO LINES BEFORE THIS ___________
if ( 298 != `wc -c <'clearfp_sun.f'` ) then
echo "cshar: 'clearfp_sun.f' unpacked with wrong size"
endif
# end of 'clearfp_sun.f'
endif
if ( -e 'cryst1-2-cif.awk' && ${1} != "-c" ) then
echo "cshar: Will not clobber existing file 'cryst1-2-cif.awk'"
sed "s/^X//" >'/dev/null' <<"CUT_HERE_CUT_HERE_CUT_HERE"
else
echo "cshar: Extracting 'cryst1-2-cif.awk' (1788 characters)"
#___________ THE FILE 'cryst1-2-cif.awk' BEGINS TWO LINES AFTER THIS __________
sed "s/^X//" >'cryst1-2-cif.awk' <<"CUT_HERE_CUT_HERE_CUT_HERE"
#
# cryst1-2-cif.awk
#
# simple awk script to convert a file with PDB CRYST1 records
# to a CIF for input to iterate. Extracted from pdb2cif
# by Bourne, Bernstein and Bernstein. See
# http://ndbserver.rutgers.edu/software/pdb2cif
#
# Warning -- this script assumes the entry id is in columns 73-76
# If this is not true, then you make get a blank entry_id
#
BEGIN {
printf ("\n data_cells\n\n loop_\n")
printf ("_cell.entry_id\n")
printf ("_cell.length_a\n")
printf ("_cell.length_b\n")
printf ("_cell.length_c\n")
printf ("_cell.angle_alpha\n")
printf ("_cell.angle_beta\n")
printf ("_cell.angle_gamma\n")
printf ("_cell.volume\n")
printf ("_cell.Z_PDB\n")
printf ("_cell.space_group_name_H-M # not defined in dictionary\n\n\n")
}
#==========================================================================
# keyword CRYST1
#
#
{
if ($1 == "CRYST1") {
#
# Contains a b c alpha beta gamma SG Z
#
# calculate cell volume
{
ca = cos(substr( ($0),34, 7) * 0.0174532)
cb = cos(substr( ($0),41, 7) * 0.0174532)
cc = cos(substr( ($0),48, 7) * 0.0174532)
cz = (1.0 - (ca*ca - cb*cb - cc*cc) + (2.0*ca*cb*cc))
vol = (substr( ($0), 7, 9) *\
substr( ($0),16, 9) *\
substr( ($0),25, 9) * (sqrt(cz)))
}
# localize space group and Z
{
sg = substr( ($0), 56, 11)
Z = substr( ($0), 67, 4 )
}
if (vol-1 > .01) {
printf (" %s",substr( ($0), 73, 4))
printf (" %9.3f", substr( ($0), 7, 9))
printf (" %9.3f", substr( ($0),16, 9))
printf (" %9.3f", substr( ($0),25, 9))
printf (" %7.2f", substr( ($0),34, 7))
printf (" %7.2f", substr( ($0),41, 7))
printf (" %7.2f\n", substr( ($0),48, 7))
printf (" %10.1f", vol)
printf (" %3d ", Z)
printf (" '%11s'\n", sg)
}
}
}
END {}
"CUT_HERE_CUT_HERE_CUT_HERE"
#___________ THE FILE 'cryst1-2-cif.awk' ENDS TWO LINES BEFORE THIS ___________
if ( 1788 != `wc -c <'cryst1-2-cif.awk'` ) then
echo "cshar: 'cryst1-2-cif.awk' unpacked with wrong size"
endif
# end of 'cryst1-2-cif.awk'
endif
if ( -e 'hash_funcs.f' && ${1} != "-c" ) then
echo "cshar: Will not clobber existing file 'hash_funcs.f'"
sed "s/^X//" >'/dev/null' <<"CUT_HERE_CUT_HERE_CUT_HERE"
else
echo "cshar: Extracting 'hash_funcs.f' (7064 characters)"
#___________ THE FILE 'hash_funcs.f' BEGINS TWO LINES AFTER THIS __________
sed "s/^X//" >'hash_funcs.f' <<"CUT_HERE_CUT_HERE_CUT_HERE"
C
C hash_funcs.f -- a library of hash table management routines
C
C by
C
C Herbert J. Bernstein
C Bernstein + Sons
C P.O. Box 177, Bellport, NY 11713-0177, USA
C Phone/Fax 1-516-286-1999
C email: yaya@aip.org
C
C work on these routines done in part at Brookhaven National
C Laboratory, under contract to the U.S. Department of Energy
C
C-------------------------------------------------------------------------------
C
C Routines
C
C hash_init Initializes a hash table controlled list
C call hash_init(data_structure_args)
C
C hash_find Searches for a string in a list
C call hash_find(name,data_structure_args,ifind)
C
C hash_store Inserts as new string in a list
C call hash_store(name,data_structure_args,ifind)
C
C hash_value Integer function returns index into hash_list
C ih = hash_value(name,hash_length)
C
C The necessary data_structure_args for these routines are
C name_list -- an array of character strings
C character*(*) name_list(list_length)
C chain_list -- chain pointers for searches
C integer chain_list(list_length)
C list_length -- the size of the list arrays
C integer list_length
C num_list -- number of entries in the list
C integer num_list
C hash_table -- the initial hashed pointers
C integer hash_table
C hash_length -- the size of the hash table
C integer hash_length
C
C
C The two remaining arguments are
C name -- string to search for
C character*(*) name
C ifind -- return value, 0 for not found (hash_find)
C or list full (hash_store), otherwise
C the index in name_list of the entry
C
C The relationship among the arrays used is:
C
C hash_table is an array (preferably of a modest prime
C dimension) which starts containing all zeros, which are
C replaced by pointers to entries in name_list, based
C values returned by hash_value ranging from 1 to hash_length.
C Each name is placed in name_list. A initial zero is placed
C in the matching entry in chain_list, when the first entry
C is made. When a new entry with the same hash_value must be
C placed a pointer is inserted into chain_list to hook the
C values together.
C
subroutine hash_init(name_list,chain_list,list_length,num_list,
* hash_table,hash_length)
C
C initialization routine for a hash table controlled list
C name_list -- a list of character strings
C chain_list -- chain pointers for searches
C list_length -- the size of the list arrays
C num_list -- number of entries in the list
C hash_table -- the initial hashed pointers
C hash_length -- the size of the hash table
C
character*(*) name_list(list_length)
integer hash_length,list_length,num_list,i
integer chain_list(list_length)
integer hash_table(hash_length)
num_list=0
do i = 1,hash_length
hash_table(i)=0
enddo
return
end
subroutine
* hash_find(name,name_list,chain_list,list_length,num_list,
* hash_table,hash_length,ifind)
C
C search routine for a hash table controlled list
C name -- string to find
C name_list -- a list of character strings
C chain_list -- chain pointers for searches
C list_length -- the size of the list arrays
C num_list -- number of entries in the list
C hash_table -- the initial hashed pointers
C hash_length -- the size of the hash table
C ifind -- returned index or 0
C
character*(*) name
integer hash_length
character*(*) name_list(list_length)
integer chain_list(list_length)
integer hash_table(hash_length)
integer hash_value
integer ifind,list_length,num_list,ih,ip
ifind=0
ih=hash_value(name,hash_length)
ip=hash_table(ih)
100 if (ip.eq.0) return
if (name_list(ip).eq.name) then
ifind=ip
return
else
ip=chain_list(ip)
go to 100
endif
end
subroutine
* hash_store(name,name_list,chain_list,list_length,num_list,
* hash_table,hash_length,ifind)
C
C store routine for a hash table controlled list
C name -- string to find
C name_list -- a list of character strings
C chain_list -- chain pointers for searches
C list_length -- the size of the list arrays
C num_list -- number of entries in list
C hash_table -- the initial hashed pointers
C hash_length -- the size of the hash table
C ifind -- index of entry or 0 (table full)
C
character*(*) name
character*(*) name_list(list_length)
integer hash_length
integer chain_list(list_length)
integer hash_table(hash_length)
integer hash_value
integer ifind,list_length,num_list,ih,ip,iq
ifind=0
ih = hash_value(name,hash_length)
ip=hash_table(ih)
iq=0
100 if (ip.eq.0) go to 200
if (name_list(ip).eq.name) then
ifind=ip
return
else
iq=ip
ip=chain_list(ip)
go to 100
endif
200 if (num_list.lt.list_length) then
num_list=num_list+1
name_list(num_list)=name
chain_list(num_list)=0
if (iq.eq.0) then
hash_table(ih)=num_list
else
chain_list(iq)=num_list
endif
ifind=num_list
return
else
ifind = 0
return
endif
end
integer function hash_value(name,hash_length)
C
C function to return a hash value of string name to fit
C a hash table of length hash_length
character*(*) name
integer hash_length,id,ii,i,ic,lenn
lenn = len(name)
hash_value=1
id = 0
do ii = 1,lenn
i = 1+lenn-ii
ic = ichar(name(i:i))
if (ic.ge.65) then
hash_value=mod(hash_value*(ic-64),hash_length)+1
id = id+1
if (id.gt.3) return
endif
enddo
return
end
"CUT_HERE_CUT_HERE_CUT_HERE"
#___________ THE FILE 'hash_funcs.f' ENDS TWO LINES BEFORE THIS ___________
if ( 7064 != `wc -c <'hash_funcs.f'` ) then
echo "cshar: 'hash_funcs.f' unpacked with wrong size"
endif
# end of 'hash_funcs.f'
endif
if ( -e 'iterate.csh.m4' && ${1} != "-c" ) then
echo "cshar: Will not clobber existing file 'iterate.csh.m4'"
sed "s/^X//" >'/dev/null' <<"CUT_HERE_CUT_HERE_CUT_HERE"
else
echo "cshar: Extracting 'iterate.csh.m4' (1292 characters)"
#___________ THE FILE 'iterate.csh.m4' BEGINS TWO LINES AFTER THIS __________
sed "s/^X//" >'iterate.csh.m4' <<"CUT_HERE_CUT_HERE_CUT_HERE"
#!/bin/csh
# iterate.csh
#
# Herbert J. Bernstein, Bernstein + Sons
# Lawrence C. Andrews, Thuridion, Inc.
#
# 29 September 1996
#
# This is a service script for the iterate.html web page
# It must be placed in an appropriate cgi-bin directory on
# the server pointed to by iterate.html
#
#
# To operate correctly, the programs tr and sed must be in the
# default path and the /bin/echo version of echo must follow
# system V conventions sufficiently to produce an empty line
# call, below
#
/bin/echo "Content-type: text/html"
/bin/echo
echo ""
echo "G6 Lattice Identification"
echo ""
echo ""
echo ""
tr '\&' '\n' |sed "s/^./set &/" > /tmp/outstr$$
#cat /tmp/outstr$$
source /tmp/outstr$$
rm /tmp/outstr$$
echo "# G6 Lattice Identification
"
echo $Centering > /tmp/instr$$
echo "# Centering: " $Centering
echo $A $B $C $Alpha $Beta $Gamma >>/tmp/instr$$
echo "
# Cell: " $A $B $C $Alpha $Beta $Gamma
echo $sigA $sigB $sigC $sigAlpha $sigBeta $sigGamma >>/tmp/instr$$
echo "
# Sigmas: " $sigA $sigB $sigC $sigAlpha $sigBeta $sigGamma
echo "
# Results of ITERATE Run
"
setenv ITERATE_QUERY NO
setenv OUTPUT_STYLE $OutputStyle
echo ""
BINPATH() < /tmp/instr$$
rm /tmp/instr$$
#cat /tmp/instr$$
echo "
"
echo ""
"CUT_HERE_CUT_HERE_CUT_HERE"
#___________ THE FILE 'iterate.csh.m4' ENDS TWO LINES BEFORE THIS ___________
if ( 1292 != `wc -c <'iterate.csh.m4'` ) then
echo "cshar: 'iterate.csh.m4' unpacked with wrong size"
endif
chmod +x 'iterate.csh.m4'
# end of 'iterate.csh.m4'
endif
if ( -e 'iterate.f' && ${1} != "-c" ) then
echo "cshar: Will not clobber existing file 'iterate.f'"
sed "s/^X//" >'/dev/null' <<"CUT_HERE_CUT_HERE_CUT_HERE"
else
echo "cshar: Extracting 'iterate.f' (101777 characters)"
#___________ THE FILE 'iterate.f' BEGINS TWO LINES AFTER THIS __________
sed "s/^X//" >'iterate.f' <<"CUT_HERE_CUT_HERE_CUT_HERE"
C ITERATE -- Program for G6 Bravais Lattice Determination
C
C by
C
C Lawrence C. Andrews, Thuridion, Inc.,
C andrews@thuridion.com
C
C and
C
C Herbert J. Bernstein, Bernstein+Sons,
C yaya@bernstein-plus-sons.com
C
C This program finds the cells which are "close" to the cell given,
C in order to help find the Bravais lattice of highest symmetry
C consistent with the cell.
C
C A central problem in the solution of every crystal structure
C is to determine the correct Bravais lattice of the crystal.
C The Bravais lattices as they are usually listed are:
C
C aP triclinic (anorthic) primitive
C mP monoclinic primitive
C mS monoclinic side-centered (usually C-centered)
C oP orthorhombic primitive
C oS orthorhombic side-centered
C oF orthorhombic face-centered
C oI orthorhombic body-centered
C hP hexagonal primitive
C hR hexagonal rhombohedrally-centered
C tP tetragonal primitive
C tI tetragonal body-centered
C cP cubic primitive
C cF cubic face-centered
C cI cubic body-centered
C
C
C Failure to find the highest correct symmetry has several consequences,
C the worst of which is that the structure may not be solved. The
C least of the consequences is that Richard Marsh may publish a paper
C that points out the error, corrects it, and finds a better solution
C to the structure. Many methods have been described for finding the
C correct Bravais lattice. A summary of the published methods was
C published in the paper that described the G6 formalism (which is used
C in this program).
C
C "Lattices and Reduced Cells as Points in 6-Space and Selection of
C Bravais Lattice Type by Projections", Lawrence C. Andrews and
C Herbert J. Bernstein, Acta Crystallographica, A44, 1009-1018 (1988).
C
C This program accepts cell parameters and esd's and produces a list
C of cells "close" to the cell given
C
C**********************************************************************C
SUBROUTINE BADCAL (A,B)
CHARACTER *6 A,B
C----------------------------------------------------------------------C
include 'ITERATE.cmn'
WRITE (*,*) ' '//hm//' BAD SUBROUTINE CALL TO ',B
WRITE (*,*) ' '//hm//' CALLING NAME =',A
STOP
END
CC***********************************************************************
C SUBROUTINE CPYVN (N,V1,V2)
C REAL V1(N),V2(N)
C
C-----------------------------------------------------------------------
C DO 1000 I=1,N
C 1000 V2(I) = V1(I)
C END
C***********************************************************************
subroutine bldprj (maxprj,nproj,itdesg,chrlat,pjn,prj,test)
C This function builds the projectors of Paciorek and Bonin, J. Appl.
C Cryst., 25, (5) pp 632-637. Internal checks of the correctness of
C the projectors are made before it exits.
C
C The output values are just transferred from the stored parameters
C in data statements. Note that if ngtype is negative, then the
C projector is not output. This is because those are ones that
C are easy to find later by searching. Either they are cases where
C two projectors are exactly the same subspace (and therefore the
C same projector) or else they are just simple exchange of axes
C (which are dealt with in mkrefl in program iterate). If the routine
C is extracted for other uses, then it may be correct to set the
C values all positive.
C maxprj sets the maximum number of projectors to build
C nproj is the actual number that bldprj constructs
C itdesg is the numeric designation of the Niggli type in the
C International Tables for Crystallography
C chrlat is the returned 2 character designators for the lattice
C type (also call Pearson symbols)
C pjn is the normalizer of the projector -- the integer values
C stored in the returned matrix need to be divided by the
C normalizer to make the actual projector
C prj the integer part of the projector -- divide by the normalizer
C (pjn) to get the actual projector
include 'ITERATE.cmn'
real pjn(maxprj)
real prj(36,MAXPRJ)
integer itdesg(maxprj)
character *2 chrlat(maxprj)
integer ngtype(42)
real projct(36,42)
real zprj(6,6,42)
equivalence (projct,zprj)
real pjnorm(42)
character *2 lattyp(42)
character *6 test
data ngtype(1) /3/
data lattyp(1) /'cP'/
data pjnorm(1) /3./
data( projct(i,1),i=1,36) /
1 1,1,1,0,0,0,
2 1,1,1,0,0,0,
3 1,1,1,0,0,0,
4 0,0,0,0,0,0,
5 0,0,0,0,0,0,
6 0,0,0,0,0,0 /
data ngtype(2) /5/
data lattyp(2) /'cI'/
data pjnorm(2) /39./
data (projct(i,2),i=1,36) /
1 9, 9, 9,-6,-6,-6,
2 9, 9, 9,-6,-6,-6,
3 9, 9, 9,-6,-6,-6,
4 -6,-6,-6, 4, 4, 4,
5 -6,-6,-6, 4, 4, 4,
6 -6,-6,-6, 4, 4, 4 /
data ngtype(3) /1/
data lattyp(3) /'cF'/
data pjnorm(3) /6./
data (projct(i,3),i=1,36) /
1 1,1,1,1,1,1,
2 1,1,1,1,1,1,
3 1,1,1,1,1,1,
4 1,1,1,1,1,1,
5 1,1,1,1,1,1,
6 1,1,1,1,1,1 /
data ngtype(4) /11/
data lattyp(4) /'tP'/
data pjnorm(4) /2./
data (projct(i,4),i=1,36) /
1 1,1,0,0,0,0,
2 1,1,0,0,0,0,
3 0,0,2,0,0,0,
4 0,0,0,0,0,0,
5 0,0,0,0,0,0,
6 0,0,0,0,0,0 /
data ngtype(5) /-21/
data lattyp(5) /'tP'/
data pjnorm(5) /2./
data (projct(i,5),i=1,36) /
1 2,0,0,0,0,0,
2 0,1,1,0,0,0,
3 0,1,1,0,0,0,
4 0,0,0,0,0,0,
5 0,0,0,0,0,0,
6 0,0,0,0,0,0 /
data ngtype(6) /15/
data lattyp(6) /'tI'/
data pjnorm(6) /4./
data (projct(i,6),i=1,36) /
1 1, 1, 0,-1,-1, 0,
2 1, 1, 0,-1,-1, 0,
3 0, 0, 4, 0, 0, 0,
4 -1,-1, 0, 1, 1, 0,
5 -1,-1, 0, 1, 1, 0,
6 0, 0, 0, 0, 0, 0 /
data ngtype(7) /6/
data lattyp(7) /'tI'/
data pjnorm(7) /26./
data (projct(i,7),i=1,36) /
1 6, 6, 6,-4,-4,-4,
2 6, 6, 6,-4,-4,-4,
3 6, 6, 6,-4,-4,-4,
4 -4,-4,-4, 7, 7,-6,
5 -4,-4,-4, 7, 7,-6,
6 -4,-4,-4,-6,-6,20 /
data ngtype(8) /-7/
data lattyp(8) /'tI'/
data pjnorm(8) /26./
data (projct(i,8),i=1,36) /
1 6, 6, 6,-4,-4,-4,
2 6, 6, 6,-4,-4,-4,
3 6, 6, 6,-4,-4,-4,
4 -4,-4,-4,20,-6,-6,
5 -4,-4,-4,-6, 7, 7,
6 -4,-4,-4,-6, 7, 7 /
data ngtype(9) /18/
data lattyp(9) /'tI'/
data pjnorm(9) /26./
data (projct(i,9),i=1,36) /
1 8, 0, 0, 4, 8, 8,
2 0,13,13, 0, 0, 0,
3 0,13,13, 0, 0, 0,
4 4, 0, 0, 2, 4, 4,
5 8, 0, 0, 4, 8, 8,
6 8, 0, 0, 4, 8, 8 /
data ngtype(10) /12/
data lattyp(10) /'hP'/
data pjnorm(10) /3./
data (projct(i,10),i=1,36) /
1 1, 1, 0, 0, 0,-1,
2 1, 1, 0, 0, 0,-1,
3 0, 0, 3, 0, 0, 0,
4 0, 0, 0, 0, 0, 0,
5 0, 0, 0, 0, 0, 0,
6 -1,-1, 0, 0, 0, 1 /
data ngtype(11) /-22/
data lattyp(11) /'hP'/
data pjnorm(11) /3./
data (projct(i,11),i=1,36) /
1 3, 0, 0, 0, 0, 0,
2 0, 1, 1,-1, 0, 0,
3 0, 1, 1,-1, 0, 0,
4 0,-1,-1, 1, 0, 0,
5 0, 0, 0, 0, 0, 0,
6 0, 0, 0, 0, 0, 0 /
data ngtype(12) /9/
data lattyp(12) /'hR'/
data pjnorm(12) /5./
data (projct(i,12),i=1,36) /
1 1,1,0,1,1,1,
2 1,1,0,1,1,1,
3 0,0,5,0,0,0,
4 1,1,0,1,1,1,
5 1,1,0,1,1,1,
6 1,1,0,1,1,1 /
data ngtype(13) /2/
data lattyp(13) /'hR'/
data pjnorm(13) /3./
data (projct(i,13),i=1,36) /
1 1,1,1,0,0,0,
2 1,1,1,0,0,0,
3 1,1,1,0,0,0,
4 0,0,0,1,1,1,
5 0,0,0,1,1,1,
6 0,0,0,1,1,1 /
data ngtype(14) /-4/
data lattyp(14) /'hR'/
data pjnorm(14) /3./
data (projct(i,14),i=1,36) /
1 1,1,1,0,0,0,
2 1,1,1,0,0,0,
3 1,1,1,0,0,0,
4 0,0,0,1,1,1,
5 0,0,0,1,1,1,
6 0,0,0,1,1,1 /
data ngtype(15) /24/
data lattyp(15) /'hR'/
data pjnorm(15) /53./
data (projct(i,15),i=1,36) /
1 27, 3, 3, 6,-18,-18,
2 3, 18, 18,-17, -2, -2,
3 3, 18, 18,-17, -2, -2,
4 6,-17,-17, 19, -4, -4,
5 -18, -2, -2, -4, 12, 12,
6 -18, -2, -2, -4, 12, 12 /
data ngtype(16) /32/
data lattyp(16) /'oP'/
data pjnorm(16) /1./
data (projct(i,16),i=1,36) /
1 1,0,0,0,0,0,
2 0,1,0,0,0,0,
3 0,0,1,0,0,0,
4 0,0,0,0,0,0,
5 0,0,0,0,0,0,
6 0,0,0,0,0,0 /
data ngtype(17) /36/
data lattyp(17) /'oS'/
data pjnorm(17) /2./
data (projct(i,17),i=1,36) /
1 1, 0, 0, 0,-1, 0,
2 0, 2, 0, 0, 0, 0,
3 0, 0, 2, 0, 0, 0,
4 0, 0, 0, 0, 0, 0,
5 -1, 0, 0, 0, 1, 0,
6 0, 0, 0, 0, 0, 0 /
data ngtype(18) /-38/
data lattyp(18) /'oS'/
data pjnorm(18) /2./
data (projct(i,18),i=1,36) /
1 1, 0, 0, 0, 0,-1,
2 0, 2, 0, 0, 0, 0,
3 0, 0, 2, 0, 0, 0,
4 0, 0, 0, 0, 0, 0,
5 0, 0, 0, 0, 0, 0,
6 -1, 0, 0, 0, 0, 1 /
data ngtype(19) /13/
data lattyp(19) /'oS'/
data pjnorm(19) /2./
data (projct(i,19),i=1,36) /
1 1,1,0,0,0,0,
2 1,1,0,0,0,0,
3 0,0,2,0,0,0,
4 0,0,0,0,0,0,
5 0,0,0,0,0,0,
6 0,0,0,0,0,2 /
data ngtype(20) /-23/
data lattyp(20) /'oS'/
data pjnorm(20) /2./
data (projct(i,20),i=1,36) /
1 2,0,0,0,0,0,
2 0,1,1,0,0,0,
3 0,1,1,0,0,0,
4 0,0,0,2,0,0,
5 0,0,0,0,0,0,
6 0,0,0,0,0,0 /
data ngtype(21) /-40/
data lattyp(21) /'oS'/
data pjnorm(21) /2./
data (projct(i,21),i=1,36) /
1 2, 0, 0, 0, 0, 0,
2 0, 1, 0,-1, 0, 0,
3 0, 0, 2, 0, 0, 0,
4 0,-1, 0, 1, 0, 0,
5 0, 0, 0, 0, 0, 0,
6 0, 0, 0, 0, 0, 0 /
data ngtype(22) /16/
data lattyp(22) /'oF'/
data pjnorm(22) /10./
data (projct(i,22),i=1,36) /
1 3, 3, 0,-2,-2,-2,
2 3, 3, 0,-2,-2,-2,
3 0, 0,10, 0, 0, 0,
4 -2,-2, 0, 3, 3,-2,
5 -2,-2, 0, 3, 3,-2,
6 -2,-2, 0,-2,-2, 8 /
data ngtype(23) /26/
data lattyp(23) /'oF'/
data pjnorm(23) /13./
data (projct(i,23),i=1,36) /
1 4, 0, 0, 2, 4, 4,
2 0,13, 0, 0, 0, 0,
3 0, 0,13, 0, 0, 0,
4 2, 0, 0, 1, 2, 2,
5 4, 0, 0, 2, 4, 4,
6 4, 0, 0, 2, 4, 4 /
data ngtype(24) /8/
data lattyp(24) /'oI'/
data pjnorm(24) /13./
data (projct(i,24),i=1,36) /
1 3, 3, 3,-2,-2,-2,
2 3, 3, 3,-2,-2,-2,
3 3, 3, 3,-2,-2,-2,
4 -2,-2,-2,10,-3,-3,
5 -2,-2,-2,-3,10,-3,
6 -2,-2,-2,-3,-3,10 /
data ngtype(25) /19/
data lattyp(25) /'oI'/
data pjnorm(25) /6./
data (projct(i,25),i=1,36) /
1 2,0,0,0,2,2,
2 0,3,3,0,0,0,
3 0,3,3,0,0,0,
4 0,0,0,6,0,0,
5 2,0,0,0,2,2,
6 2,0,0,0,2,2 /
data ngtype(26) /42/
data lattyp(26) /'oI'/
data pjnorm(26) /2./
data (projct(i,26),i=1,36) /
1 1, 0, 0, 0,-1, 0,
2 0, 1, 0,-1, 0, 0,
3 0, 0, 2, 0, 0, 0,
4 0,-1, 0, 1, 0, 0,
5 -1, 0, 0, 0, 1, 0,
6 0, 0, 0, 0, 0, 0 /
data ngtype(27) /33/
data lattyp(27) /'mP'/
data pjnorm(27) /1./
data (projct(i,27),i=1,36) /
1 1,0,0,0,0,0,
2 0,1,0,0,0,0,
3 0,0,1,0,0,0,
4 0,0,0,0,0,0,
5 0,0,0,0,1,0,
6 0,0,0,0,0,0 /
data ngtype(28) /-35/
data lattyp(28) /'mP'/
data pjnorm(28) /1./
data (projct(i,28),i=1,36) /
1 1,0,0,0,0,0,
2 0,1,0,0,0,0,
3 0,0,1,0,0,0,
4 0,0,0,1,0,0,
5 0,0,0,0,0,0,
6 0,0,0,0,0,0 /
data ngtype(29) /-34/
data lattyp(29) /'mP'/
data pjnorm(29) /1./
data (projct(i,29),i=1,36) /
1 1,0,0,0,0,0,
2 0,1,0,0,0,0,
3 0,0,1,0,0,0,
4 0,0,0,0,0,0,
5 0,0,0,0,0,0,
6 0,0,0,0,0,1 /
data ngtype(30) /39/
data lattyp(30) /'mS'/
data pjnorm(30) /2./
data (projct(i,30),i=1,36) /
1 1, 0, 0, 0, 0,-1,
2 0, 2, 0, 0, 0, 0,
3 0, 0, 2, 0, 0, 0,
4 0, 0, 0, 2, 0, 0,
5 0, 0, 0, 0, 0, 0,
6 -1, 0, 0, 0, 0, 1 /
data ngtype(31) /-41/
data lattyp(31) /'mS'/
data pjnorm(31) /2./
data (projct(i,31),i=1,36) /
1 2, 0, 0, 0, 0, 0,
2 0, 1, 0,-1, 0, 0,
3 0, 0, 2, 0, 0, 0,
4 0,-1, 0, 1, 0, 0,
5 0, 0, 0, 0, 2, 0,
6 0, 0, 0, 0, 0, 0 /
data ngtype(32) /-37/
data lattyp(32) /'mS'/
data pjnorm(32) /2./
data (projct(i,32),i=1,36) /
1 1, 0, 0, 0,-1, 0,
2 0, 2, 0, 0, 0, 0,
3 0, 0, 2, 0, 0, 0,
4 0, 0, 0, 2, 0, 0,
5 -1, 0, 0, 0, 1, 0,
6 0, 0, 0, 0, 0, 0 /
data ngtype(33) /10/
data lattyp(33) /'mS'/
data pjnorm(33) /2./
data (projct(i,33),i=1,36) /
1 1,1,0,0,0,0,
2 1,1,0,0,0,0,
3 0,0,2,0,0,0,
4 0,0,0,1,1,0,
5 0,0,0,1,1,0,
6 0,0,0,0,0,2 /
data ngtype(34) /-14/
data lattyp(34) /'mS'/
data pjnorm(34) /2./
data (projct(i,34),i=1,36) /
1 1,1,0,0,0,0,
2 1,1,0,0,0,0,
3 0,0,2,0,0,0,
4 0,0,0,1,1,0,
5 0,0,0,1,1,0,
6 0,0,0,0,0,2 /
data ngtype(35) /-20/
data lattyp(35) /'mS'/
data pjnorm(35) /2./
data (projct(i,35),i=1,36) /
1 2,0,0,0,0,0,
2 0,1,1,0,0,0,
3 0,1,1,0,0,0,
4 0,0,0,2,0,0,
5 0,0,0,0,1,1,
6 0,0,0,0,1,1 /
data ngtype(36) /-25/
data lattyp(36) /'mS'/
data pjnorm(36) /2./
data (projct(i,36),i=1,36) /
1 2,0,0,0,0,0,
2 0,1,1,0,0,0,
3 0,1,1,0,0,0,
4 0,0,0,2,0,0,
5 0,0,0,0,1,1,
6 0,0,0,0,1,1 /
data ngtype(37) /28/
data lattyp(37) /'mS'/
data pjnorm(37) /10./
data (projct(i,37),i=1,36) /
1 5, 0, 0, 0, 5, 0,
2 0,10, 0, 0, 0, 0,
3 0, 0,10, 0, 0, 0,
4 0, 0, 0, 2, 0, 4,
5 5, 0, 0, 0, 5, 0,
6 0, 0, 0, 4, 0, 8 /
data ngtype(38) /-30/
data lattyp(38) /'mS'/
data pjnorm(38) /10./
data (projct(i,38),i=1,36) /
1 10, 0, 0, 0, 0, 0,
2 0, 5, 0, 5, 0, 0,
3 0, 0,10, 0, 0, 0,
4 0, 5, 0, 5, 0, 0,
5 0, 0, 0, 0, 2, 4,
6 0, 0, 0, 0, 4, 8 /
data ngtype(39) /-29/
data lattyp(39) /'mS'/
data pjnorm(39) /10./
data (projct(i,39),i=1,36) /
1 5, 0, 0, 0, 0, 5,
2 0,10, 0, 0, 0, 0,
3 0, 0,10, 0, 0, 0,
4 0, 0, 0, 2, 4, 0,
5 0, 0, 0, 4, 8, 0,
6 5, 0, 0, 0, 0, 5 /
data ngtype(40) /43/
data lattyp(40) /'mI'/
data pjnorm(40) /20./
data (projct(i,40),i=1,36) /
1 11, 1, 0, 1,-9,-4,
2 1,11, 0,-9, 1,-4,
3 0, 0,20, 0, 0, 0,
4 1,-9, 0,11, 1,-4,
5 -9, 1, 0, 1,11,-4,
6 -4,-4, 0,-4,-4,16 /
data ngtype(41) /17/
data lattyp(41) /'mI'/
data pjnorm(41) /10./
data (projct(i,41),i=1,36) /
1 3, 3, 0,-2,-2,-2,
2 3, 3, 0,-2,-2,-2,
3 0, 0,10, 0, 0, 0,
4 -2,-2, 0, 8,-2,-2,
5 -2,-2, 0,-2, 8,-2,
6 -2,-2, 0,-2,-2, 8 /
data ngtype(42) /27/
data lattyp(42) /'mI'/
data pjnorm(42) /3./
data (projct(i,42),i=1,36) /
1 1,0,0,0,1,1,
2 0,3,0,0,0,0,
3 0,0,3,0,0,0,
4 0,0,0,3,0,0,
5 1,0,0,0,1,1,
6 1,0,0,0,1,1 /
C-----------------------------------------------------------------------
nprob = 0
if (test .ne. 'BLDPRJ') then
write (*,*)
* ' '//hm//' test string was not BLDPRJ in that routine'
stop
endif
C check the projectors for internal correctness
do 4000 iproj=1,42
C there are only 42 niggli lattice types (ignoring triclinic)
if (ngtype(iproj).gt.43) then
write (*,*)
* ' '//hm//' bad ngtype ',iproj,ngtype(iproj)
nprob = nprob + 1
endif
C check that the lattice type is present
if (lattyp(iproj).eq.' ') then
write (*,*) ' '//hm//' blank lattyp ',iproj
nprob = nprob + 1
endif
C check that the normalizers are in the range in Paciorek and Bonin
if (pjnorm(iproj) .le. 0 .or. pjnorm(iproj).gt. 60) then
write (*,*) ' '//hm//' bad pjnorm ',iproj,' ',pjnorm(iproj)
nprob = nprob + 1
endif
C check that the projectors are symmetrical matrices
do 1000 i=1,5
do 1000 j=i+1,6
if (zprj(j,i,iproj) .ne. zprj(i,j,iproj)) then
write (*,*) ' '//hm//' bad projector ',iproj,' ',i,' ',j
write (*,*) zprj(j,i,iproj),' ',zprj(i,j,iproj)
nprob = nprob + 1
endif
1000 continue
do 2000 i=1,36
C check that the actual projector has no value greater than 1.0
if (abs(projct(i,iproj)/(pjnorm(iproj))).gt. 1.00001) then
write (*,*) ' '//hm//' bad projector, value > 1.0 '
write (*,*) ' '//hm,iproj,' pjnorm ',pjnorm(iproj),
2 ' i ',i,' ',projct(i,iproj)
write (*,*)
nprob = nprob + 1
endif
2000 continue
C check that the projector is positive definite
do 3500 i=1,3
sum = 0.0
do 3000 j=1,6
sum = sum + abs(zprj(i,j,iproj))
3000 continue
if (sum .le. 0) then
write (*,*) ' '//hm//' zero xyz row, ',iproj
nprob = nprob + 1
endif
3500 continue
4000 continue
C actually output the projectors
NPROJ = 0
do 5100 i=1,MIN(MAXPRJ,42)
IF (NGTYPE(I) .GT. 0) THEN
NPROJ = NPROJ + 1
itdesg(NPROJ) = ngtype(i)
chrlat(NPROJ) = lattyp(i)
pjn(NPROJ) = pjnorm(i)
do 5000 ii=1,36
prj(ii,nproj) = projct(ii,i)
5000 continue
ENDIF
5100 continue
if (nprob .gt. 0) then
write (*,*) ' '//hm//' ',
* nprob,' problems were found with projectors'
stop
endif
itemp = projct(1,1)
projct(1,1) = itemp
end
C**********************************************************************C
SUBROUTINE BLDTRE (MXTREE,NVEC,X,IDIN,TREE,TEST)
C BLDTRE is called once for each point to be loaded into its internal
C data structure (TREE). It builds the tree structure of Kalantari
C and McDonald (IEEE Transactions on Software Engineering, v. SE-9,
C pp. 631-634,1983) for the extremely fast retrieval of coordinates.
C MXTREE is the largest index that is allowed in the array TREE.
C TREE must be the order of 9-10 times the number of points to be
C included.
C X is an input point's coordinates.
C IDIN is an arbitrary integer input, which will often be an array
C index to be retrieved later. BLDTRE does not examine IDIN.
C TEST must be the string 'BLDTRE' -- it is used to make sure
C that the number of formal parameters is correct.
C TREE is used by NEARST and by INSPHR to find the nearest neighbor
C to a probe point. To initialize (or reinitialize) a TREE, set
C TREE(1) equal to 0.0
include 'ITERATE.cmn'
CHARACTER*6 TEST
LOGICAL DEBUG
REAL TREE(MXTREE)
INTEGER RMAX
DATA LINK,RMAX,ID,ICHILD /1,2,3,4/
DATA DEBUG /.false./
C----------------------------------------------------------------------C
IF (TEST .NE. 'BLDTRE' .AND. TEST .NE. 'bldtre')
2 CALL BADCAL (TEST,'BLDTRE')
C THE NODE SIZE IS 4 PLUS THE SIZE OF THE VECTOR
NODSIZ = 4+NVEC
IPOINT = 2
IF (TREE(1) .GT. 0) THEN
IFREE = TREE(1)
ELSE
IFREE = 2
ENDIF
tree(ifree) = 0
1000 CONTINUE
IF (DEBUG) WRITE (*,*)
2 ' '//hm//' AFTER 1000 IN BLDTRE, IPOINT,TREE(IPOINT) ',
3 IPOINT,TREE(IPOINT)
IF (DEBUG) WRITE (*,*) ' '//hm//' IFREE,TREE(1),TREE(2) ',
2 IFREE,INT(TREE(1)),INT(TREE(2))
IF (TREE(IPOINT) .EQ. 0) THEN
IF (DEBUG) WRITE (*,*)
* ' '//hm//' A NEW NODE IS BEING ALLOCATED'
IPOINT = IFREE
TREE(IPOINT) = -1
TREE(IPOINT+LINK) = -1
TREE(IPOINT+ID) = IDIN
CALL CPYVN (NVEC,X,TREE(IPOINT+ICHILD))
TREE(1) = IFREE + NODSIZ
RETURN
ELSEIF (TREE(IPOINT) .EQ. -1) THEN
IF (DEBUG) WRITE (*,*)
* ' '//hm//' RIGHT CHILD OF NODE IS BEING FILLED'
TREE(IPOINT) = IFREE
TREE(IFREE+LINK) = -1
TREE (IFREE+RMAX) = -1.0
TREE(IFREE+ID) = IDIN
CALL CPYVN (NVEC,X,TREE(IFREE+ICHILD))
IFREE = IFREE + NODSIZ
TREE(1) = IFREE
RETURN
ELSE
IRIGHT = TREE(IPOINT)
DL = TREELN (NVEC,X,TREE(IPOINT+ICHILD))
DR = TREELN (NVEC,X,TREE(IRIGHT+ICHILD))
IF (DEBUG) WRITE (*,*) ' '//hm//' DL,DR ',DL,DR
IF (DR .GT. DL) THEN
IF (DEBUG) WRITE (*,*)
* ' '//hm//' ',DR,DL,IPOINT,LINK,RMAX
IF (TREE(IPOINT+LINK) .LE. 0) THEN
TREE(IPOINT+RMAX) = DL
TREE(IPOINT+LINK) = IFREE
IPOINT = IFREE
ELSE
TREE(IPOINT+RMAX) =
2 MAX(DL,TREE(IPOINT+RMAX))
IPOINT = TREE(IPOINT+LINK)
ENDIF
ELSE
IF (TREE(IRIGHT+LINK) .LE. 0) THEN
TREE(IRIGHT+RMAX) = DR
TREE(IRIGHT+LINK) = IFREE
IPOINT = IFREE
ELSE
TREE(IRIGHT+RMAX) =
2 MAX(DR,TREE(IRIGHT+RMAX))
IPOINT = TREE(IRIGHT+LINK)
ENDIF
ENDIF
GO TO 1000
ENDIF
END
C***********************************************************************
SUBROUTINE CHKVEC(V)
C Check that a g6 vector represents a valid cell. Currently, it only
C checks that the cell edges are real
include 'ITERATE.cmn'
REAL V(6)
C-----------------------------------------------------------------------
NBAD = 0
DO 1000 I=1,3
IF(V(I).LE. 0.0) THEN
WRITE (*,*) ' '//hm//' BAD VECTOR, I=',I,' ',V(I)
NBAD = NBAD + 1
ENDIF
1000 CONTINUE
IF (NBAD .GT. 0) STOP
END
C**********************************************************************C
SUBROUTINE CPYVN (NVEC,X,Y)
C----COPY A VECTOR X INTO A VECTOR Y
DIMENSION X(NVEC), Y(NVEC)
C----------------------------------------------------------------------C
DO 1000 I=1,NVEC
Y(I) = X(I)
1000 CONTINUE
END
C***********************************************************************
SUBROUTINE CTOG6 (CV,CVE,G,GE,SIZE,ERRSIZ,RATIO,TEST)
C Convert from a unit cell (edge lengths and angles) to a g6 vector,
C also computing the errors in the g6 vector. The length of the vector
C and the error in the length of the vector is what iterate is really
C going to use.
include 'ITERATE.cmn'
CHARACTER *6 TEST
REAL COSI(4:6)
REAL C(6),CV(6),CE(6),CVE(6),G(6),GE(6)
C-----------------------------------------------------------------------
IF (TEST .NE. 'CTOG6 ') THEN
WRITE (*,*) ' '//hm//' TEST IS WRONG IN CTOG6'
STOP
ENDIF
RAD = ATAN2(0.0,-1.0) / 180.0
DO 1000 I=1,6
C(I) = CV(I)
CE(I) = CVE(I)
1000 CONTINUE
DO 1100 I=4,6
C(I) = C(I) * RAD
CE(I) = CE(I) * RAD
1100 CONTINUE
DO 1500 I=1,3
J = I + 3
G(I) = C(I)*C(I)
COSI(J) = COS(C(J))
IF (C(I).NE.0.0) THEN
G(J)=2.0*C(1)*C(2)*C(3)*COSI(J)/C(I)
ELSE
G(J) = 0.0
ENDIF
GE(I) = 2.*ABS(C(I)*CE(I))
C WRITE (*,*) ' '//hm//' GE(I) ',I,' ',GE(I)
1500 CONTINUE
DO 2000 I=1,3
J = I + 3
IF (J.EQ.4) THEN
I1 = 2
I2 = 3
ELSEIF (J.EQ.5) THEN
I1 = 1
I2 = 3
ELSE
I1 = 1
I2 = 2
ENDIF
GE(J) = 2.0*SQRT(G(I1)*(COSI(J)*CE(I2))**2 +
2 G(I2)*(COSI(J)*CE(I1))**2 +
3 G(I1)*G(I2)*(SIN(C(J))*CE(J))**2)
C
C NOTE THE UNITS IMBALANCE ABOVE
C
C WRITE (*,*) ' '//hm//' GE(J) ',J,' ',GE(J)
2000 CONTINUE
SIZE = 0.0
ERRSIZ = 0.0
DO 3000 I=1,6
SIZE = SIZE + G(I)*G(I)
ERRSIZ = ERRSIZ + GE(I)*GE(I)
3000 CONTINUE
SIZE = SQRT(SIZE)
ERRSIZ = SQRT(ERRSIZ)
RATIO = ERRSIZ / SIZE
END
C***********************************************************************
FUNCTION DOTVN (N,V1,V2)
C compute a dot product
REAL V1(N),V2(N)
C-----------------------------------------------------------------------
DOTVN = 0.0
DO 1000 I=1,N
DOTVN = DOTVN + V1(I)*V2(I)
1000 CONTINUE
END
C***********************************************************************
SUBROUTINE DRMV6 (V1,M,V2)
REAL V1(6),V2(6)
REAL M(36)
DOUBLE PRECISION SUM
C-----------------------------------------------------------------------
DO 3000 I=1,6
SUM = 0.0D0
DO 2000 J=1,6
SUM = SUM + DBLE(M(6*(I-1)+J))*DBLE(V1(J))
2000 CONTINUE
V2(I) = SUM
3000 CONTINUE
END
C***********************************************************************
logical function G6TOC (G,C,TEST)
C compute the normal unit cell parameters from a given g6 vector
include 'ITERATE.cmn'
CHARACTER *6 TEST
REAL G(6),C(6)
C-----------------------------------------------------------------------
IF (TEST .NE. 'G6TOC ') THEN
WRITE (*,*) ' '//hm//' TEST WAS WRONG IN G6TOC'
STOP
ENDIF
g6toc = .true.
DO 900 I=1,3
IF (G(I) .LE. 0.0) THEN
WRITE (*,*) ' '//hm//' G(I)<=0, I=',I,' ',G(I)
g6toc = .false.
ENDIF
900 CONTINUE
DO 1000 I=1,3
1000 C(I) = SQRT(G(I))
AC = 0.5*G(4)/C(2)/C(3)
IF (ABS(AC) .LE. 1.0) THEN
C(4) = 57.296*ACOS(AC)
ELSE
C(4) = 0.0
WRITE (*,*) ' '//hm//' ARG>1.0 C(4) ',AC,1.0-ABS(AC)
g6toc = .false.
ENDIF
AC = 0.5*G(5)/C(1)/C(3)
IF (ABS(AC) .LE. 1.0) THEN
C(5) = 57.296*ACOS(AC)
ELSE
C(5) = 0.0
WRITE (*,*) ' '//hm//' ARG>1.0 C(5) ',AC,1.0-ABS(AC)
g6toc = .false.
ENDIF
AC = 0.5*G(6)/C(1)/C(2)
IF (ABS(AC) .LE. 1.0) THEN
C(6) = 57.296*ACOS(AC)
ELSE
C(6) = 0.0
WRITE (*,*) ' '//hm//' ARG>1.0 C(6) ',AC,1.0-ABS(AC)
g6toc = .false.
ENDIF
END
C***********************************************************************
subroutine g6tor3 (g6,m3)
C compute the normal unit cell parameters from a given g6 vector
include 'ITERATE.cmn'
logical pcmnt_
integer i,j
real g6(6,6),m3(3,3)
C-----------------------------------------------------------------------
do 2000 i=1,3
do 1000 j=1,3
if (g6(i,j) .lt. -1.0e-6) then
if (ostyle.ne.'CIF ') then
write (*,*)
* ' '//hm//' negative element in upper left of g6'
else
cifres=pcmnt_(' negative square in g6 matrix')
endif
elseif (g6(i,j) .lt. 1.0e-6) then
m3(i,j) = 0.0
else
m3(i,j) = sqrt(g6(i,j))
endif
1000 continue
2000 continue
do 3000 i=1,3
call gtr3sn(m3(i,1),m3(i,2),m3(i,3), g6(i,4),g6(i,5),g6(i,6))
3000 continue
if (abs(g6(5,5)) .gt. 1.0e-6) then
if (g6(5,5)*(m3(1,1)*m3(3,3)+m3(1,3)*m3(3,1)) .lt. 0.0) then
do 4000 i= 1,3
m3(3,i) = -m3(3,i)
4000 continue
endif
elseif (abs(g6(5,4)) .gt. 1.0e-6) then
if (g6(5,4)*(m3(1,2)*m3(3,3)+m3(1,3)*m3(3,2)) .lt. 0.0) then
do 4100 i= 1,3
m3(3,i) = -m3(3,i)
4100 continue
endif
elseif (abs(g6(5,6)) .gt. 1.0e-6) then
if (g6(5,6)*(m3(1,1)*m3(3,2)+m3(1,2)*m3(3,1)) .lt. 0.0) then
do 4200 i= 1,3
m3(3,i) = -m3(3,i)
4200 continue
endif
endif
if (abs(g6(6,6)) .gt. 1.0e-6) then
if (g6(6,6)*(m3(1,1)*m3(2,2)+m3(1,2)*m3(2,1)) .lt. 0.0) then
do 5000 i= 1,3
m3(2,i) = -m3(2,i)
5000 continue
endif
elseif (abs(g6(6,5)) .gt. 1.0e-6) then
if (g6(6,5)*(m3(1,1)*m3(2,3)+m3(1,3)*m3(2,1)) .lt. 0.0) then
do 5100 i= 1,3
m3(2,i) = -m3(2,i)
5100 continue
endif
elseif (abs(g6(6,4)) .gt. 0.0) then
if (g6(6,4)*(m3(1,2)*m3(2,3)+m3(1,3)*m3(2,2)) .lt. 0.0) then
do 5200 i= 1,3
m3(2,i) = -m3(2,i)
5200 continue
endif
endif
end
subroutine gtr3sn(e1,e2,e3, g4,g5,g6)
if (e1 .ne. 0.0) then
e2 = unitsn(g6)*e2
e3 = unitsn(g5)*e3
elseif (e3 .ne. 0.0) then
e3 = unitsn(g4)*e3
endif
end
C***********************************************************************
LOGICAL FUNCTION INPCEL (LATSYM,CV,CE,eof)
C get the input lattice type, cell, and errors in the cell parameters
include 'ITERATE.cmn'
logical char_
logical numb_
EXTERNAL OKCELL
LOGICAL OKCELL
CHARACTER *1 LATSYM
REAL CV(6),CE(6)
logical eof
C-----------------------------------------------------------------------
eof = .false.
inpcel = .true.
cifeid = '.'
cifsgs = 'P'
1000 CONTINUE
IF (querst.ne.'NO')
* WRITE (*,*) ' '//hm//' Input Xtal Lattice Centering '
if (istyle.ne.'CIF ') then
READ (*,'(A1)',end=9000) LATSYM
else
cifres = char_('_cell.entry_id',cifeid)
if (.not.cifres) cifeid = '.'
cifres = char_('_cell.space_group_name_H-M',cifsgs)
LATSYM = 'P'
if(cifres) LATSYM=cifsgs(1:1)
if(cifsgs.eq.' ') cifsgs = 'P'
endif
IF (LATSYM .GE. 'a' .AND. LATSYM .LE. 'z')
2 LATSYM = CHAR(ICHAR(LATSYM)-ICHAR('a')+ICHAR('A'))
IF (querst.ne.'NO')
* WRITE (*,*) ' '//hm//' Input Cell Parameters'
if (istyle.ne.'CIF ') then
READ (*,*,end=9000) CV
else
do ii = 1,6
CE(II) = 0.
enddo
cifres = numb_('_cell.length_a',cv(1),ce(1))
if (.not.cifres) goto 9000
cifres = numb_('_cell.length_b',cv(2),ce(2))
cifres = numb_('_cell.length_c',cv(3),ce(3))
cifres = numb_('_cell.angle_alpha',cv(4),ce(4))
cifres = numb_('_cell.angle_beta',cv(5),ce(5))
cifres = numb_('_cell.angle_gamma',cv(6),ce(6))
endif
INPCEL = OKCELL(LATSYM,CV,' TALK')
IF (.NOT. INPCEL) GO TO 1000
IF (querst.ne.'NO')
* WRITE (*,*)
* ' '//hm//' Input Standard Deviations of Cell Parameters'
if (istyle.ne.'CIF ') then
READ (*,*,end=9000) CE
endif
DO II = 1,6
CE(II) = MAX(CE(II),ABS(CV(II))*5.E-7,1.E-4)
ENDDO
return
9000 eof=.true.
inpcel = .false.
END
C**********************************************************************C
FUNCTION INSPHR (MXTREE,NVEC,X,RADMAX,TREE,
2 MXLIST,NLIST,LIST,IDLIST,TEST)
C After the TREE has been constructed using BLDTRE,
C INSPHR is used to retrieve all of the points within
C RADMAX of the point X. MXTREE is the maximum size of
C TREE. The indices of the found points are returned in
C the array LIST; NLIST are returned, up to a maximum of
C MXLIST. For instance, TREE(LIST(3)) is the vector of
C the third point found in the list. IDLIST contains the
C corresponding list of the input ID's. If no points are
C found within RADMAX of X, then INSPHR and NLIST are
C returned as 0; otherwise they are returned as the
C number of points found. If more than MXLIST points were
C found, then NLIST is returned equal to MXLIST, and
C INSPHR is returned equal to -MXLIST. TEST must be the
C string 'INSPHR'.
C See also BLDTRE and NEARST.
include 'ITERATE.cmn'
CHARACTER*6 TEST
LOGICAL DEBUG
INTEGER ISTAK(1000)
REAL TREE(MXTREE)
INTEGER LIST(MXLIST),IDLIST(MXLIST)
INTEGER RMAX
DATA DEBUG /.FALSE./
DATA LINK,RMAX,ID,ICHILD /1,2,3,4/
DATA RIGHT,LEFT,END /111,112,113/
C----------------------------------------------------------------------C
IF (TEST .NE. 'INSPHR' .AND. TEST .NE. 'insphr')
2 CALL BADCAL (TEST,'INSPHR')
ISTKP = 0
NLIST = 0
IPOINT = 2
CURMIN = RADMAX
DIR = LEFT
DIRPRV = RIGHT
IPREV = IPOINT -1
if (tree(1).le.0.0) go to 8000
1000 CONTINUE
IF (IPREV .EQ. IPOINT .AND. DIRPRV .EQ. DIR) THEN
WRITE (*,*) ' '//hm//' INTERNAL ERROR IN INSPHR '
WRITE (*,*)
* ' '//hm//' TREE POINTER DIDN''T CHANGE',IPOINT,' ',DIR
STOP
ELSEIF (IPOINT .EQ. 0) THEN
WRITE (*,*)
* ' '//hm//' INTERNAL ERROR IN INSPHR, IPOINT = 0'
STOP
ENDIF
IPREV = IPOINT
DIRPRV = DIR
IF (DEBUG) WRITE (*,*)
* ' '//hm//' IN INSPHR 1000, IPOINT = ',IPOINT
IF (TREE(IPOINT) .EQ. 0) THEN
IF (DEBUG) WRITE (*,*)
* ' '//hm//' AT AN END WITH IPOINT = ',IPOINT
DIR = END
ELSEIF (DIR .EQ. RIGHT) THEN
IRIGHT = TREE(IPOINT)
IF (DEBUG) WRITE (*,*)
* ' '//hm//' WENT RIGHT WITH IPOINT ', IPOINT
DR = TREELN (NVEC,TREE(IRIGHT+ICHILD),X)
IF (DR .LT. CURMIN) THEN
NLIST = NLIST + 1
IF (NLIST .GT. MXLIST) GO TO 8000
LIST(NLIST) = IRIGHT + ICHILD
IDLIST(NLIST) = TREE(IRIGHT+ID)
ENDIF
IF (TREE(IRIGHT+LINK) .LE. 0) THEN
IF (DEBUG) WRITE (*,*)
* ' '//hm//' ON RIGHT BRANCH, UNSTACK A POINT'
DIR = END
ELSEIF (TREE(IRIGHT+RMAX)+CURMIN .GT. DR) THEN
IPOINT = TREE(IRIGHT+LINK)
DIR = LEFT
ELSE
DIR = END
ENDIF
ELSE
IF (DEBUG) WRITE (*,*) ' '//hm//' WENT LEFT, IPOINT ',IPOINT
DIR = LEFT
IF (TREE(IPOINT) .GT. 0) THEN
IF (DEBUG) WRITE (*,*) ' '//hm//' STACK ONE '
CALL TRSTCK (IPOINT,ISTAK,ISTKP)
ENDIF
DL = TREELN (NVEC,TREE(IPOINT+ICHILD),X)
IF (DL .LT. CURMIN) THEN
NLIST = NLIST + 1
IF (NLIST .GT. MXLIST) GO TO 8000
LIST(NLIST) = IPOINT+ICHILD
IDLIST(NLIST) = TREE(IPOINT+ID)
ENDIF
IF (TREE(IPOINT+LINK) .LE. 0) THEN
IF (DEBUG) WRITE (*,*) ' '//hm//' NO LEFT LINK, GO BACK'
DIR = END
ELSEIF (TREE(IPOINT+RMAX) .LT. 0.0) THEN
IF (DEBUG) WRITE (*,*) ' '//hm//' NO DESCENDING LEFT TREE'
DIR = END
ELSEIF (TREE(IPOINT+RMAX)+CURMIN .GT. DL) THEN
IF (DEBUG) WRITE (*,*) ' '//hm//
2 ' GOING TO GO DOWN ONE LEVEL ',
3 ' IPOINT AND UPDATE ',IPOINT,' ',INT(TREE(IPOINT+LINK))
IPOINT = TREE(IPOINT+LINK)
ELSE
IF (DEBUG) WRITE (*,*) ' '//hm//' NO CLOSER POINTS ON LEFT '
IF (DEBUG) WRITE (*,*) ' '//hm//' CURMIN,TREE(IPOINT+RMAX),
2 DL ',CURMIN,TREE(IPOINT+RMAX),DL
DIR = END
ENDIF
ENDIF
IF (DIR .EQ. END) THEN
IF (IUNSTK(IPOINT,ISTAK,ISTKP) .LE. 0) GO TO 8000
DIR = RIGHT
ENDIF
GO TO 1000
8000 CONTINUE
IF (NLIST .LE. MXLIST) THEN
INSPHR = NLIST
ELSE
INSPHR = -MXLIST
NLIST = MXLIST
ENDIF
END
C**********************************************************************C
FUNCTION IUNSTK(NEXT,ISTAK,ISTKP)
C helper function for INSPHR and NEARST
include 'ITERATE.cmn'
INTEGER ISTAK(1000)
LOGICAL DEBUG
DATA DEBUG /.FALSE./
C----------------------------------------------------------------------C
IF (DEBUG) WRITE (*,*) ' '//hm//' IUNSTK,NEXT ',NEXT
IF (ISTKP .GT. 0) THEN
NEXT = ISTAK(ISTKP)
ISTKP = ISTKP-1
ELSE
NEXT = 0
ENDIF
IUNSTK = NEXT
END
C***********************************************************************
SUBROUTINE MKNORM (VI,Mnorm,VOUT,TEST)
C converts an input g6 vector to "normalized" form (Gruber's
C terminology) or "standard presentation" (Andrews and Bernstein's
C terminology) and the corresponding transformation matrix
include 'ITERATE.cmn'
CHARACTER *6 TEST
REAL VIN(6),VI(6),VOUT(6)
real Mnorm(36),MAT(36),M1(36),MAT3(9)
LOGICAL AGAIN
C-----------------------------------------------------------------------
IF (TEST .NE. 'MKNORM') THEN
WRITE (*,*) ' '//hm//' TEST WAS WRONG IN MKNORM'
STOP
ENDIF
CALL CPYVN(6,VI,VIN)
CALL RUNTMN(6,Mnorm)
NCYCLE = 0
1000 CONTINUE
NCYCLE = NCYCLE + 1
AGAIN =.FALSE.
CALL ZEROS (36,MAT)
IF ( (ABS(VIN(1)).GT.ABS(VIN(2))) .OR.
2 (VIN(1).EQ.VIN(2) .AND. ABS(VIN(4)).GT.ABS(VIN(5))) ) THEN
MAT(2) = 1
MAT(7) = 1
MAT(15) = 1
MAT(23) = 1
MAT(28) = 1
MAT(36) = 1
AGAIN = .TRUE.
ELSEIF ( (ABS(VIN(2)).GT.ABS(VIN(3))) .OR.
2 (VIN(2).EQ.VIN(3) .AND. ABS(VIN(5)).GT.ABS(VIN(6))) ) THEN
MAT(1) = 1
MAT(14) = 1
MAT(9) = 1
MAT(22) = 1
MAT(35) = 1
MAT(30) = 1
AGAIN = .TRUE.
ENDIF
IF (AGAIN) THEN
CALL mm6(MAT,Mnorm,M1)
CALL CPYVN(36,M1,Mnorm)
CALL RMV6(VIN,MAT,VOUT)
CALL CPYVN(6,VOUT,VIN)
ENDIF
IF (AGAIN .AND. NCYCLE.LT.4) GO TO 1000
NUMNEG = 0
DO 2000 I=4,6
IF (VIN(I).LT.0.0) NUMNEG = NUMNEG + 1
2000 CONTINUE
CALL RUNTMN(3,MAT3)
DO 4000 I=4,6
IF (NUMNEG.EQ.1) THEN
IF(VIN(I).GE.0.) MAT(3*(I-4)+I-3) = -1.
C MAT(6*(I-1)+I) = -SIGN(1.0,VIN(I))
ELSEIF (NUMNEG.EQ.2) THEN
IF(VIN(I).LT.0.) MAT(3*(I-4)+I-3) = -1.
C MAT(6*(I-1)+I) = SIGN(1.0,VIN(I))
ENDIF
4000 CONTINUE
call r3tog6(MAT3,MAT)
CALL mm6(MAT,Mnorm,M1)
CALL CPYVN(36,M1,Mnorm)
CALL RMV6(VIN,MAT,VOUT)
CALL CPYVN(6,VOUT,VIN)
END
C***********************************************************************
SUBROUTINE MKPRIM (LATSYM,GIN,M,GOUT,TEST)
C converts and input g6 vector to one corresponding to a primitive
C lattice and the corresponding transformation matrix
include 'ITERATE.cmn'
CHARACTER *6 TEST
CHARACTER LATSYM
REAL GIN(6),GOUT(6)
REAL M(36)
C-----------------------------------------------------------------------
IF (TEST .NE. 'MKPRIM') THEN
WRITE (*,*) ' '//hm//' TEST WAS WRONG IN MKPRIM'
STOP
ENDIF
CALL ZEROS (36,M)
IF (LATSYM .EQ. 'P') THEN
CALL RUNTMN(6,M)
ELSEIF (LATSYM .EQ. 'I') THEN
M(1) = 1
M(8) = 1
DO 1000 I=13,18
1000 M(I) = 0.25
M(20) = 1
M(22) = 0.5
M(24) = 0.5
M(25) = 1
M(29) = 0.5
M(30) = 0.5
M(36) = 1
ELSEIF (LATSYM .EQ. 'F') THEN
M(1) = 0.25
M(2) = 0.25
M(6) = 0.25
M(7) = 0.25
M(9) = 0.25
M(11) = 0.25
M(14) = 0.25
M(15) = 0.25
M(16) = 0.25
M(21) = 0.5
DO 2000 I=22,24
2000 M(I)= 0.25
M(26) = 0.5
DO 2100 I=28,30
2100 M(I) = 0.25
M(31) = 0.5
DO 2200 I=34,36
2200 M(I) = 0.25
ELSEIF (LATSYM .EQ. 'A') THEN
M(1) = 1
M(8) = 1
M(14) = 0.25
M(15) = 0.25
M(16) = 0.25
M(20) = 1
M(22) = 0.5
M(29) = 0.5
M(30) = 0.5
M(36) = 1
ELSEIF (LATSYM .EQ. 'B') THEN
M(1) = 1
M(8) = 1
M(13) = 0.25
M(15) = 0.25
M(17) = 0.25
M(22) = 0.5
M(24) = 0.5
M(25) = 1
M(29) = 0.5
M(36) = 1
ELSEIF (LATSYM .EQ. 'C') THEN
M(1) = 1
M(7) = 0.25
M(8) = 0.25
M(12) = 0.25
M(15) = 1
M(22) = 0.5
M(23) = 0.5
M(29) = 1
M(31) = 1
M(36) = 0.5
ELSEIF (LATSYM .EQ. 'R') THEN
DO 2300 I=1,36
2300 M(I) = 1./9.
M(1) = 4./9.
M(5) = 2./9.
M(6) = 2./9.
M(11) = -1./9.
M(12) = -1./9.
M(14) = 4./9.
M(16) = -2./9.
M(17) = -1./9.
M(18) = 2./9.
M(19) = 2./9.
M(20) = -4./9.
M(21) = 2./9.
M(22) = -1./9.
M(23) = -2./9.
M(25) = -4./9.
M(26) = -4./9.
M(27) = 2./9.
M(28) = -1./9.
M(30) = -5./9.
M(31) = -4./9.
M(32) = 2./9.
M(33) = 2./9.
M(34) = 2./9.
ELSE
WRITE (*,*) ' '//hm//' DID NOT FIND LATTICE SYMBOL ',LATSYM
STOP
ENDIF
CALL RMV6(GIN,M,GOUT)
END
C CALL MKREFL (RATIO,MXTREE,TREE,NVMAX,V,MATREF,NV,GRED,'MKREFL')
C*********************************************************************** SUBROUTINE mm6(M1,M2,M3)
SUBROUTINE MKREFL
2 (DEBUG,RATIO,MXTREE,TREE,NVMAX,V,MATREF,NV,GRED,TEST)
C MKREFL performs the iterations to search out the various
C representations of a single lattice by a set of different unit cells.
C There are many ways to do this iteration, and this one may well
C not be optimal. It has been found to work well in practice as long
C as cutoffs are not too strict. It is clear that, in general, it does
C not find all of the possible unit cells within a particular radius
C in g6.
include 'ITERATE.cmn'
PARAMETER (MXSWTC=2)
CHARACTER *6 TEST
REAL TREE(MXTREE), V(6,NVMAX), GRED(6)
REAL VT(6),VTT(6)
real mt2k(36),mti2k(36),mt1k(36)
real REFL(36,24)
real SWTCH(36,MXSWTC)
real MATREF(36,NVMAX)
LOGICAL DEBUG
C-----------------------------------------------------------------------
DATA (SWTCH(I,1),I=1,36)/
1 1,0,0, 0,0,0,
2 0,1,0, 0,0,0,
3 1,0,1, 0,1,0,
4 0,0,0, 1,0,1,
5 2,0,0, 0,1,0,
6 0,0,0, 0,0,1 /
DATA (SWTCH(I,2),I=1,36)/
1 1,0,0, 0,0,0,
2 0,1,0, 0,0,0,
3 1,1,1, 1,1,1,
4 0,2,0, 1,0,1,
5 2,0,0, 0,1,1,
6 0,0,0, 0,0,1 /
DATA (REFL(I,1),I=1,36) /
2 1,0,0,0,0,0,
3 0,1,0,0,0,0,
4 0,0,1,0,0,0,
5 0,0,0,1,0,0,
6 0,0,0,0,1,0,
7 0,0,0,0,0,1 /
DATA (REFL(I,2),I=1,36) /
2 1,0,0,0,0,0,
3 0,1,0,0,0,0,
4 0,0,1,0,0,0,
5 0,0,0,-1,0,0,
6 0,0,0,0,-1,0,
7 0,0,0,0,0,1 /
DATA (REFL(I,3),I=1,36) /
2 1,0,0,0,0,0,
3 0,1,0,0,0,0,
4 0,0,1,0,0,0,
5 0,0,0,1,0,0,
6 0,0,0,0,-1,0,
7 0,0,0,0,0,-1 /
DATA (REFL(I,4),I=1,36) /
2 1,0,0,0,0,0,
3 0,1,0,0,0,0,
4 0,0,1,0,0,0,
5 0,0,0,-1,0,0,
6 0,0,0,0,1,0,
7 0,0,0,0,0,-1 /
DATA (REFL(I,5),I=1,36) /
2 0,1,0,0,0,0,
3 1,0,0,0,0,0,
4 0,0,1,0,0,0,
5 0,0,0,0,1,0,
6 0,0,0,1,0,0,
7 0,0,0,0,0,1 /
DATA (REFL(I,6),I=1,36) /
2 0,1,0,0,0,0,
3 1,0,0,0,0,0,
4 0,0,1,0,0,0,
5 0,0,0,0,-1,0,
6 0,0,0,-1,0,0,
7 0,0,0,0,0,1 /
DATA (REFL(I,7),I=1,36) /
2 0,1,0,0,0,0,
3 1,0,0,0,0,0,
4 0,0,1,0,0,0,
5 0,0,0,0,-1,0,
6 0,0,0,1,0,0,
7 0,0,0,0,0,-1 /
DATA (REFL(I,8),I=1,36) /
2 0,1,0,0,0,0,
3 1,0,0,0,0,0,
4 0,0,1,0,0,0,
5 0,0,0,0,1,0,
6 0,0,0,-1,0,0,
7 0,0,0,0,0,-1 /
DATA (REFL(I,9),I=1,36) /
2 1,0,0,0,0,0,
3 0,0,1,0,0,0,
4 0,1,0,0,0,0,
5 0,0,0,1,0,0,
6 0,0,0,0,0,1,
7 0,0,0,0,1,0 /
DATA (REFL(I,10),I=1,36) /
2 1,0,0,0,0,0,
3 0,0,1,0,0,0,
4 0,1,0,0,0,0,
5 0,0,0,-1,0,0,
6 0,0,0,0,0,-1,
7 0,0,0,0,1,0 /
DATA (REFL(I,11),I=1,36) /
2 1,0,0,0,0,0,
3 0,0,1,0,0,0,
4 0,1,0,0,0,0,
5 0,0,0,-1,0,0,
6 0,0,0,0,0,1,
7 0,0,0,0,-1,0 /
DATA (REFL(I,12),I=1,36) /
2 1,0,0,0,0,0,
3 0,0,1,0,0,0,
4 0,1,0,0,0,0,
5 0,0,0,1,0,0,
6 0,0,0,0,0,-1,
7 0,0,0,0,-1,0 /
DATA (REFL(I,13),I=1,36) /
2 0,0,1,0,0,0,
3 0,1,0,0,0,0,
4 1,0,0,0,0,0,
5 0,0,0,0,0,1,
6 0,0,0,0,1,0,
7 0,0,0,1,0,0 /
DATA (REFL(I,14),I=1,36) /
2 0,0,1,0,0,0,
3 0,1,0,0,0,0,
4 1,0,0,0,0,0,
5 0,0,0,0,0,-1,
6 0,0,0,0,-1,0,
7 0,0,0,1,0,0 /
DATA (REFL(I,15),I=1,36) /
2 0,0,1,0,0,0,
3 0,1,0,0,0,0,
4 1,0,0,0,0,0,
5 0,0,0,0,0,-1,
6 0,0,0,0,1,0,
7 0,0,0,-1,0,0 /
DATA (REFL(I,16),I=1,36) /
2 0,0,1,0,0,0,
3 0,1,0,0,0,0,
4 1,0,0,0,0,0,
5 0,0,0,0,0,1,
6 0,0,0,0,-1,0,
7 0,0,0,-1,0,0 /
DATA (REFL(I,17),I=1,36) /
2 0,1,0,0,0,0,
3 0,0,1,0,0,0,
4 1,0,0,0,0,0,
5 0,0,0,0,1,0,
6 0,0,0,0,0,1,
7 0,0,0,1,0,0 /
DATA (REFL(I,18),I=1,36) /
2 0,1,0,0,0,0,
3 0,0,1,0,0,0,
4 1,0,0,0,0,0,
5 0,0,0,0,-1,0,
6 0,0,0,0,0,-1,
7 0,0,0,1,0,0 /
DATA (REFL(I,19),I=1,36) /
2 0,1,0,0,0,0,
3 0,0,1,0,0,0,
4 1,0,0,0,0,0,
5 0,0,0,0,-1,0,
6 0,0,0,0,0,1,
7 0,0,0,-1,0,0 /
DATA (REFL(I,20),I=1,36) /
2 0,1,0,0,0,0,
3 0,0,1,0,0,0,
4 1,0,0,0,0,0,
5 0,0,0,0,1,0,
6 0,0,0,0,0,-1,
7 0,0,0,-1,0,0 /
DATA (REFL(I,21),I=1,36) /
2 0,0,1,0,0,0,
3 1,0,0,0,0,0,
4 0,1,0,0,0,0,
5 0,0,0,0,0,1,
6 0,0,0,1,0,0,
7 0,0,0,0,1,0 /
DATA (REFL(I,22),I=1,36) /
2 0,0,1,0,0,0,
3 1,0,0,0,0,0,
4 0,1,0,0,0,0,
5 0,0,0,0,0,-1,
6 0,0,0,-1,0,0,
7 0,0,0,0,1,0 /
DATA (REFL(I,23),I=1,36) /
2 0,0,1,0,0,0,
3 1,0,0,0,0,0,
4 0,1,0,0,0,0,
5 0,0,0,0,0,-1,
6 0,0,0,1,0,0,
7 0,0,0,0,-1,0 /
DATA (REFL(I,24),I=1,36) /
2 0,0,1,0,0,0,
3 1,0,0,0,0,0,
4 0,1,0,0,0,0,
5 0,0,0,0,0,1,
6 0,0,0,-1,0,0,
7 0,0,0,0,-1,0 /
IF (TEST .NE. 'MKREFL') THEN
WRITE (*,*) ' '//hm//' TEST WAS WRONG IN MKREFL'
STOP
ENDIF
TREE(1) = 0
TREE(2) = 0
NV = 0
DMIN = 0.1*SQRT(DOTVN(6,GRED,GRED))*RATIO
GRMIN = GRED(1)**2+GRED(2)**2+GRED(3)**2
DO 4000 ICYCLE=1,MXSWTC+1
DO 2000 IREFL=1,24
CALL RMV6(GRED,REFL(1,IREFL),VT)
C#
C#
C accumulate matrices starting from the reduced vector
call cpyvn (36,refl(1,irefl),mt2k)
C#
C#
C There is a theoretical upper bound on how much the sum of the
C edge lengths can change (3.0) in a single transformation. So
C by allowing their square to only change by a factor of 9.0
C (incremented to allow for error), we can limit the search.
IF (VT(1)**2+VT(2)**2+VT(3)**2.GT. 10.0*GRMIN) THEN
IF (DEBUG)WRITE (*,*) ' '//hm//' REJECT 1 ',GRMIN,VT
ELSEIF (ABS(VT(4)**2/VT(2)/VT(3)) .GT. 1.2 .OR.
2 ABS(VT(5)**2/VT(1)/VT(3)) .GT. 1.2 .OR.
3 ABS(VT(6)**2/VT(1)/VT(2)) .GT. 1.2 ) THEN
IF (DEBUG) WRITE (*,*)
* ' '//hm//' REJECT 2 ',VT(4)**2/VT(2)/VT(3),
2 VT(5)**2/VT(1)/VT(3),
3 VT(6)**2/VT(1)/VT(2)
ELSE
C set things up so that on the first cycle, the original vector
C is stored, and nothing else happens. after that, the rest of the
C switch matrices are used (from the 4000 loop).
IF (ICYCLE.EQ.1) THEN
MXINNR = 1
C#
C#
call cpyvn (36,mt2k,mti2k)
C#
C#
ELSE
MXINNR = 24
CALL RMV6(VT,SWTCH(1,ICYCLE-1),VTT)
C#
C#
call mm6(swtch(1,icycle-1),mt2k,mti2k)
C#
C#
CALL CPYVN(6,VTT,VT)
ENDIF
DO 1000 INNER=1,MXINNR
IF (ICYCLE.GT.1) THEN
CALL RMV6(VT,REFL(1,INNER),VTT)
CALL CPYVN(6,VTT,VT)
C#
C#
call mm6(refl(1,inner),mti2k,mt1k)
call cpyvn (36,mt1k,mti2k)
C#
C#
else
call cpyvn (36,mti2k,mt1k)
ENDIF
NV1 = NV + 1
IF (NEARST(DEBUG,
2 MXTREE,6,VT,DMIN,TREE,NV1,ID,'NEARST') .EQ. 0)
3 THEN
NV = NV + 1
CALL BLDTRE (MXTREE,6,VT,NV,TREE,'BLDTRE')
CALL CPYVN(6,VT,V(1,NV))
C#
C#
call cpyvn(36,mt1k,matref(1,nv))
C#
C#
IF (NV .EQ. NVMAX) GO TO 8000
ENDIF
1000 CONTINUE
ENDIF
2000 CONTINUE
4000 CONTINUE
8000 CONTINUE
if (debug) then
WRITE (*,*) ' '//hm//' NV IN MKREFL ',NV
endif
END
C***********************************************************************
SUBROUTINE mm6(M1,M2,M3)
C multiply two matrices (6x6), both in a linear array
REAL M1(36),M2(36),M3(36)
C-----------------------------------------------------------------------
CALL ZEROS (36,M3)
DO 3000 I36=1,36
IROW = (I36+5)/6
ICOL = MOD(I36-1,6)+1
K = ICOL-6
DO 2000 J=6*IROW-5,6*IROW
K = K + 6
M3(I36) = M3(I36) + M1(J)*M2(K)
1000 CONTINUE
2000 CONTINUE
3000 CONTINUE
END
C***********************************************************************
LOGICAL FUNCTION NEARRD (A,R,SIG)
C test whether a cell (really a g6 vector) is nearly Buerger reduced
include 'ITERATE.cmn'
REAL A(6),R(6), B(6)
LOGICAL DEBUG
DATA DEBUG /.FALSE./
C R IS THE REDUCED CELL VECTOR
C check that the vector is near to reduced
C-----------------------------------------------------------------------
DO 1000 I=1,3
IF (A(I) .LT. 1.0) THEN
NEARRD = .FALSE.
RETURN
ENDIF
1000 CONTINUE
DO 2000 I=1,6
B(I) = A(I)
2000 CONTINUE
IF (B(1) .GT. B(2)) THEN
BT = B(1)
B(1) = B(2)
B(2) = BT
BT = B(4)
B(4) = B(5)
B(5) = BT
ENDIF
IF (B(2) .GT. B(3)) THEN
BT = B(2)
B(2) = B(3)
B(3) = BT
BT = B(5)
B(5) = B(6)
B(6) = BT
ENDIF
IF (B(1) .GT. B(2)) THEN
BT = B(1)
B(1) = B(2)
B(2) = BT
BT = B(4)
B(4) = B(5)
B(5) = BT
ENDIF
NEARRD = .TRUE.
DO 4000 I=1,3
IF (B(I)-R(I) .GT. 5.0*SIG) THEN
IF (DEBUG) WRITE (*,*) ' '//hm//
2 ' NOT NEAR BUERGER ',I,B(I),R(I)
NEARRD = .FALSE.
RETURN
ENDIF
4000 CONTINUE
IF (ABS(B(4))/SQRT(B(2)*B(3)) .GT. 1.0+5.0*SIG) THEN
IF (DEBUG) WRITE (*,*)' '//hm//' ALPHA BAD '
NEARRD = .FALSE.
RETURN
ELSEIF (ABS(B(5))/SQRT(B(1)*B(3)) .GT. 1.0+5.0*SIG) THEN
IF (DEBUG) WRITE (*,*)' '//hm//' BETA BAD '
NEARRD = .FALSE.
RETURN
ELSEIF (ABS(B(6))/SQRT(B(1)*B(2)) .GT. 1.0+5.0*SIG) THEN
IF (DEBUG) WRITE (*,*)' '//hm//' GAMMA BAD '
NEARRD = .FALSE.
RETURN
ENDIF
END
C**********************************************************************C
FUNCTION NEARST (debug,MXTREE,NVEC,X,RADMAX,TREE,IP,IDOUT,TEST)
C The parameters are the same as those of INSPHR, except
C that only the nearest point is found. If no points are
C found within RADMAX of X, then IP and NEARST are
C returned as 0; otherwise they are returned as the index
C in TREE of the coordinates of the vector of the point
C closest to X. TEST must be the string 'NEARST'. ID is
C returned as the ID of the nearest point (if the value
C input as IDIN in BLDTRE was the index)
C See also BLDTRE and INSPHR
include 'ITERATE.cmn'
CHARACTER*6 TEST
LOGICAL DEBUG
parameter (maxstk=1000)
INTEGER ISTAK(maxstk)
REAL TREE(MXTREE)
C DATA DEBUG /.FALSE./
DATA LINK,RMAX,ID,ICHILD /1,2,3,4/
DATA RIGHT,LEFT,END /111,112,113/
C----------------------------------------------------------------------C
IF (TEST .NE. 'NEARST' .AND. TEST .NE. 'nearst')
2 CALL BADCAL (TEST,'NEARST')
ISTKP = 0
IP = 0
IPOINT = 2
CURMIN = RADMAX
DIR = LEFT
DIRPRV = RIGHT
IPREV = IPOINT -1
if (tree(1).le.0.0) go to 8000
1000 CONTINUE
IPREV = IPOINT
DIRPRV = DIR
IF (DEBUG) WRITE (*,*) ' '//hm//
2 ' IN NEARST 1000, IPOINT = ',IPOINT
IF (TREE(IPOINT) .EQ. 0) THEN
IF (DEBUG) WRITE (*,*) ' '//hm//
2 ' AT AN END WITH IPOINT = ',IPOINT
DIR = END
ELSEIF (DIR .EQ. RIGHT) THEN
IRIGHT = TREE(IPOINT)
IF (DEBUG) WRITE (*,*) ' '//hm//
2 ' WENT RIGHT WITH IPOINT ', IPOINT
DR = TREELN (NVEC,TREE(IRIGHT+ICHILD),X)
IF (DR .LT. CURMIN) THEN
CURMIN = DR
IP = IRIGHT + ICHILD
IDOUT = TREE(IRIGHT+ID)
ENDIF
IF (TREE(IRIGHT+LINK) .LE. 0) THEN
IF (DEBUG) WRITE (*,*) ' '//hm//
2 ' ON RIGHT BRANCH, UNSTACK A POINT'
DIR = END
ELSEIF (TREE(IRIGHT+RMAX)+CURMIN .GT. DR) THEN
IPOINT = TREE(IRIGHT+LINK)
DIR = LEFT
ELSE
DIR = END
ENDIF
ELSE
IF (DEBUG) WRITE (*,*) ' '//hm//' WENT LEFT, IPOINT ',IPOINT
DIR = LEFT
IF (TREE(IPOINT) .GT. 0) THEN
IF (DEBUG) WRITE (*,*) ' '//hm//' STACK ONE '
CALL TRSTCK(IPOINT,ISTAK,ISTKP)
ENDIF
DL = TREELN (NVEC,TREE(IPOINT+ICHILD),X)
IF (DL .LT. CURMIN) THEN
CURMIN = DL
IP = IPOINT+ICHILD
IDOUT = TREE(IPOINT+ID)
ENDIF
IF (TREE(IPOINT+LINK) .LE. 0) THEN
IF (DEBUG) WRITE (*,*) ' '//hm//' NO LEFT LINK, GO BACK'
DIR = END
ELSEIF (TREE(IPOINT+RMAX) .LT. 0.0) THEN
IF (DEBUG) WRITE (*,*) ' '//hm//' NO DESCENDING LEFT TREE'
DIR = END
ELSEIF (TREE(IPOINT+RMAX)+CURMIN .GT. DL) THEN
IF (DEBUG) WRITE (*,*) ' '//hm//
2 ' GOING TO GO DOWN ONE LEVEL ',
3 ' IPOINT AND UPDATE ',IPOINT,' ',INT(TREE(IPOINT+LINK))
IPOINT = TREE(IPOINT+LINK)
ELSE
IF (DEBUG) WRITE (*,*) ' '//hm//' NO CLOSER POINTS ON LEFT '
IF (DEBUG) WRITE (*,*) ' '//hm//' CURMIN,TREE(IPOINT+RMAX),
2 DL ',CURMIN,TREE(IPOINT+RMAX),DL
DIR = END
ENDIF
ENDIF
IF (DIR .EQ. END) THEN
if (debug) write (*,*) ' '//hm//
2 ' call unstk ',ipoint,istkp,maxstk
IF (IUNSTK(IPOINT,ISTAK,ISTKP) .LE. 0) GO TO 8000
DIR = RIGHT
ENDIF
GO TO 1000
8000 CONTINUE
NEARST = IP
END
C***********************************************************************
LOGICAL FUNCTION OKCELL (LATSYM,CV,talk)
C okcell determines if a particular cell makes sense. The variable
C talk is used to determine if output is printed.
include 'ITERATE.cmn'
logical pcmnt_
PARAMETER (NLATT=16)
CHARACTER *1 SYMLST(NLATT)
CHARACTER *1 LATSYM
REAL CV(6)
character *(*) talk
DATA SYMLST /'P','A','B','C','I','F','R','H',
2 'p','a','b','c','i','f','r','h'/
C-----------------------------------------------------------------------
OKCELL = .TRUE.
IF (CV(4).GT.175.0 .OR. CV(5).GT.175.0 .OR. CV(6).GT.175.0) THEN
if (talk .eq. ' TALK')
2 WRITE (*,*)
* ' '//hm//' THE LATTICE ANGLES MUST BE LESS THAN 175 DEGREES'
OKCELL = .FALSE.
ENDIF
IF (CV(4).LT.5.0 .OR. CV(5).LT.5.0 .OR. CV(6).LT.5.0) THEN
if (talk .eq. ' TALK')
2 WRITE (*,*)
* ' '//hm//' THE LATTICE ANGLES MUST EXCEED 5.0 DEGREES'
OKCELL = .FALSE.
ENDIF
DO 1000 I=1,NLATT
IF (LATSYM.EQ.SYMLST(I)) THEN
GO TO 1100
ENDIF
1000 CONTINUE
if (talk .eq. ' TALK')
2 WRITE (*,*)
* ' '//hm//' XTAL CENTERING SYMBOL ',LATSYM,' IS NOT IMPLEMENTED'
OKCELL = .FALSE.
1100 CONTINUE
IF (CV(4).GE.CV(5)+CV(6)) THEN
if (talk .eq. ' TALK')
2 WRITE (*,*)
* ' '//hm//' ERROR, ALPHA EXCEEDS BETA PLUS GAMMA'
OKCELL = .FALSE.
ENDIF
IF (CV(5).GE.CV(4)+CV(6)) THEN
if (talk .eq. ' TALK')
2 WRITE (*,*) ' '//hm//' ERROR, BETA EXCEEDS ALPHA PLUS GAMMA'
OKCELL = .FALSE.
ENDIF
IF (CV(6).GE.CV(4)+CV(5)) THEN
if (talk .eq. ' TALK')
2 WRITE (*,*) ' '//hm//' ERROR, GAMMA EXCEEDS ALPHA PLUS BETA'
OKCELL = .FALSE.
ENDIF
IF (CV(4)+CV(5)+CV(6) .GT. 355.0) THEN
if (talk .eq. ' TALK') then
WRITE (*,*) ' '//hm//
2 ' THE SUM OF THE CELL ANGLES MUST BE LESS THAN'
WRITE (*,*) ' '//hm//' 355.0 DEGREES'
endif
OKCELL = .FALSE.
ENDIF
DO 2000 I=1,3
IF (CV(I) .LT. 1.0) THEN
WRITE (*,*) ' '//hm//' THE AXIAL LENGTHS MUST EXCEED 1.0'
OKCELL = .FALSE.
ENDIF
2000 CONTINUE
IF (LATSYM .EQ. 'R' .OR. LATSYM .EQ. 'r' .OR.
* LATSYM .EQ. 'H' .OR. LATSYM .EQ. 'h') THEN
IF (ABS(CV(4)-90.).GT.1. .OR.
* ABS(CV(5)-90.).GT.1. .OR.
* ABS(CV(6)-120.).GT.1. .OR.
* ABS(CV(1)-CV(2)).GT..5) THEN
IF (LATSYM.EQ.'R' .OR. LATSYM.EQ.'r') THEN
if (talk .eq. ' TALK') then
if (ostyle.ne.'CIF ') then
WRITE (*,*) ' '//hm//' RHOMBOHEDRAL CENTERING REQUIRES'
WRITE (*,*) ' '//hm//
2 ' A HEXAGONAL CELL (A=B, ALPHA=BETA=90,'
WRITE (*,*) ' '//hm//' GAMMA =120)'
else
cifres=pcmnt_(' Non-hexagonal R lattice '//cifeid)
endif
endif
IF (ABS(CV(1)-CV(2)) .LT. .5 .AND.
* ABS(CV(2)-CV(3)) .LT. .5 .AND.
* ABS(CV(1)-CV(3)) .LT. .5 .AND.
* ABS(CV(4)-CV(5)) .LT. 1. .AND.
* ABS(CV(5)-CV(6)) .LT. 1. .AND.
* ABS(CV(4)-CV(6)) .LT. 1.) THEN
LATSYM = 'P'
if (ostyle.ne.'CIF ') then
WRITE (*,*) ' '//hm//
2 ' PROCESSING AS PRMITIVE RHOMBOHEDRAL'
else
cifres=pcmnt_(' Processing as primitve rhombohedral')
endif
ELSE
OKCELL = .FALSE.
ENDIF
ELSE
if (talk .eq. ' TALK') then
if (ostyle.ne.'CIF ') then
WRITE (*,*) ' '//hm//' A HEXAGONAL CELL REQUIRES'
WRITE (*,*) ' '//hm//' A=B, ALPHA=BETA=90,GAMMA =120'
else
cifres=pcmnt_(' non-hexagonal cell '//cifeid)
endif
endif
OKCELL = .FALSE.
ENDIF
ENDIF
IF (LATSYM .EQ. 'H' .OR. LATSYM .EQ. 'h') THEN
WRITE(*,*) ' '//hm//' PROCESSING H AS P'
LATSYM = 'P'
ENDIF
ENDIF
RETURN
END
subroutine r3tog6(e3,g6)
C give a 3-space transformation matrix, determine the corresponding one
C in g6
real e3(9),g6(36)
C-----------------------------------------------------------------------
C------------------------
C upper left 3x3 block
C------------------------
g6(1)=e3(1)**2
g6(2)=e3(2)**2
g6(3)=e3(3)**2
g6(7)=e3(4)**2
g6(8)=e3(5)**2
g6(9)=e3(6)**2
g6(13)=e3(7)**2
g6(14)=e3(8)**2
g6(15)=e3(9)**2
C------------------------
C upper right 3x3 block
C------------------------
g6(4)=e3(2)*e3(3)
g6(5)=e3(1)*e3(3)
g6(6)=e3(1)*e3(2)
g6(10)=e3(5)*e3(6)
g6(11)=e3(4)*e3(6)
g6(12)=e3(4)*e3(5)
g6(16)=e3(8)*e3(9)
g6(17)=e3(7)*e3(9)
g6(18)=e3(7)*e3(8)
C------------------------
C lower left 3x3 block
C------------------------
g6(19)=2*e3(4)*e3(7)
g6(20)=2*e3(5)*e3(8)
g6(21)=2*e3(6)*e3(9)
g6(25)=2*e3(1)*e3(7)
g6(26)=2*e3(2)*e3(8)
g6(27)=2*e3(3)*e3(9)
g6(31)=2*e3(1)*e3(4)
g6(32)=2*e3(2)*e3(5)
g6(33)=2*e3(3)*e3(6)
C------------------------
C lower right 3x3 block
C------------------------
g6(22)=e3(5)*e3(9) + e3(8)*e3(6)
g6(23)=e3(4)*e3(9) + e3(7)*e3(6)
g6(24)=e3(4)*e3(8) + e3(4)*e3(5)
g6(28)=e3(2)*e3(9) + e3(8)*e3(3)
g6(29)=e3(1)*e3(9) + e3(7)*e3(3)
g6(30)=e3(1)*e3(8) + e3(7)*e3(2)
g6(34)=e3(2)*e3(6) + e3(5)*e3(3)
g6(35)=e3(1)*e3(6) + e3(4)*e3(3)
g6(36)=e3(1)*e3(5) + e3(4)*e3(2)
end
C***********************************************************************
PROGRAM RED
include 'ITERATE.cmn'
logical init_
logical data_
logical ocif_
logical pdata_
logical pfile_
logical ploop_
logical pnumb_
logical pcmnt_
logical pchar_
INTEGER MAXPRJ
PARAMETER (MAXPRJ=42)
PARAMETER (NVMAX=1000)
PARAMETER (MXTREE=11*NVMAX)
INTEGER ITDESG(MAXPRJ)
CHARACTER *2 CHRLAT(MAXPRJ)
REAL PJNORM(MAXPRJ)
C common /cmprjn/ pjnorm
real PRJ(36,MAXPRJ)
C common /cmprj/ prj
REAL P(36),AP(36)
REAL CV(6),CE(6),G(6),GE(6),TG(6),AG(6),COUT(6)
REAL MPRIM(36)
real m3ptrd(9)
real MRED(36)
REAL GOUT(6),GRED(6),CRED(6)
LOGICAL INPCEL
LOGICAL NEARRD
EXTERNAL INPCEL
CHARACTER LATSYM
REAL TREE(MXTREE)
C common /cmtree/ tree
REAL V(6,NVMAX)
C common /cmv/ v
real MATREF(36,NVMAX),m1(36),m2(36)
real retcel(6),m3(9),m3t(9),m3ti(9)
REAL VBEST(6),AVBEST(6)
REAL DOTVN
EXTERNAL DOTVN
INTEGER NPROJ,I,J
REAL RATIO,SIZE,ERRSIZ
character*8 cellst(6)
logical eof,dbgmkr
logical g6toc, test
logical debug
data debug /.FALSE./
C-----------------------------------------------------------------------
CALL GETENV('ITERATE_QUERY',querst)
CALL GETENV('OUTPUT_STYLE',ostyle)
CALL GETENV('INPUT_STYLE',istyle)
if(istyle.eq.'CIF '.or.ostyle.eq.'CIF ') then
cifres = init_(5,6,21,0)
endif
if(istyle.eq.'CIF ') then
querst = 'NO'
cifres = ocif_(' ')
cifres = data_(' ')
endif
hm = ' '
if (ostyle.eq.'CIF ') then
hm = '#'
tabl_ = .false.
cifres = pfile_(' ')
cifres = pdata_('G6_SEARCH')
cifres = ploop_('_cell.entry_id')
cifres = ploop_('_cell.id')
cifres = ploop_('_cell.space_group_name_H-M')
cifres = ploop_('_cell.Bravais_lattice_symbol')
cifres = ploop_('_cell.length_a')
cifres = ploop_('_cell.length_b')
cifres = ploop_('_cell.length_c')
cifres = ploop_('_cell.angle_alpha')
cifres = ploop_('_cell.angle_beta')
cifres = ploop_('_cell.angle_gamma')
cifres = ploop_('_cell.unreduced_length_a')
cifres = ploop_('_cell.unreduced_length_b')
cifres = ploop_('_cell.unreduced_length_c')
cifres = ploop_('_cell.unreduced_angle_alpha')
cifres = ploop_('_cell.unreduced_angle_beta')
cifres = ploop_('_cell.unreduced_angle_gamma')
endif
do 1 i=2,mxtree
tree(i) = 2**30
1 continue
iunt0 = 0
iunt1 = 1
iunt2 = 2
iunt3 = 3
iunt10 = 10
dbgmkr = .false.
IF (querst.ne.'NO') THEN
WRITE (*,*) ' '//hm//' BEFORE BLDPRJ'
ENDIF
CALL BLDPRJ (MAXPRJ,NPROJ,ITDESG,CHRLAT,PJNORM,PRJ,'BLDPRJ')
C write (*,*) ' '//hm//' nproj = ',nproj
C loop until the input contains no more data
1000 continue
IF (INPCEL(LATSYM,CV,CE,eof)) THEN
CALL CTOG6(CV,CE,G,GE,SIZE,ERRSIZ,RATIO,'CTOG6 ')
if (ostyle.ne.'CIF ') then
WRITE (*,'(A,A)') ' '//hm//' Input Lattice Symbol ',
2 latsym
WRITE (*,*)
WRITE (*,'('' '//hm//' INPUT CELL AND ERRORS '',
2 '' INPUT VECTOR AND ERRORS'')')
WRITE (*,*)
else
cifres = pchar_(' ',cifeid)
cifres = pchar_(' ','.')
cifres = pchar_(' ',cifsgs)
cifres = pchar_(' ',latsym)
endif
CALL WRCLV6 (CV,CE,G,GE,'WRCLV6')
if (ostyle.ne.'CIF ') then
WRITE (*,*)
else
do ii = 1,6
cifres = pchar_(' ','.')
enddo
endif
CALL WRSIZE (SIZE,ERRSIZ,RATIO,'WRSIZE')
CALL MKPRIM (LATSYM,G,MPRIM,GOUT,'MKPRIM')
call cpyvn (36,mprim,m1)
CALL CHKVEC(GOUT)
CALL RUNTMN(6,MRED)
CALL REDUCE (GOUT,MRED,GRED,'REDUCE')
call mm6 (mred,m1,m2)
CALL CHKVEC(GRED)
test = G6TOC (GRED,CRED,'G6TOC ')
SIZE = SQRT(DOTVN(6,GRED,GRED))
ERRSIZ = RATIO * SIZE
if (ostyle.ne.'CIF ') then
WRITE (*,*)
WRITE (*,'('' '//hm//' REDUCED CELL '')')
WRITE (*,'('' '//hm//' Red. Cell '',6F10.3)')CRED
WRITE (*,'('' '//hm//' Red. Vector '',6F9.2)') GRED
endif
CALL MKREFL (dbgmkr,RATIO,MXTREE,TREE,NVMAX,
2 V,MATREF,NV,GRED,'MKREFL')
if (querst .ne. 'NO') then
if(ostyle.ne.'CIF ')write (*,*) ' '//hm//' nv=',nv
endif
do 1200 iv=1,nv
call mm6(matref(1,iv),m2,m1)
call cpyvn (36,m1,matref(1,iv))
1200 continue
if (debug) then
do 1300 iv=1,nv
write (*,'(1x,a3,i5,6f8.3)')
* ' '//hm//' ',iv,(v(ip,iv),ip=1,6)
1300 continue
endif
C after reducing the input cell and iterating to find various
C nearly reduced cells, test the found cells using the projectors
C of Paciorek and Bonin
nmatch = 0
DO 4000 I=1,NPROJ
nrej = 0
DBEST = 1.0E20
DBESTO = DBEST
nbest = 0
C subtract the projector from the unit matrix (to give the "prep")
C which when multiplied times a vector gives the vector component
C not in the subspace defined by the projector
DO 3100 J=1,36
P(J) = REAL(PRJ(J,I))/PJNORM(I)
3100 AP(J) = -P(J)
DO 3200 J=1,36,7
AP(J) = 1.0 + AP(J)
3200 CONTINUE
DO 3300 IV=1,NV
CALL RMV6 (V(1,IV),AP,AG)
CALL RMV6 (V(1,IV),P,TG)
DTEST = DOTVN(6,AG,AG)/dotvn(6,tg,tg)
IF (DTEST .GE. DBEST) THEN
ELSEIF (TG(1) .LT. 1.0) THEN
ELSEIF (TG(2) .LT. 1.0) THEN
ELSEIF (TG(3) .LT. 1.0) THEN
ELSEIF (TG(1)*TG(2)*TG(3) + 0.125*TG(6)*TG(4)*TG(5)
2 + 0.125*TG(5)*TG(6)*TG(4) - 0.25*TG(5)*TG(2)*TG(5)
3 - 0.25*TG(1)*TG(4)*TG(4) - 0.25*TG(6)*TG(6)*TG(3)
4 .LE. 0.0) THEN
C if the metric tensor is negative, the cell is bad
NREJ = NREJ + 1
IF (DEBUG) WRITE (*,*) ' '//hm//' FAILED METRIC TENSOR'
ELSEIF (.NOT. NEARRD(V(1,IV),GRED,ERRSIZ)) THEN
IF (DEBUG) WRITE (*,*) ' '//hm//' FAILED NEARRD'
ELSE
DBEST = DTEST
CALL CPYVN(6,TG,VBEST)
CALL CPYVN(6,AG,AVBEST)
NBEST = IV
IF (DEBUG) WRITE (*,*) ' '//hm//
2 ' NBEST,DBEST ',NBEST,DBEST
ENDIF
3300 CONTINUE
3400 CONTINUE
C output those cases that are acceptable
cutoff = amax1(10.0,amin1(errsiz,999.0))
if (nbest .ne. 0) dbest = sqrt(DOTVN(6,AVBEST,AVBEST))
IF (nbest .ne. 0 .and. DBEST .LE. cutoff) THEN
nmatch = nmatch + 1
if(ostyle.ne.'CIF ') then
WRITE (*,*)
WRITE (*,*)
write (*,'(1x,I3,A,A,A,F7.2,A,A,A,I2,A)')
2 i,' ',chrlat(i),' '//hm//' ',dbest,
3 ' = Distance Projected',' ',
4 'Internat. Tables#(',itdesg(i),')'
else
cifres = pchar_(' ',cifeid)
cifres = pnumb_(' ',float(i),0.)
cifres = pchar_(' ',cifsgs)
cifres = pchar_(' ',chrlat(i))
endif
if(ostyle.ne.'CIF') then
WRITE (*,'('' '//hm//' Vector '',6F8.1)')
2 VBEST
endif
if (G6TOC (VBEST,COUT,'G6TOC ')) then
if(ostyle.ne.'CIF ') then
WRITE (*,'(2x,'''//hm//
2 ' cell'',/12x,2(3F8.2,'' ''))') COUT
else
do ii = 1,6
if(cout(ii).ge.0.)itemp = cout(ii)*100.+.5
if(cout(ii).lt.0.)itemp = cout(ii)*100.-.5
cifres = pnumb_(' ',float(itemp)/100.,0.)
enddo
endif
call g6tor3 (matref(1,nbest),m3t)
call inver (m3t,m3ti)
call unredc (itdesg(i),' ',vbest,retcel,m3,cellst)
if(ostyle.eq.'CIF ') then
do ii = 1,6
cifres = pchar_(' ',cellst(ii))
enddo
endif
call trnspz (m3,m3t)
call wrcent(chrlat(i),retcel)
if (debug) then
write (*,'(3(1x,'' '//hm//' '',10x,3i3,/))')
* (int(m3(im3)),im3=1,9)
write (*,*)
write (*,*) ' '//hm//' matref ,nbest = ',
2 nbest,' of ',nv
write (*,'(6(1x,'' '//hm//
2 ' '',3x,3f8.3,2x,3f8.3/))')
3 (matref(im,nbest),im=1,36)
endif
call g6tor3 (matref(1,nbest),m3ptrd)
call matmul (m3t,m3ptrd,m3ti)
if(ostyle.ne.'CIF ') then
write (*,*)' '//hm//' transformation from ',
2 '3-space original cell'
if (index('PR',chrlat(i)(2:2)) .ne. 0) then
write (*,*) ' '//hm//' to final primitive cell'
else
write (*,*) ' '//hm//' to final centered cell'
endif
write (*,'(3(1x,'' '//hm//' '',8x,3f7.3/))') m3ti
else
cifres= pcmnt_(char(0))
endif
endif
ENDIF
4000 CONTINUE
ENDIF
if (nmatch .eq. 0) then
if(ostyle.ne.'CIF ') then
WRITE (*,*) ' '//hm//
2 '*****************************************'
3 //'*************'
write (*,*) ' '//hm//
2 ' NO MATCHES WERE FOUND WITH THE SPECIFIED '
3 //'UNCERTAINTIES'
WRITE (*,*) ' '//hm//
2 '*****************************************'
3 //'*************'
else
cifres = pcmnt_(char(0))
cifres = pcmnt_(
* ' NO MATCHES WERE FOUND WITH THE SPECIFIED UNCERTAINTIES')
endif
endif
IF (querst.ne.'NO') THEN
if (.not. eof) THEN
if(ostyle.ne.'CIF ') then
WRITE (*,*)
WRITE (*,*)
else
cifres = pcmnt_(char(0))
cifres = pcmnt_(' ')
endif
go to 1000
endif
ENDIF
if (istyle.eq.'CIF ') then
if(loop_) goto 1000
endif
if (ostyle.eq.'CIF ' .or. istyle.eq.'CIF ') then
call close_
endif
END
C***********************************************************************
SUBROUTINE REDUCE (VI,M,VOUT,TEST)
C for a given input vector, determine the reduced vector and the
C transformation from one to the other
include 'ITERATE.cmn'
CHARACTER *6 TEST
REAL VI(6),VIN(6), VOUT(6)
LOGICAL AGAIN
real M(36),M1(36),M2(36),mnorm(36)
REAL ZEROS(6)
DATA ZEROS /6*0.0/
C-----------------------------------------------------------------------
IF (TEST .NE. 'REDUCE') THEN
WRITE (*,*) ' '//hm//' TEST WAS WRONG IN REDUCE'
STOP
ENDIF
CALL CPYVN(6,VI,VIN)
NCYCLE = 0
1000 CONTINUE
LAST = 0
CALL MKNORM (VIN,Mnorm,VOUT,'MKNORM')
call mm6 (mnorm,m,m2)
CALL CPYVN(36,m2,m)
CALL CPYVN(6,VOUT,VIN)
CALL RUNTMN (6,M1)
IF (ABS(VIN(4)) .GT. ABS(VIN(2))) THEN
M1(14) = 1.0
M1(16) = -SIGN (1.0,VIN(4))
M1(20) = -2.0*SIGN(1.0,VIN(4))
M1(30) = M1(16)
AGAIN = .TRUE.
CALL mm6(M1,M,M2)
CALL CPYVN (36,M2,M)
CALL RMV6(VIN,M1,VOUT)
LAST = 5
ELSEIF (ABS(VIN(5)) .GT. ABS(VIN(1))) THEN
M1(13) = 1
M1(17) = -SIGN (1.0,VIN(5))
M1(24) = M1(17)
M1(25) = 2*M1(17)
AGAIN = .TRUE.
CALL mm6(M1,M,M2)
CALL CPYVN (36,M2,M)
CALL RMV6(VIN,M1,VOUT)
LAST = 6
ELSEIF (ABS(VIN(6)) .GT. ABS(VIN(1))) THEN
M1(7) = 1
M1(12) = -SIGN(1.0,VIN(6))
M1(23) = M1(12)
M1(31) = 2*M1(12)
AGAIN = .TRUE.
CALL mm6(M1,M,M2)
CALL CPYVN (36,M2,M)
CALL RMV6(VIN,M1,VOUT)
LAST = 7
ELSEIF (VIN(4)+VIN(5)+VIN(6)+ABS(VIN(1))+ABS(VIN(2)) .LT. 0.0)
2 THEN
DO 2000 I=13,18
2000 M1(I) = 1
M1(20) = 2
M1(24) = 1
M1(25) = 2
M1(30) = 1
AGAIN = .TRUE.
CALL mm6(M1,M,M2)
CALL CPYVN (36,M2,M)
CALL RMV6(VIN,M1,VOUT)
LAST = 8
ELSEIF ( (VIN(4).EQ.VIN(2) .AND. 2.0*VIN(5).LT.VIN(6)) .OR.
2 (VIN(4).EQ.-ABS(VIN(2)) .AND. VIN(6).LT. 0.0) )
3 THEN
M1(14) = 1
M1(16) = -SIGN(1.0,VIN(4))
M1(20) = 2*M1(16)
M1(30) = M1(16)
AGAIN = .TRUE.
CALL mm6(M1,M,M2)
CALL CPYVN (36,M2,M)
CALL RMV6(VIN,M1,VOUT)
LAST = 15
ELSEIF ( (VIN(5).EQ.VIN(1) .AND. 2.0*VIN(4).LT.VIN(6)) .OR.
2 (VIN(5).EQ.-ABS(VIN(1)) .AND. VIN(6).LT.0.0) )
3 THEN
M1(13) = 1
M1(17) = -SIGN(1.0,VIN(5))
M1(24) = M1(17)
M1(25) = 2*M1(17)
AGAIN = .TRUE.
CALL mm6(M1,M,M2)
CALL CPYVN (36,M2,M)
CALL RMV6(VIN,M1,VOUT)
LAST = 16
ELSEIF ( (VIN(6).EQ.VIN(1) .AND. 2.0*VIN(4).LT.VIN(5)) .OR.
2 (VIN(6).EQ.-VIN(1) .AND. VIN(5).LT.0.0) ) THEN
M1(7) = 1
M1(12) = -SIGN(1.0,VIN(6))
M1(23) = M1(12)
M1(31) = 2*M1(12)
AGAIN = .TRUE.
CALL mm6(M1,M,M2)
CALL CPYVN (36,M2,M)
CALL RMV6(VIN,M1,VOUT)
LAST = 17
ELSEIF ( (VIN(4)+VIN(5)+VIN(6)+ABS(VIN(1))+ABS(VIN(2)).EQ.0.0)
2 .AND. ( 2.0*(ABS(VIN(1))+VIN(5))+VIN(6).GT.0.0) ) THEN
DO 3000 I=13,18
3000 M1(I) = 1
M1(20) = 2
M1(24) = 1
M1(25) = 2
M1(30) = 1
AGAIN = .TRUE.
CALL mm6(M1,M,M2)
CALL CPYVN (36,M2,M)
CALL RMV6(VIN,M1,VOUT)
LAST = 18
ELSE
AGAIN = .FALSE.
CALL CPYVN (6,VIN,VOUT)
call cpyvn (36,m2,m)
ENDIF
CALL MKNORM (VOUT,Mnorm,VIN,'MKNORM')
CALL mm6 (mnorm,m,m2)
CALL CPYVN(36,M2,M)
CALL CPYVN(6,VIN,VOUT)
IF (VIN(1).LT. 0.0 .OR. VIN(2).LT.0.0 .OR. VIN(3).LT.0.0) THEN
WRITE (*,*) ' '//hm//' NEG. SQ. AXIS ',NCYCLE
CALL WRVEC6(VIN,ZEROS,'WRVEC6')
CALL WRVEC6(VOUT,ZEROS,'WRVEC6')
if (istyle.ne.'CIF ') READ (*,*)
ENDIF
NCYCLE = NCYCLE + 1
IF (NCYCLE .LT. 25 .AND. AGAIN) GO TO 1000
END
C***********************************************************************
SUBROUTINE RMV6 (V1,M,V2)
REAL V1(6),V2(6)
REAL M(36)
C-----------------------------------------------------------------------
DO 3000 I=1,6
SUM = 0.0
DO 2000 J=1,6
SUM = SUM + M(6*(I-1)+J)*V1(J)
2000 CONTINUE
V2(I) = SUM
3000 CONTINUE
END
function root (a)
if (a.ge.0.0) then
root = sqrt(a)
else
root = 0
endif
end
C***********************************************************************
SUBROUTINE RUNTMN (N,M)
INTEGER N
REAL M(N,N)
C-----------------------------------------------------------------------
DO 1000 I=1,N
DO 1000 J=1,N
1000 M(I,J) = 0.0
DO 2000 I=1,N
2000 M(I,I) = 1.0
END
function sqr (a)
sqr = a*a
end
C**********************************************************************C
FUNCTION TREELN (NVEC,A,B)
C-----GET THE SEPARATION BETWEEN THE ENDS OF TWO VECTORS
C used by bldtre, insphr and nearst
DIMENSION A(NVEC),B(NVEC)
C----------------------------------------------------------------------C
SUM = 0.0
DO 1000 I=1,NVEC
SUM = SUM + (A(I)-B(I))**2
1000 CONTINUE
TREELN = SQRT(SUM)
END
C**********************************************************************C
SUBROUTINE TRSTCK (NEXT,ISTAK,ISTKP)
C helper routine for NEARST and INSPHR
include 'ITERATE.cmn'
INTEGER ISTAK(1000)
LOGICAL DEBUG
DATA DEBUG /.FALSE./
C----------------------------------------------------------------------C
IF (DEBUG) WRITE (*,*) ' '//hm//
2 ' IN TRSTCK, STACK POINTER,IPOINT ',
3 ISTKP,NEXT
ISTKP = ISTKP + 1
ISTAK(ISTKP) = NEXT
END
C**********************************************************************C
subroutine unitmx(n,a)
C produce a unit matrix of order n
real a(n,n)
C-----------------------------------------------------------------------
do 1000 i=1,n
do 1000 j=1,n
1000 a(i,j) = 0.0
do 2000 i=1,n
2000 a(i,i) = 1.0
end
C**********************************************************************C
function unitsn (a)
C-----------------------------------------------------------------------
if (a.ge.0.0) then
unitsn = 1
else
unitsn = -1
endif
end
C**********************************************************************C
subroutine unredc (itcase,rfcase,v,cell,m,cellst)
C given a g6 vector and a particular bravais lattice (designated either
C by itcase or rfcase, that is by the International tables or by Niggli
C and Roof's designations), compute the standard (often non-reduced)
C unit cell
real m(9),v(6),cell(6)
character *3 rfcase
character *8 cellst(6),ctemp
integer itcase
real redv(6),altcel(6),v1(3),v2(3),v3(3),a(3),b(3),c(3),vtemp(3)
C-----------------------------------------------------------------------
PI = 4.0*ATAN(1.0)
MAXINT = 32768.0
do 1000 I=1,6
Cell(I) = Maxint
altcel(i) = maxint
cellst(i) = '.'
1000 continue
do 1100 I = 1,9
M(I) = 0.0
1100 continue
do 1200 I=1,6
redv(I) = root(v(I))
1200 continue
redv(4) = v(4)/(redv(2)*redv(3))
redv(5) = v(5)/(redv(1)*redv(3))
redv(6) = v(6)/(redv(1)*redv(2))
IF (itcase .eq. 3 .or. rfcase .EQ. '44A') THEN
Cell(1) = redv(1)
CALL UNITMX(3,M)
ELSEIF (itcase .eq. 5 .or. rfcase .EQ. '44B') THEN
Cell(1) = root(4.0/3.0*v(1))
CALL UNITMX(3,M)
M(7) = 1.0
M(2) = 1.0
M(6) = 1.0
ELSEIF (itcase .eq. 1 .or. rfcase .EQ. '44C') THEN
Cell(1) = root(2*v(1))
do 4400 i = 1,9
M(I) = 1.0
4400 CONTINUE
M(4) = -1.0
M(8) = -1.0
M(3) = -1.0
ELSEIF (itcase .eq. 11 .or. rfcase .EQ. '45A') THEN
Cell(1) = redv(1)
Cell(3) = redv(3)
CALL UNITMX(3,M)
ELSEIF (itcase .eq. 21 .or. rfcase .eq. '45B') THEN
Cell(1) = redv(2)
Cell(3) = redv(1)
M(4) = 1.0
M(8) = 1.0
M(3) = 1.0
ELSEIF (itcase .eq. 15 .or. rfcase .eq. '45C') THEN
Cell(1) = redv(1)
Cell(3) = root(4*v(3)-2*v(1))
M(1) = 1.0
M(5) = 1.0
M(9) = 2.0
M(3) = 1.0
M(6) = 1.0
ELSEIF (itcase .eq. 6 .or. rfcase .eq. '45D') THEN
Cell(1) = root(2*v(1)+v(6))
Cell(3) = root(2*v(1)+v(4))
do 4500 I= 1,9
M(I) = 1.0
4500 CONTINUE
do 4510 I = 1,9,4
M(I) = 0.0
4510 CONTINUE
ELSEIF (itcase .eq. 7 .or. rfcase .eq. '45d') THEN
Cell(1) = root(2*v(1)+v(6))
Cell(3) = root(2*v(1)+v(4))
do 4550 I=1,9,4
M(I) = 1.0
4550 CONTINUE
M(4) = 0.0
M(8) = 0.0
M(3) = 0.0
ELSEIF (itcase .eq. 18 .or. rfcase .eq. '45E') THEN
Cell(1) = redv(1)
Cell(3) = root( 2*v(3)-0.5*v(1))
M(4) = -1.0
M(7) = 1.0
M(2) = 1.0
M(5) = -1.0
M(8) = -1.0
M(3) = 1.0
ELSEIF (itcase .eq. 12 .or. rfcase .eq. '48A') THEN
Cell(1) = redv(1)
Cell(3) = redv(3)
CALL UNITMX(3,M)
ELSEIF (itcase .eq. 22 .or. rfcase .eq. '48B') THEN
Cell(1) = redv(3)
Cell(3) = redv(1)
M(4) = 1.0
M(8) = 1.0
M(3) = 1.0
ELSEIF (itcase .eq. 9 .or. rfcase .eq. '49B') THEN
Cell(1) = redv(3)
Cell(4) = 1-2*Sqr(redv(1)/(2*redv(3)))
altcel(1) = redv(1)
altcel(3) = root(9*v(3)-3*v(1))
M(1) = 1.0
M(2) = -1.0
M(5) = 1.0
M(3) = -1.0
M(6) = -1.0
M(9) = 3.0
ELSEIF (itcase .eq. 2 .or. itcase .eq. 4
2 .or. rfcase .eq. '49x') THEN
Cell(1) = redv(1)
Cell(4) = v(4)/(2*v(1))
altcel(1) = root(0.5*v(3)-0.25*v(4))
altcel(3) = root(9*v(1)-3*Sqr(altcel(1)))
M(1) = 1.0
M(4) = -1.0
M(2) = -1.0
M(8) = 1.0
do 4900 I= 3,9,3
M(I) = -1.0
4900 CONTINUE
ELSEIF (itcase .eq. 24 .or. rfcase .eq. '49E') THEN
Cell(1) = redv(3)
Cell(4) = v(4)/(2*v(2))
altcel(1) = root(0.5*v(3)-0.25*v(4))
altcel(3) = redv(1)
M(1) = 1.0
M(4) = 2.0
M(7) = 1.0
M(5) = -1.0
M(8) = 1.0
M(3) = 1.0
ELSEIF (itcase .eq. 32 .or. rfcase .eq. '50C') THEN
do 5000 I=1,3
Cell(I) = redv(I)
5000 CONTINUE
CALL UNITMX(3,M)
ELSEIF (itcase .eq. 36 .or. rfcase .eq. '50A') THEN
Cell(1) = redv(1)
Cell(2) = root(4*v(3)-v(1))
Cell(3) = redv(2)
M(1) = 1.0
M(2) = -1.0
M(8) = -2.0
M(6) = 1.0
ELSEIF (itcase .eq. 38 .or. rfcase .eq. '50B') THEN
Cell(1) = redv(1)
Cell(3) = redv(3)
Cell(2) = root(4*v(2)-v(1))
M(1) = 1.0
M(2) = -1.0
M(5) = -2.0
M(9) = 1.0
ELSEIF (itcase .eq. 13 .or. rfcase .eq. '50D') THEN
Cell(1) = root(2*v(1)+v(6))
Cell(2) = root(2*v(1)-v(6))
Cell(3) = redv(3)
M(1) = 1.0
M(4) = 1.0
M(2) = -1.0
M(5) = 1.0
M(9) = 1.0
ELSEIF (itcase .eq. 23 .or. rfcase .eq. '50E') THEN
Cell(1) = root(2*v(3)+v(4))
Cell(2) = root(2*v(3)-v(4))
Cell(3) = redv(1)
M(4) = 1.0
M(7) = 1.0
M(5) = -1.0
M(8) = 1.0
M(3) = 1.0
ELSEIF (itcase .eq. 40 .or. rfcase .eq. '50F') THEN
Cell(1) = redv(2)
Cell(2) = root(4*v(3)-v(2))
Cell(3) = redv(1)
M(4) = 1.0
M(5) = -1.0
M(8) = -2.0
M(3) = -1.0
ELSEIF (itcase .eq. 16 .or. rfcase .eq. '51A') THEN
Cell(1) = root(2*v(1)+v(6))
Cell(2) = root(2*v(1)-v(6))
Cell(3) = root(4*v(3)-Sqr(Cell(1)))
M(1) = 1.0
M(4) = -1.0
M(5) = 1.0
M(5) = 1.0
M(8) = 2.0
M(3) = -1.0
M(6) = -1.0
ELSEIF (itcase .eq. 26 .or. rfcase .eq. '51B') THEN
Cell(1) = redv(1)
Cell(2) = root(4*v(2)-v(1))
Cell(3) = root(4*v(3)-v(1))
M(1) = -1.0
M(4) = 2.0
M(2) = -1.0
M(8) = 2.0
M(3) = 1.0
ELSEIF (itcase .eq. 8 .or. rfcase .eq. '52A') THEN
Cell(1) = root(2*v(1)+v(6))
Cell(2) = root(2*v(1)+v(5))
Cell(3) = root(2*v(1)+v(4))
M(1) = 1.0
M(7) = 1.0
M(2) = 1.0
M(5) = 1.0
M(6) = 1.0
M(9) = 1.0
ELSEIF (itcase .eq. 19 .or. rfcase .eq. '52B') THEN
Cell(1) = redv(1)
Cell(2) = root(2*v(2)-v(4))
Cell(3) = root(2*v(3)+v(4)-v(1))
M(1) = -1.0
M(2) = -1.0
M(5) = 1.0
M(8) = 1.0
M(6) = -1.0
M(9) = 1.0
ELSEIF (itcase .eq. 42 .or. rfcase .eq. '52C') THEN
Cell(1) = redv(1)
Cell(2) = redv(2)
Cell(3) = root(4*v(3)-v(1)-v(2))
M(1) = 1.0
M(5) = 1.0
M(3) = -1.0
M(6) = -1.0
M(9) = -2.0
ELSEIF (itcase .eq. 33 .or. rfcase .eq. '53A') THEN
Cell(1) = redv(1)
Cell(2) = redv(2)
Cell(3) = redv(3)
Cell(5) = v(5)/(2*Cell(1)*Cell(3))
CALL UNITMX(3,M)
ELSEIF (itcase .eq. 35 .or. rfcase .eq. '53B') THEN
Cell(1) = redv(2)
Cell(2) = redv(1)
Cell(3) = redv(3)
Cell(5) = v(4)/(2*Cell(1)*Cell(3))
M(4) = 1.0
M(2) = 1.0
M(9) = 1.0
ELSEIF (itcase .eq. 34 .or. rfcase .eq. '53C') THEN
Cell(1) = redv(1)
Cell(2) = redv(3)
Cell(3) = redv(2)
Cell(5) = v(6)/(2*Cell(1)*Cell(3))
M(1) = 1.0
M(8) = 1.0
M(6) = 1.0
ELSEIF (itcase .eq. 39 .or. rfcase .eq. '54A') THEN
Cell(1) = root(4*v(2)-v(1))
Cell(2) = redv(1)
Cell(3) = redv(3)
Cell(5) = v(4)/(Cell(1)*Cell(3))
M(1) = 1.0
M(4) = 2.0
M(2) = 1.0
M(9) = 1.0
ELSEIF (itcase .eq. 41 .or. rfcase .eq. '54B') THEN
Cell(1) = root(4*v(3)-v(2))
Cell(2) = redv(2)
Cell(3) = redv(1)
Cell(5) = v(5)/(Cell(1)*Cell(3))
M(4) = 1.0
M(7) = 2.0
M(5) = 1.0
M(3) = 1.0
ELSEIF (itcase .eq. 37 .or. rfcase .eq. '54C') THEN
Cell(1) = root(4*v(3)-v(1))
Cell(2) = redv(1)
Cell(3) = redv(2)
Cell(5) = v(4)/(Cell(1)*Cell(3))
M(1) = 1.0
M(7) = 2.0
M(2) = 1.0
M(6) = 1.0
ELSEIF (itcase .eq. 10 .or. itcase .eq. 14
2 .or. rfcase .eq. '55A') THEN
call v2cart(v,v1,v2,v3)
call vecsum(v1,v2,A)
call vecdif(v2,v1,B)
call cpyvec(v3,c)
Cell(1) = Sqrt(Dot(A,A))
Cell(2) = Sqrt(Dot(B,B))
Cell(3) = Sqrt(Dot(C,C))
Cell(5) = -ABS(Dot(A,C)/Cell(1)/Cell(3))
IF (Cell(1) .EQ. Cell(2)) THEN
CALL UNITMX(3,M)
M(4) = 1.0
M(2) = -1.0
ELSE
M(4) = 1.0
M(7) = 1.0
M(5) = -1.0
M(8) = 1.0
M(3) = 1.0
ENDIF
ELSEIF (itcase .eq. 20 .or. itcase .eq. 25
2 .or. rfcase .eq. '55B') THEN
call v2cart(v,v1,v2,v3)
call vecsum(v3,v2,A)
call vecdif(v3,v2,B)
call cpyvec (v1,c)
Cell(1) = Sqrt(Dot(A,A))
Cell(2) = Sqrt(Dot(B,B))
Cell(3) = Sqrt(Dot(C,C))
Cell(5) = -ABS(Dot(A,C)/Cell(1)/Cell(3))
IF ((v(4).LT.0.0) .AND. (v(5).GE.0.0) .AND. (V(6).LT.0.0)) THEN
M(4) = 1.0
M(7) = 1.0
M(5) = 1.0
M(8) = -1.0
M(3) = -1.0
ELSE
M(4) = 1.0
M(7) = 1.0
M(5) = -1.0
M(8) = 1.0
M(3) = 1.0
ENDIF
ELSEIF (itcase .eq. 28 .or. rfcase .eq. '56A') THEN
call v2cart(v,A,C,v3)
do 5600 I=1,3
v3(I)=2*v3(I)
5600 continue
call vecdif(v3,A,B)
Cell(1)=Sqrt(Dot(A,A))
Cell(2)=Sqrt(Dot(B,B))
Cell(3)=Sqrt(Dot(C,C))
Cell(5)=-ABS(Dot(A,C)/Cell(1)/Cell(3))
M(1) = -1.0
M(2) = -1.0
M(8) = 2.0
M(6) = 1.0
ELSEIF (itcase .eq. 30 .or. rfcase .eq. '56B') THEN
call v2cart (v,C,A,v3)
do 5620 I=1,3
v3(I)=2*v3(I)
5620 continue
call vecdif(v3,A,B)
Cell(1)=Sqrt(Dot(A,A))
Cell(2)=Sqrt(Dot(B,B))
Cell(3)=Sqrt(Dot(C,C))
Cell(5)=-ABS(Dot(A,C)/Cell(1)/Cell(3))
M(4) = -1.0
M(5) = -1.0
M(8) = 2.0
M(3) = 1.0
ELSEIF (itcase .eq. 29 .or. rfcase .eq. '56C') THEN
call v2cart (v,A,v2,C)
do 5640 I=1,3
v2(I)=2*v2(I)
5640 continue
call vecdif(v2,A,B)
Cell(1)=Sqrt(Dot(A,A))
Cell(2)=Sqrt(Dot(B,B))
Cell(3)=Sqrt(Dot(C,C))
Cell(5)=-ABS(Dot(A,C)/Cell(1)/Cell(3))
M(1) = -1.0
M(2) = -1.0
M(5) = 2.0
M(9) = 1.0
ELSEIF (itcase .eq. 43 .or. rfcase .eq. '57A') THEN
call v2cart(v,A,C,v3)
do 5700 I=1,3
v3(I)=2*v3(I)
5700 continue
call vecsum(v3,C,v3)
call vecsum(v3,A,B)
Cell(1)=Sqrt(Dot(A,A))
Cell(2)=Sqrt(Dot(B,B))
Cell(3)=Sqrt(Dot(C,C))
Cell(5)=-ABS(Dot(A,C)/Cell(1)/Cell(3))
M(1) = 1.0
M(2) = 1.0
M(5) = 1.0
M(8) = 2.0
M(6) = 1.0
ELSEIF (itcase .eq. 17 .or. rfcase .eq. '57B') THEN
call v2cart(v,v1,v2,v3)
call vecsum(v1,v2,B)
call vecsum(v3,v2,C)
call vecsum(v3,v1,A)
Cell(1)=Sqrt(Dot(A,A))
Cell(2)=Sqrt(Dot(B,B))
Cell(3)=Sqrt(Dot(C,C))
Cell(5)=-ABS(Dot(A,C)/Cell(1)/Cell(3))
M(4) = 1.0
M(7) = 1.0
M(2) = 1.0
M(5) = 1.0
M(3) = 1.0
M(9) = 1.0
ELSEIF (itcase .eq. 27 .or. rfcase .eq. '57C') THEN
call v2cart(v,B,v2,v3)
call vecsum(v2,v3,vtemp)
call vecdif(vtemp,B,C)
call vecdif(v2,v3,A)
Cell(1)=Sqrt(Dot(A,A))
Cell(2)=Sqrt(Dot(B,B))
Cell(3)=Sqrt(Dot(C,C))
Cell(5)=-ABS(Dot(A,C)/Cell(1)/Cell(3))
M(4) = 1.0
M(7) = -1.0
M(2) = 1.0
M(3) = -1.0
M(6) = 1.0
M(9) = 1.0
ENDIF
do 8000 I = 4,6
IF (Cell(I) .NE. Maxint) THEN
Cell(I) = 180.0/PI*acos(Cell(I))
write(cellst(i),'(f8.2)') cell(i)
ENDIF
8000 CONTINUE
do 8100 i = 1,3
if(cell(i).ne.Maxint)
* write(cellst(i),'(f8.2)') cell(i)
8100 CONTINUE
do 8500 I = 4,6
IF (altcel(I) .NE. Maxint) THEN
altcel(I) = 180.0/PI*acos(altcel(I))
ENDIF
8500 CONTINUE
do 8600 I = 1,6
do 8700 II = 1,8
if (cellst(I)(II:II).ne.' ') goto 8800
8700 continue
cellst(I)='.'
return
8800 ctemp = cellst(I)(II:8)
cellst(I)=ctemp
8600 continue
return
END
C**********************************************************************C
subroutine v2Cart(v ,v1,v2,v3 )
C Compute the 3-space Cartesianizing transformation matrix corresponding
C to a particular g6 vector. The base vectors of the transformation
C are returned
real v(6),v1(3),v2(3),v3(3),cell(6)
real mat8
C-----------------------------------------------------------------------
Cell(1) = Root(v(1))
Cell(2) = Root(v(2))
Cell(3) = Root(v(3))
Cell(4) = 0.5*v(4)/(Cell(2)*Cell(3))
Cell(5) = 0.5*v(5)/(Cell(1)*Cell(3))
Cell(6) = 0.5*v(6)/(Cell(1)*Cell(2))
SinAl = Sqrt(1-Sqr(Cell(4)))
SinBe = Sqrt(1-Sqr(Cell(5)))
SinGa = Sqrt(1-Sqr(Cell(6)))
v1(1) = Cell(1)
v1(2) = 0
v1(3) = 0
v2(1) = Cell(2) * Cell(6)
v2(2) = Cell(2) * SinGa
v2(3) = 0
v3(1) = Cell(3) * Cell(5)
Mat8 = (Cell(4)-Cell(5)*Cell(6)) / SinGa
v3(2) = Mat8 * Cell(3)
v3(3) = Cell(3)*Sqrt( Sqr(SinBe)-Sqr(Mat8) )
end
C***********************************************************************
SUBROUTINE WRCELL (C,CE,TEST)
include 'ITERATE.cmn'
REAL C(6),CE(6)
CHARACTER *6 TEST
C-----------------------------------------------------------------------
IF (TEST .NE. 'WRCELL') THEN
WRITE (*,*) ' '//hm//' TEST WAS WRONG IN WRCELL'
STOP
ENDIF
SUM = 0.0
DO 1000 I=1,6
1000 SUM = SUM + CE(I)
if (ostyle.ne.'CIF ') then
IF (SUM .EQ. 0) THEN
WRITE (*,'(1x,a3,1X,6F10.3)') ' '//hm//' ',C
ELSE
DO 3000 I=1,6
IF (CE(I) .GT. 0.0) THEN
WRITE (*,'(1x,a3,1X,F10.3,3X,F10.3)')
* ' '//hm//' ',C(I),CE(I)
ELSE
WRITE (*,'(1x,a3,1X,F10.3,3X,F10.3)') ' '//hm//' ',C(I)
ENDIF
3000 CONTINUE
ENDIF
WRITE (*,*)
endif
END
C***********************************************************************
subroutine wrcent(lat,retcel)
include 'ITERATE.cmn'
character *2 lat
real retcel(6)
character *79 line
C-----------------------------------------------------------------------
if (lat(2:2) .eq. 'P') then
line = ' primitive'
elseif (lat(2:2) .eq. 'S') then
line = ' side-centered '
elseif (lat(2:2) .eq. 'I') then
line = ' body-centered '
elseif (lat(2:2) .eq. 'F') then
line = ' face-centered '
elseif (lat(2:2) .eq. 'R') then
line = ' as rhomboh. '
endif
write (line(16:),'(f10.3)') retcel(1)
if (lat(1:1) .eq. 'm') then
write (line(26:),'(f10.3)') retcel(2)
write (line(36:),'(f10.3)') retcel(3)
write (line(49:),'(a)') 'beta'
write (line(56:),'(f10.3)') retcel(5)
elseif (lat(1:1) .eq. 'o') then
write (line(26:),'(f10.3)') retcel(2)
write (line(36:),'(f10.3)') retcel(3)
elseif (lat(1:1) .eq. 't') then
write (line(26:),'(f10.3)') retcel(3)
elseif (lat .eq. 'hR') then
write (line(29:),'(a)') 'alpha'
write (line(36:),'(f10.3)') retcel(4)
elseif (lat(1:1) .eq. 'h') then
write (line(26:),'(f10.3)') retcel(3)
elseif (lat(1:1) .eq. 'c') then
endif
if (ostyle.ne.'CIF ') write (*,'(a)') ' '//hm//' '//line
end
C***********************************************************************
SUBROUTINE WRMATR (N,A)
include 'ITERATE.cmn'
REAL A(N,N)
C-----------------------------------------------------------------------
DO 1000 I=1,N
WRITE (*,*) ' '//hm,(A(J,I),J=1,N)
1000 CONTINUE
END
C***********************************************************************
SUBROUTINE WRSIZE (SIZE,ERRSIZ,RATIO,TEST)
include 'ITERATE.cmn'
CHARACTER *6 TEST
C-----------------------------------------------------------------------
IF (TEST .NE. 'WRSIZE') THEN
WRITE (*,*) ' '//hm//' TEST WAS WRONG IN WRSIZE'
STOP
ENDIF
if (ostyle.ne.'CIF ') WRITE
2 (*,'('' '//hm//
3 ' INPUT VECTOR SIZE AND ERROR AND RATIO '',3F10.2)')
4 SIZE,ERRSIZ,RATIO
END
C***********************************************************************
SUBROUTINE WRVEC6(V,VE,TEST)
include 'ITERATE.cmn'
REAL V(6),VE(6)
CHARACTER *6 TEST
C-----------------------------------------------------------------------
IF (TEST .NE. 'WRVEC6') THEN
WRITE (*,*) ' '//hm//' TEST WAS WRONG IN WRVEC6'
STOP
ENDIF
SUM = 0.0
DO 1000 I=1,6
1000 SUM = SUM + ABS(VE(I))
if (ostyle.ne.'CIF ') then
IF (SUM .EQ. 0) THEN
WRITE (*,'(1x,a3,1X,6F10.2)') ' '//hm//' ',V
ELSE
DO 3000 I=1,6
WRITE (*,'(1x,a3,1X,F10.3,3X,F10.3)') ' '//hm//' '
2 ,V(I),VE(I)
3000 CONTINUE
ENDIF
WRITE (*,*)
endif
END
C**********************************************************************C
SUBROUTINE ZEROS (N,V)
REAL V(N)
C----------------------------------------------------------------------C
DO 1000 I=1,N
1000 V(I) = 0.0
END
C**********************************************************************C
SUBROUTINE INVER (A,B)
C----INVERT A THREE BY THREE MATRIX
REAL X(9)
DIMENSION A(9),B(9)
DIMENSION IDATA1(3),IDATA2(3)
DATA IDATA1 /4,7,1/
DATA IDATA2 /7,1,4/
C----------------------------------------------------------------------C
J = 0
DO 1000 I=1,9,3
J = J + 1
ID1 = IDATA1(J)
ID2 = IDATA2(J)
CALL CROSS(A(ID1),A(ID2),X(I))
1000 CONTINUE
DETA = DET(A)
IF (ABS(DETA) .LE. 1.0E-20) THEN
DETA = SIGN(1.0E-20,DETA)
ELSE
DETA = 1.0 / DETA
ENDIF
CALL CONMAT (X,DETA,X)
CALL TRNSPZ(X,B)
END
C**********************************************************************C
SUBROUTINE VECSUM (X,Y,Z)
C----ADD TWO VECTORS AND RETURN THE SUM IN Z
DIMENSION X(3), Y(3), Z(3)
C----------------------------------------------------------------------C
DO 1000 I=1,3
Z(I) = X(I) + Y(I)
1000 CONTINUE
END
C**********************************************************************C
SUBROUTINE VECDIF (X,Y,Z)
C----SUBTRACT TWO VECTORS AND RETURN THE RESULT IN Z
DIMENSION X(3), Y(3), Z(3)
C----------------------------------------------------------------------C
DO 1000 I=1,3
Z(I) = X(I) - Y(I)
1000 CONTINUE
END
C**********************************************************************C
SUBROUTINE CPYVEC (X,Y)
C----COPY A VECTOR X INTO A VECTOR Y
DIMENSION X(3), Y(3)
C----------------------------------------------------------------------C
DO 1000 I=1,3
Y(I) = X(I)
1000 CONTINUE
END
C**********************************************************************C
FUNCTION DOT (X,Y)
C----COMPUTE AND RETURN THE DOT PRODUCT OF X AND Y
DIMENSION X(3),Y(3)
C----------------------------------------------------------------------C
DOT = 0.0
DO 1000 I=1,3
DOT = DOT + X(I) * Y(I)
1000 CONTINUE
END
C**********************************************************************C
SUBROUTINE TRNSPZ (A,B)
C----PUT THE TRANSPOSE OF A INTO B
DIMENSION A(9), B(9)
C----------------------------------------------------------------------C
J = 0
DO 1000 I=1,9,3
J = J + 1
CALL UNVEC (A(I),B(J),B(J+3),B(J+6))
1000 CONTINUE
END
C**********************************************************************C
SUBROUTINE MATMUL (A,B,C)
C----MULTIPLY TWO MATRICIES
REAL X(9)
DIMENSION A(9), B(9), C(9)
C----GET THE TRANSPOSE OF B INTO X
C----------------------------------------------------------------------C
CALL TRNSPZ (B,X)
IJ = 0
DO 2000 I=1,9,3
DO 2000 J=1,9,3
IJ = IJ + 1
C(IJ) = DOT(A(I),X(J))
2000 CONTINUE
END
C**********************************************************************C
SUBROUTINE CROSS (X,Y,Z)
C----COMPUTE Z = X CROSS Y
DIMENSION X(3),Y(3),Z(3)
C----------------------------------------------------------------------C
Z(1) = X(2)*Y(3) - Y(2)*X(3)
Z(2) =-X(1)*Y(3) + Y(1)*X(3)
Z(3) = X(1)*Y(2) - Y(1)*X(2)
END
C**********************************************************************C
FUNCTION DET(A)
C----RETURN THE VALUE OF THE DETERMINANT OF A MATRIX
REAL X(3)
DIMENSION A(9)
C----------------------------------------------------------------------C
CALL CROSS (A(1),A(4),X)
DET = DOT (X,A(7))
END
C**********************************************************************C
SUBROUTINE CONMAT(AMAT,X,BMAT)
DIMENSION AMAT(9),BMAT(9)
C----------------------------------------------------------------------C
DO 1000 I=1,9
BMAT(I) = X * AMAT(I)
1000 CONTINUE
END
C**********************************************************************C
SUBROUTINE UNVEC (X,F,G,H)
C----RETURN THE VECTOR COMPONENTS AS SCALARS
DIMENSION X(3)
C----------------------------------------------------------------------C
F = X(1)
G = X(2)
H = X(3)
END
C**********************************************************************C
SUBROUTINE WRCLV6 (C,CE,V,VE,TEST)
include 'ITERATE.cmn'
logical pnumb_
CHARACTER *6 TEST
REAL C(6),CE(6),V(6),VE(6)
C----------------------------------------------------------------------C
if (test .ne. 'WRCLV6') then
write (*,*) ' '//hm//
2 ' test string was not WRCLV6 in that routine'
stop
endif
if (ostyle.ne.'CIF ') then
DO 3000 I=1,6
WRITE (*,'(1x,a3,1X,F10.3,3X,F8.3,8X,f10.3,3x,f8.2))')
2 ' '//hm//' ',C(I),CE(I),V(I),VE(I)
3000 CONTINUE
WRITE (*,*)
else
do ii = 1,6
cifres = pnumb_(' ',c(ii),0.)
enddo
endif
END
"CUT_HERE_CUT_HERE_CUT_HERE"
#___________ THE FILE 'iterate.f' ENDS TWO LINES BEFORE THIS ___________
if ( 101777 != `wc -c <'iterate.f'` ) then
echo "cshar: 'iterate.f' unpacked with wrong size"
endif
# end of 'iterate.f'
endif
if ( -e 'iterate.html.m4' && ${1} != "-c" ) then
echo "cshar: Will not clobber existing file 'iterate.html.m4'"
sed "s/^X//" >'/dev/null' <<"CUT_HERE_CUT_HERE_CUT_HERE"
else
echo "cshar: Extracting 'iterate.html.m4' (7843 characters)"
#___________ THE FILE 'iterate.html.m4' BEGINS TWO LINES AFTER THIS __________
sed "s/^X//" >'iterate.html.m4' <<"CUT_HERE_CUT_HERE_CUT_HERE"
WWW G6 Bravais Lattice Determination
G6 Bravais Lattice
Determination Interface
by
Lawrence C. Andrews, Thuridion, Inc.,
andrews@thuridion.com and
Herbert J. Bernstein, Bernstein+Sons,
yaya@bernstein-plus-sons.com
NOTICE
Some of the software and documents included within this
software package are the intellectual property of various
parties, and placement in this package does not in anyway
imply that any such rights have in any way been waived
or diminished.
With respect to any software or documents for which a
copyright exists, ALL RIGHTS ARE RESERVED TO THE OWNERS
OF SUCH COPYRIGHT.
Even though the authors of the various documents
and software found here have made a good faith
effort to ensure that the documents are correct and
that the software performs according to its
documentation, and we would greatly appreciate
hearing of any problems you
may encounter, the programs and documents
any files created by the programs are provided
**AS IS** without any warrantee as to
correctness, merchantability or fitness
for any particular or general use.
THE RESPONSIBILITY FOR ANY ADVERSE
CONSEQUENCES FROM THE USE OF PROGRAMS
OR DOCUMENTS OR ANY FILE OR FILES CREATED BY
USE OF THE PROGRAMS OR DOCUMENTS
LIES SOLELY WITH THE USERS OF THE PROGRAMS
OR DOCUMENTS
OR FILE OR FILES AND NOT WITH AUTHORS OF
THE PROGRAMS OR DOCUMENTS.
Access to the source of ITERATE
This program and related scripts are available as
a self-extracting shell-script archive or as
a self-extracting C-shell-script archive.
What Does This Web Page Do?
In simple terms, what this page does is to find the cells which are
"close" to the cell given, in order to help find the
Bravais lattice of highest symmetry consistent with the cell.
A central problem in the solution of every crystal structure
is to determine the correct Bravais lattice of the crystal.
The Bravais lattices as they are usually listed
are:
| aP | triclinic (anorthic) primitive
|
| mP | monoclinic primitive
|
| mS | monoclinic side-centered (usually C-centered)
|
| oP | orthorhombic primitive
|
| oS | orthorhombic side-centered
|
| oF | orthorhombic face-centered
|
| oI | orthorhombic body-centered
|
| hP | hexagonal primitive
|
| hR | hexagonal rhombohedrally-centered
|
| tP | tetragonal primitive
|
| tI | tetragonal body-centered
|
| cP | cubic primitive
|
| cF | cubic face-centered
|
| cI | cubic body-centered
|
Failure to find the highest correct symmetry has several consequences, the
worst of which is that the structure may not be solved. The least of the
consequences is that Richard Marsh may
publish a paper that points out the error,
corrects it, and finds a better solution
to the structure. Many methods have been
described for finding the correct Bravais lattice. A summary of the published
methods was published in the paper that
described the G6 formalism (which is used
in the program on this web page).
"Lattices and Reduced Cells as
Points in 6-Space and Selection of Bravais
Lattice Type by Projections."
Lawrence C. Andrews and Herbert J. Bernstein, Acta
Crystallographica, A44, 1009-1018 (1988).
The program on this Web page implements a
search in G6 for the various Bravais
lattices that the user's cell may fit. For each lattice type, the best metric
match is reported. If the higher symmetry type is actually correct, then that is
likely to be the best cell from which to start further refinement. However, the
possibility exists that one of the rejected cells (which did not match as well)
was actually the correct one to use. The reason for this ambiguity is
experimental error and its propagation in the transformations of the lattices in
the program. Fortunately, the rejected cells are usually quite similar to the
accepted one.
A note on standard deviations: First, even in the best of circumstances,
standard deviations of unit cell dimensions
from 4-circle diffractometer data are
always underestimated (by at least a factor
of 2). In addition, the points chosen
for the determination are often not well distributed (for example all in the
first octant of orthorhombic lattices). These less than optimal choices cause
substantial systematic error. The experimental errors are amplified in the
mathematical conversions between various
lattices that any lattice search program
must perform. It is not a rare occurrence for angles to be incorrect by 0.5
degrees in initial unit cell determinations.
Note: Even in most well determined unit
cells, the actual errors in the edge lengths is 0.2 to 0.5 parts per thousand.
(Note that reproducibility of the measurements is substantially better, leading
to the illusion that diffractometers
produce excellent unit cell parameters). Use
of standard deviations that are too small is a common reason for failure of
Bravais lattice searches. For small molecules, 0.1 Angstroms is a reasonable
error for the edge lengths, for proteins, 0.4 to 0.5 (or even more for
preliminary measurements). Accurate unit cell parameters must by determined by a
number of more complex methods and must `include' extrapolation to remove
systematic effects. For an excellent summary,
see "Xray Structure Determination",
G.H.Stout and L.H.Jensen, Wiley, 1989.
"CUT_HERE_CUT_HERE_CUT_HERE"
#___________ THE FILE 'iterate.html.m4' ENDS TWO LINES BEFORE THIS ___________
if ( 7843 != `wc -c <'iterate.html.m4'` ) then
echo "cshar: 'iterate.html.m4' unpacked with wrong size"
endif
# end of 'iterate.html.m4'
endif
echo "cshar: End of archive 1 (of 1)."
cp /dev/null ark1isdone
set MISSING=""
foreach I ( 1 )
if ( ! -e ark${I}isdone ) then
set MISSING="${MISSING} ${I}"
endif
end
if ( "${MISSING}" == "" ) then
echo You have the archive.
rm -f ark{1,2,3,4,5,6,7,8,9}isdone
else
echo You still need to unpack the following archives:
echo " " ${MISSING}
endif
## End of C-shell archive.
exit 0