1c***********************************************************************
2c
3c
4c   VERY IMPORTANT: Note that the same data structures for amatrix,
5c                   amatrixinv, gmatrix, etc. are used for the 2
6c                   and 3-dimensional cases. The last row and column
7c                   of these matrices have been set equal to zero EXCEPT
8c                   element (3,3) which is equal to 1. This allows a large
9c                   amount of code reuse for the various cases. However,
10c                   Beware of doing mindless matrix algebra with
11c                   these quantities since only the (1,1),(1,2),(2,1) and
12c                   (2,2) elements have any physical meaning.
13c
14c
15c         Direct space geometry stuff for the 2-d Bravias lattices
16c
17c 1) Computes the a-matrix used to transform from crystallographic to
18c    cartesian coordinates
19c    note the following defn. of the amatrix to reduce confusion
20c
21c      A = a11 i + a21 j
22c      B = 0   i + a22 j
23c
24c      thus: Amatrix =[ a11 0   ]
25c                     [ a21 a22 ]
26c
27c     where: {A,B} are crystallographic axis, {i,j} unit Cartesian
28c            vectors. There is no concept of Crystallographic C-axis
29c
30c     So Amatrix premultiplies a fractional vector (a vector defined in
31c     the crystal basis). Each successive row of a (when mult. times
32c     the vector) will give the particular cartesian component.
33c     Atranspose*A=metrix matrix. An alternative way's to think about
34c     its definition are:
35c
36c**    Premultiplying a fractional vector by amatrix will result in a
37c**    Cartesian vector in atomic units.
38c
39c    It represents a chosen convention of two axis systems.Where the
40c    oblique crystal system is oriented relative to an external
41c    rectilinear laboratory coordinate system such that: the b-axis is
42c    along j (unit rectilinear vector) ab-plane contains i (rect)
43c    ,i.e., j along b, i in  ab-plane. This convention will result in
44c    the z-direction of Cartesian space being perpendicular to the
45c    periodic repeat direction of the 2-d slab.
46c
47c
48c 2) computes the metric matrix for the lattice.
49c
50c 3) computes the direct space cell volume
51c
52c 4) computes reciprocal space lattice constants & volume
53c
54c 5) computes the transformation matrix from reciprocal lattice vectors
55c    to Cartesians (i.e, finds cartesian components of a recip lat vec)
56c
57c 6) computes some constants needed for Ewald routines that depend on
58c    the geometry
59c
60c --> Important Conventions
61c
62c     The lattice angles are input in Degrees
63c
64c                                                   AC Hess
65c***********************************************************************
66      subroutine geom_2d(geom,c2au)
67*
68* $Id$
69*
70      implicit none
71#include "inp.fh"
72#include "geom.fh"
73#include "nwc_const.fh"
74#include "geomP.fh"
75c
76      integer geom,i,j
77      double precision c2au
78      double precision cdist(3),cang(3),gmat(3,3),amat(3,3),vol
79      double precision c3,s3,rad,ainv(3,3)
80      double precision pi,bmat(3,3)
81c
82      parameter(rad=57.295779513082343d0)
83c
84      pi=acos(0.0d0)*2.0d0
85c
86c--> get direct space lattice vectors
87c
88      do 100 i=1,3
89        cdist(i)=lattice_vectors(i,geom)*c2au
90  100 continue
91        cang(3) =lattice_angles(3,geom)/rad
92c
93c        write(*,*) 'lattice vectors and angles'
94c        write(*,*) 'vectors',(lattice_vectors(i,geom),i=1,3)
95c        write(*,*) 'angles',(lattice_angles(i,geom),i=1,3)
96c
97c--------> build the metrical matrix (atomic units)
98c
99      do 200 i=1,2
100        gmat(i,i)=cdist(i)**2
101  200 continue
102      gmat(1,2)=cdist(1)*cdist(2)*cos(cang(3))
103      gmat(2,1)=gmat(1,2)
104      gmat(1,3)=0.0d+00
105      gmat(2,3)=0.0d+00
106      gmat(3,1)=0.0d+00
107      gmat(3,2)=0.0d+00
108      gmat(3,3)=0.0d+00
109c
110      do 230 i=1,3
111        do 240 j=1,3
112          metric_matrix(i,j,geom)=gmat(i,j)
113  240   continue
114  230 continue
115c
116c      write(*,*) 'metrix matrix'
117c      do i=1,3
118c         write(*,*) (gmat(i,j),j=1,3)
119c      enddo
120c
121c-----------> build a-matrix
122c
123c***j along b, i in  ab-plane
124c
125      s3=dsin(cang(3))
126      c3=dcos(cang(3))
127      amat(1,1) = cdist(1)*s3
128      amat(1,2) = 0.0d+00
129      amat(1,3) = 0.0d+00
130      amat(2,1) = cdist(1)*c3
131      amat(2,2) = cdist(2)
132      amat(2,3) = 0.0d+00
133      amat(3,1) = 0.0d+00
134      amat(3,2) = 0.0d+00
135      amat(3,3) = cdist(3)
136c
137c      write(*,*) 'a- matrix'
138c      do i=1,3
139c         write(*,*) (amat(i,j),j=1,3)
140c      enddo
141c
142c     load it into common
143c
144      do 250 i=1,3
145        do 260 j=1,3
146          amatrix(i,j,geom)=amat(i,j)
147  260   continue
148  250 continue
149c
150c compute direct space volume in atomic bohr**3
151c
152      vol = abs(amat(1,1)*amat(2,2)-amat(1,2)*amat(2,1))
153c
154      volume_direct(geom)=vol
155c
156c--> build amatrix_inv for 2-d systems
157c
158      ainv(1,1)= amat(2,2)/vol
159      ainv(1,2)=-amat(1,2)/vol
160      ainv(1,3)= 0.0d0
161      ainv(2,1)=-amat(2,1)/vol
162      ainv(2,2)= amat(1,1)/vol
163      ainv(2,3)= 0.0d0
164      ainv(3,1)= 0.0d0
165      ainv(3,2)= 0.0d0
166      ainv(3,3)= 1.0d0/cdist(3)
167c
168c      write(*,*) 'ainv- matrix'
169c      do i=1,3
170c         write(*,*) (ainv(i,j),j=1,3)
171c      enddo
172c
173c     load it into common
174c
175      do 270 i=1,3
176        do 280 j=1,3
177          amatrix_inv(i,j,geom)=ainv(i,j)
178  280   continue
179  270 continue
180c
181c--> construct bmatrix (used to transform reciprocal vectors
182c    to Cartesian form). 2pi*(At)-1
183c
184      bmat(1,1)=(2.0d0*pi/vol)*amat(2,2)
185      bmat(1,2)=(2.0d0*pi/vol)*(-amat(2,1))
186      bmat(1,3)=0.0d0
187      bmat(2,1)=(2.0d0*pi/vol)*(-amat(1,2))
188      bmat(2,2)=(2.0d0*pi/vol)*amat(1,1)
189      bmat(2,3)=0.0d0
190      bmat(3,1)=0.0d0
191      bmat(3,2)=0.0d0
192c      bmat(3,3)=1.0d0
193      bmat(3,3)=2.0d0*pi/cdist(3)
194c
195c      write(*,*) 'b-matrix'
196c      do i=1,3
197c         write(*,*) (bmat(i,j),j=1,3)
198c      enddo
199c
200c load in into common
201c
202      do 290 i=1,3
203         do 300 j=1,3
204            bmatrix(i,j,geom)=bmat(i,j)
205 300     continue
206 290  continue
207c
208      return
209      end
210
211
212      subroutine geom_2d_amatrix(geom,scale)
213      implicit none
214      integer geom
215      double precision scale
216
217#include "errquit.fh"
218#include "nwc_const.fh"
219#include "geomP.fh"
220
221*     !**** local variables ****
222      integer i,j
223      double precision amat(3,3)
224      double precision c(3,3), vol
225c
226c
227      do i=1,3
228        do j=1,3
229          amat(i,j) = amatrix(i,j,geom)
230        end do
231      end do
232c
233c     Mmmm ... the original code only set this stuff from the input
234c     using the a,b,c,alpha,beta,gamma, but now we have changed
235c     the amatrix ... need to update ainv and also recompute the
236c     other crap ... for now just set the other crap to crap so that
237c     we'll know if it is used
238c
239      do i = 1,3
240         do j = 1,3
241            metric_matrix(i,j,geom) = 1d300
242            bmatrix(i,j,geom) = 1d300
243         end do
244         recip_lat_vectors(i,geom) = 1d300
245         recip_lat_angles(i,geom) = 1d300
246      end do
247c
248c     HERE SHOULD RECOMPUTE AMATRIX WITH STANDARD ORIENTATION
249c     SINCE IF THE GEOMETRY IS STORED AND RELOADED THE
250c     STANDARD ORIENTATION IS IMPOSED.
251c
252c     Update the amatrix inverse
253c      - Since amat=[a1,a2,a3]
254c              ainv=[b1,b2,b3]^t
255c
256      call dfill(9,0.0d0,c,1)
257      c(1,1) = amat(2,2)*amat(3,3) - amat(3,2)*amat(2,3)  ! = b(1,1)
258      c(1,2) = amat(3,2)*amat(1,3) - amat(1,2)*amat(3,3)  ! = b(2,1)
259      c(1,3) = amat(1,2)*amat(2,3) - amat(2,2)*amat(1,3)  ! = b(3,1)
260      c(2,1) = amat(2,3)*amat(3,1) - amat(3,3)*amat(2,1)  ! = b(1,2)
261      c(2,2) = amat(3,3)*amat(1,1) - amat(1,3)*amat(3,1)  ! = b(2,2)
262      c(2,3) = amat(1,3)*amat(2,1) - amat(2,3)*amat(1,1)  ! = b(3,2)
263      c(3,1) = amat(2,1)*amat(3,2) - amat(3,1)*amat(2,2)  ! = b(1,3)
264      c(3,2) = amat(3,1)*amat(1,2) - amat(1,1)*amat(3,2)  ! = b(2,3)
265      c(3,3) = amat(1,1)*amat(2,2) - amat(2,1)*amat(1,2)  ! = b(3,3)
266      vol = amat(1,1)*c(1,1)
267     >    + amat(2,1)*c(1,2)
268     >    + amat(3,1)*c(1,3)
269      volume_direct(geom) = vol
270c
271      call dscal(9,1.0d0/vol,c,1)
272c
273      call dcopy(9,c,1,amatrix_inv(1,1,geom),1)
274c
275c     Ooops ... must also update the pesky lattice parameters
276c
277      call xlattice_abc_abg(
278     $     lattice_vectors(1,geom),
279     $     lattice_vectors(2,geom),
280     $     lattice_vectors(3,geom),
281     $     lattice_angles(1,geom),
282     $     lattice_angles(2,geom),
283     $     lattice_angles(3,geom),amat)
284
285      lattice_vectors(1,geom) = lattice_vectors(1,geom)/scale
286      lattice_vectors(2,geom) = lattice_vectors(2,geom)/scale
287      lattice_vectors(3,geom) = lattice_vectors(3,geom)/scale
288c
289      return
290      end
291
292
293
294
295
296
297
298
299