1      subroutine smd_lat_init_system()
2      implicit none
3#include "errquit.fh"
4#include "inp.fh"
5#include "mafdecls.fh"
6#include "rtdb.fh"
7#include "util.fh"
8#include "global.fh"
9c
10      character*32 sp_lat
11      character*32 tag,pname
12      logical result
13
14      pname = "smd_lat_init_system"
15c
16      tag = "lattice"
17      call smd_system_get_component(sp_lat,tag,result)
18      if(.not.result)
19     >  call errquit(
20     >       pname//'no component '//tag,0,0)
21
22      call smd_lat_init(sp_lat,result)
23      if(.not.result) then
24       tag = "lattice"
25       call smd_system_unset_component(tag)
26      end if
27c
28
29      return
30      end
31
32      subroutine smd_lat_init(namespace,result)
33      implicit none
34#include "errquit.fh"
35#include "inp.fh"
36#include "mafdecls.fh"
37#include "rtdb.fh"
38#include "util.fh"
39#include "global.fh"
40c
41      character*(*) namespace
42      logical result
43c
44      integer rtdb
45      character*32 pname
46      character*80 tag
47      integer i_lc,i_lrc,i_lfc
48      double precision vol
49c
50      pname = "smd_lat_init"
51c
52c      write(*,*) "in "//pname
53c
54      call smd_rtdb_get_handle(rtdb)
55c
56c     check if there is any lattice in rtdb
57      call smd_lat_rtdb_check(rtdb,result)
58      if(.not.result) then
59        call util_warning(
60     >       pname//'no lattice found in rtdb',0,0)
61        return
62      end if
63c
64      call smd_namespace_create(namespace)
65c
66c     create lattice data structures
67c     ------------------------------
68      call smd_data_create(namespace,"lat:fconst",2,MT_DBL)
69      call smd_data_create(namespace,"lat:cell",9,MT_DBL)
70      call smd_data_create(namespace,"lat:rcell",9,MT_DBL)
71c
72c     get memory pointers
73c     -------------------
74      tag = "lat:cell"
75      call smd_data_get_index(namespace,tag,i_lc,result)
76      if(.not. result)
77     >  call errquit(
78     >       pname//'error getting index for '//tag,0, RTDB_ERR)
79      tag = "lat:rcell"
80      call smd_data_get_index(namespace,tag,i_lrc,result)
81      if(.not. result)
82     >  call errquit(
83     >       pname//'error getting index for '//tag,0, RTDB_ERR)
84      tag = "lat:fconst"
85      call smd_data_get_index(namespace,tag,i_lfc,result)
86      if(.not. result)
87     >  call errquit(
88     >       pname//'error getting index for '//tag,0, RTDB_ERR)
89
90      call smd_lat_rtdb_read(rtdb,dbl_mb(i_lc))
91      call smd_lat_invrt(dbl_mb(i_lc),dbl_mb(i_lrc))
92      call smd_latt_vol(dbl_mb(i_lc),vol)
93      dbl_mb(i_lfc) = vol
94      return
95      end
96
97      subroutine smd_lat_rtdb_check(rtdb,olatt)
98      implicit none
99#include "errquit.fh"
100#include "inp.fh"
101#include "mafdecls.fh"
102#include "rtdb.fh"
103#include "util.fh"
104#include "global.fh"
105c
106      integer rtdb
107      logical olatt
108c
109      double precision latt(3,3)
110      character*32 pname
111      character*80 tag
112      double precision a(3)
113      integer i
114c
115      pname = "smd_lat_rtdb_read"
116c
117c      write(*,*) "in "//pname
118c
119      olatt = .true.
120      tag="smd:lat_a"
121      if (.not.rtdb_get(rtdb,tag,mt_dbl,3,a(1)))
122     >      olatt=.false.
123
124      return
125      end
126
127      subroutine smd_lat_rtdb_read(rtdb,latt)
128      implicit none
129#include "errquit.fh"
130#include "inp.fh"
131#include "mafdecls.fh"
132#include "rtdb.fh"
133#include "util.fh"
134#include "global.fh"
135c
136      double precision latt(3,3)
137      integer rtdb
138c
139      character*32 pname
140      character*80 tag
141      double precision a(3)
142      integer i
143c
144      pname = "smd_lat_rtdb_read"
145c
146c      write(*,*) "in "//pname
147c
148      tag="smd:lat_a"
149      if (.not.rtdb_get(rtdb,tag,mt_dbl,3,a(1)))
150     >      call errquit(pname//'failed to get'//tag,0,
151     >       RTDB_ERR)
152      do i=1,3
153       latt(i,1)=a(i)
154      end do
155      tag="smd:lat_b"
156      if (.not.rtdb_get(rtdb,tag,mt_dbl,3,a(1)))
157     >      call errquit(pname//'failed to get'//tag,0,
158     >       RTDB_ERR)
159      do i=1,3
160       latt(i,2)=a(i)
161      end do
162      tag="smd:lat_c"
163      if (.not.rtdb_get(rtdb,tag,mt_dbl,3,a(1)))
164     >      call errquit(pname//'failed to get'//tag,0,
165     >       RTDB_ERR)
166      do i=1,3
167       latt(i,3)=a(i)
168      end do
169      return
170      end
171
172      subroutine smd_lat_invrt(latt,rlatt)
173      implicit none
174      double precision  latt(3,3),rlatt(3,3)
175c
176      double precision  det
177
178      rlatt(1,1)=latt(2,2)*latt(3,3)-latt(3,2)*latt(2,3)
179      rlatt(2,1)=latt(3,1)*latt(2,3)-latt(2,1)*latt(3,3)
180      rlatt(3,1)=latt(2,1)*latt(3,2)-latt(3,1)*latt(2,2)
181      rlatt(1,2)=latt(3,2)*latt(1,3)-latt(1,2)*latt(3,3)
182      rlatt(2,2)=latt(1,1)*latt(3,3)-latt(3,1)*latt(1,3)
183      rlatt(3,2)=latt(3,1)*latt(1,2)-latt(1,1)*latt(3,2)
184      rlatt(1,3)=latt(1,2)*latt(2,3)-latt(2,2)*latt(1,3)
185      rlatt(2,3)=latt(2,1)*latt(1,3)-latt(1,1)*latt(2,3)
186      rlatt(3,3)=latt(1,1)*latt(2,2)-latt(2,1)*latt(1,2)
187
188      det=latt(1,1)*rlatt(1,1)+latt(1,2)*rlatt(2,1)+latt(1,3)*rlatt(3,1)
189      if(abs(det).gt.0.d0)det=1.d0/det
190
191      rlatt(1,1)=det*rlatt(1,1)
192      rlatt(2,1)=det*rlatt(2,1)
193      rlatt(3,1)=det*rlatt(3,1)
194      rlatt(1,2)=det*rlatt(1,2)
195      rlatt(2,2)=det*rlatt(2,2)
196      rlatt(3,2)=det*rlatt(3,2)
197      rlatt(1,3)=det*rlatt(1,3)
198      rlatt(2,3)=det*rlatt(2,3)
199      rlatt(3,3)=det*rlatt(3,3)
200
201      return
202
203      end
204
205      subroutine smd_latt_vol(latt,vol)
206      implicit none
207      real*8 x,y,z,latt,vol
208
209      dimension latt(3,3)
210
211      x=latt(2,2)*latt(3,3)-latt(2,3)*latt(2,3)
212      y=latt(3,2)*latt(1,3)-latt(1,2)*latt(3,3)
213      z=latt(1,2)*latt(2,3)-latt(2,2)*latt(1,3)
214
215      vol=abs(latt(1,1)*x+latt(2,1)*y+latt(3,1)*z)
216
217      return
218
219      END
220
221      subroutine smd_latt_get_vol(namespace,vol)
222      implicit none
223#include "errquit.fh"
224#include "inp.fh"
225#include "mafdecls.fh"
226#include "rtdb.fh"
227#include "util.fh"
228#include "global.fh"
229c
230      character*(*) namespace
231      double precision vol
232c
233      character*72 tag
234      character*30 pname
235      integer i_fconst
236      logical result
237
238      pname = "smd_latt_vol"
239      tag = "lat:fconst"
240      call smd_data_get_index(namespace,tag,i_fconst,result)
241      if(.not. result)
242     >  call errquit(
243     >       pname//'error getting ntot '//tag,0, RTDB_ERR)
244      vol = dbl_mb(i_fconst)
245
246      return
247      end
248
249      subroutine smd_lat_rebox(n,c)
250      implicit none
251#include "errquit.fh"
252#include "inp.fh"
253#include "mafdecls.fh"
254#include "rtdb.fh"
255#include "util.fh"
256#include "global.fh"
257c
258      integer n
259      double precision c(n,3)
260c
261      character*32 sp_lattice
262c
263      character*72 tag
264      character*30 pname
265      integer na
266      integer i_c,i_lrc,i_lc
267      logical result
268
269      pname = "smd_lat_rebox"
270c
271c     get lattice params if any
272c     -------------------------
273      call smd_system_get_component(sp_lattice,"lattice",result)
274      if(.not.result) then
275        call util_warning(
276     >       pname//'skipping reboxing as there is no lattice ',0,0)
277        return
278      end if
279
280      tag = "lat:cell"
281      call smd_data_get_index(sp_lattice,tag,i_lc,result)
282      if(.not. result)
283     >  call errquit(
284     >       pname//'error getting index for '//tag,0, RTDB_ERR)
285
286      tag = "lat:rcell"
287      call smd_data_get_index(sp_lattice,tag,i_lrc,result)
288      if(.not. result)
289     >  call errquit(
290     >       pname//'error getting index for '//tag,0, RTDB_ERR)
291
292      call smd_util_rebox(n,
293     >                    dbl_mb(i_lc),
294     >                    dbl_mb(i_lrc),
295     >                    c)
296
297      return
298      end
299c $Id$
300