1c*********************************************************************** 2c Subroutine sym_nwc: a 'friend' of geom 3c 4c Driver routine:Attributes 5c 6c calls gensym chain: (creates matrix representaions of the 7c symmetry operators) 8c calls dosymops: (applies operatorst to symmetry distinct atoms 9c prints and stores expanded lists) 10c This routine has direct access to geom private include files and 11c thus common blocks. 12c 13c calls sym_map: creates atom pair mapping relations based upon the 14c symmetry operations and deposits symmetry operations 15c and mapping tables to rtdb for future retrieval. 16c 17c 18c call sym_init_inv_op: to make list of inverse operations 19c 20c A.C.H. 5/12/94 21c Interfaces to NWCHEM A.C.H. 10/10/94 + RJH 22c 23c*********************************************************************** 24 subroutine sym_nwc(geom,rtdb,nata,oprint,scale,threquiv,nops) 25C$Id$ 26 implicit real*8 (a-h,o-z) 27#include "errquit.fh" 28 parameter(maxops=192) 29 Parameter (EPS=1.D-14) 30 31#include "mafdecls.fh" 32#include "geom.fh" 33#include "nwc_const.fh" 34#include "geomP.fh" 35#include "rtdb.fh" 36 37 logical LResult, oprint 38 integer geom,rtdb 39 dimension cdist(3),cang(3) 40 double precision mass_new(max_cent), charge_new(max_cent) 41 double precision invnucexp_new(max_cent) 42C external pertbl 43c 44c--> parameters needed by sym chain 45c 46 itype = isystype(geom) 47 numgrp = group_number(geom) 48 numset = setting_number(geom) 49c 50c lattice vectors & angles (cell constants) 51c values should have been read in and initialized properly for 52c each system dimension 53c 54 do 1000 i=1,3 55 cdist(i)=lattice_vectors(i,geom) 56 cang(i) =lattice_angles(i,geom) 57 1000 continue 58c 59c generate symmetry operators for whatever group (crystals, surfaces 60c polymers or molecules) have been requested. 61c 62 if ((itype.eq.3).and.(numset.gt.2)) then 63 !write(*,*) "generating unsual groups" 64 call gensym_extra(itype,group_number(geom),numset, 65 $ sym_ops(1,1,geom), 66 $ nops,oprint, group_name(geom)) 67 else 68 call gensym(itype,numgrp,numset,sym_ops(1,1,geom), 69 $ nops,oprint, group_name(geom),geom,rtdb) 70 end if 71 72c 73c scratch space needed by dosymops: tags, coordinates 74c space on stack: pointer list, points to the assymetric (distinct) atoms 75c apply symmops to coordinate list 76c 77 LResult = MA_Push_Get(MT_Byte,nata*(nops+1)*16,'scratch labels', 78 & latags, iatags) 79 LResult = MA_Push_Get(MT_Dbl,nata*(nops+1)*3,'scratch atoms', 80 & lcoord, icoord) .and. lresult 81 Lresult = MA_Push_Get(MT_Int,nata,'distinct atoms',ldatom,idatom) 82 $ .and. lresult 83 if (.not. lresult) call errquit('sym_nwc: dosymops ma failed', 0, 84 & MA_ERR) 85c 86 call dosymops(sym_ops(1,1,geom),nops,coords(1,1,geom), 87 $ tags(1,geom),nata, 88 & itype,Byte_MB(iatags),Dbl_MB(icoord),nareal,nata*(nops+1), 89 $ cdist,cang,Int_MB(idatom), 90 $ charge(1,geom), charge_new, geom_mass(1,geom), 91 $ mass_new, geom_invnucexp(1,geom),invnucexp_new,threquiv) 92c 93 if(nareal.gt.max_cent) then 94 write(6,*)' Too many atoms, increase nw_max_atom in ', 95 & 'util/nwc_const.fh error in sym_nwc : ', nareal 96 call errquit('sym: too many atoms', 1, INPUT_ERR) 97 endif 98c 99c Allocate memory for the center map on the heap since it persists 100c 101 if (nops .gt. 0) then 102 if (.not. MA_Alloc_Get(MT_Int,nops*nareal,'atom sym map', 103 & lmscr,imscr)) call errquit 104 $ ('sym_nwc: ma_alloc_get of sym map failed', nops*nareal, 105 & MA_ERR) 106 else 107 lmscr = -1 108 imscr = 99999999 109 endif 110c 111c--> derive mapping table for atom pairs permuted by sym operations 112c--> NOTE: at this point coords will contain cartesian coordinates 113c and dbl_mb(icoord) has the fractional coordinates 114c (unless system is molecular in which case is also cartesian) 115c 116 if (nops .gt. 0) then 117 if (.not. MA_Push_Get(MT_Dbl,nops*3,'xnew scratch', 118 & lxnew,ixnew)) call errquit 119 $ ('sym_nwc: ma for xnew scratch', nops*3, MA_ERR) 120 else 121 ixnew = 0 122 endif 123 if (.not. MA_Push_Get(MT_INT,nareal,'ilbl', 124 & llbl,ilbl)) call errquit 125 $ ('sym_nwc: ma for ilbl', nareal, MA_ERR) 126c 127 call sym_map(itype, Byte_MB(iatags), sym_ops(1,1,geom), 128 $ nops, nareal, Int_Mb(imscr),Dbl_Mb(ixnew),Dbl_Mb(icoord), 129 $ oprint,INT_MB(ilbl), threquiv) 130c 131 LResult = MA_Pop_Stack(llbl) 132 if (nops .gt. 0) then 133 LResult = MA_Pop_Stack(lxnew) 134 endif 135c 136c Insert new information into the geometry data structures 137c Geom will hold the cartesian coords ... fractional ones 138c remain in icoord. 139c Note that the memory for the map (lmscr) is now freed in geom_destroy 140c 141 call sym_put_in_geom(geom, nata, nareal, 142 $ Byte_MB(iatags), int_mb(idatom), dbl_mb(icoord), 143 $ lmscr, imscr, nops, charge_new, mass_new, invnucexp_new) 144c 145c ... sloppy way for the moment to fix units for 2-d case 146c 147 if(itype.eq.2) then 148 do i=1,nareal 149 coords(3,i,geom)=coords(3,i,geom)*scale 150 enddo 151 endif 152c 153c---> print gmat,amat,vol evaluated by geom_nd routines 154c 155c 156c---> print gmat,amat,vol evaluated by geom_nd routines 157c 158c 159c--> Print Amatrix, G matrix, volume for 3-d and 2-d cells. 160c 161c if(itype.eq.3) then 162c write(*,9) 163c do 513 i=1,3 164c write(*,15) (amatrix(i,j,geom), j=1,3) 165c 513 continue 166c 167c write(*,21) 168c do 517 i=1,3 169c write(*,15) (amatrix_inv(i,j,geom), j=1,3) 170c 517 continue 171c 172c write(*,10) 173c do 130 i=1,3 174c write(*,15) (metric_matrix(i,j,geom), j=1,3) 175c 130 continue 176c 177c write(*,17) volume_direct(geom) 178c 179c elseif(itype.eq.2) then 180c write(*,9) 181c do 630 i=1,3 182c write(*,15) (amatrix(i,j,geom), j=1,3) 183c 630 continue 184c write(*,21) 185c do 631 i=1,3 186c write(*,15) (amatrix_inv(i,j,geom), j=1,3) 187c 631 continue 188c 189c write(*,10) 190c do 640 i=1,3 191c write(*,15) (metric_matrix(i,j,geom), j=1,3) 192c 640 continue 193c 194c write(*,18) volume_direct(geom) 195c endif 196c 197c--> Print cartesian coordinates of unit cell. 198c 199c if(itype.eq.3) then 200c write(*,30) 201c l=0 202c do 510 i=1,nareal 203c write(*,20) i,tags(i,geom),(coords(j,i,geom), j=1,3) 204c l=l+3 205c 510 continue 206c write(*,14) 207c convert z-coords of slab to bohr from angstroms 208c elseif(itype.eq.2) then 209c do i=1,nareal 210c coords(3,i,geom)=coords(3,i,geom)*angstrom_to_au 211c enddo 212c write(*,30) 213c l=0 214c do 512 i=1,nareal 215c write(*,20) i,tags(i,geom),(coords(j,i,geom), j=1,3) 216c l=l+3 217c 512 continue 218c write(*,14) 219c endif 220c 221c--> final clean up of memory: gensym and dosymops activities 222c 223 LResult = MA_Pop_Stack(ldatom) 224 LResult = MA_Pop_Stack(lcoord) 225 LResult = MA_Pop_Stack(latags) 226c 227 call sym_init_inv_op(geom) 228c 229c--> formats 230c 231 9 format(//20x, 232 & 'CRYSTALLOGRAPHIC TO CARTESIAN MATRIX (A.U.)',/) 233 10 format(//,23x,'METRIC MATRIX:DIRECT LATTICE (A.U.)',/) 234 14 format(/,80('-'),//) 235 15 format(20x,3(f12.6,2x)) 236 17 format(/,19x,'UNIT CELL VOLUME = ',f12.6,3x,'BOHR**3') 237 18 format(/,19x,'UNIT CELL AREA = ',f12.6,3x,'BOHR**3') 238 20 format(10x,i5,2x,a16,3(f10.6,2x)) 239 21 format(//20x, 240 & 'CARTESIAN TO CRYSTALLOGRAPHIC MATRIX (A.U.)',/) 241 30 format(/,15x,'--------- CARTESIAN COORDINATES (A.U.) ---------', 242 &/) 243 return 244 end 245 246 247 248