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