1 program SaveSFile 2c ********************************************************************** 3c program CONVERT XSF format to WIEN STRUCT FILE format 4c ********************************************************************** 5 IMPLICIT REAL*8 (A-H,O-Z) 6 include 'param.inc' 7 INTEGER NAT(NAC) 8 INTEGER C2I 9 REAL*8 X(NAC),Y(NAC),Z(NAC),DVEC(3,3),IDVEC(3,3),TMPVEC(3,3) 10 REAL*8 iRdvec(3,3) 11 REAL*8 A(NAC),B(NAC),C(NAC) 12 REAL*8 SYMMOP(1,3,3), SYMMTR(1,3) 13 CHARACTER*80 FILE1, TITLE, MODE 14 CHARACTER*4 GTYPE 15 LOGICAL l_inv_center, l_inv 16 COMMON/MULTAT/ X,Y,Z,NAT,DVEC,NATR 17 COMMON/IVEC/ IDVEC 18 COMMON/BOHRR/ BOHR,IMODE3 !this is not used, but is needed because 19 !this program uses readf1 routine 20 PARAMETER (RAD2DEG = 57.295779513084, TOLMIN=1.d-7) 21 22 DATA SYMMOP/ 23 $ 1.0d0, 0.0d0, 0.0d0, 24 $ 0.0d0, 1.0d0, 0.0d0, 25 $ 0.0d0, 0.0d0, 1.0d0/ 26 DATA SYMMTR/ 27 $ 0.0d0, 0.0d0, 0.0d0/ 28 29 narg=IArgc() 30 if (narg.ne.1) then 31 stop 'Usage: savestruct file1' 32 endif 33 34 Call GetArg(1,FILE1) 35 Open(unit=11,file=file1,status='old') 36 37 Call ReadF1('DIM-GROUP', igroup, idim) 38 if (idim.ne.3) stop 'SaveStruct:: so far I am converting to 39 $ WIEN97 STRUCT FILE just from crystals (dim=3) !!!' 40 41 if (igroup .eq. 1) GTYPE='P ' 42 if (igroup .eq. 2) GTYPE='CYZ ' ! A type 43 if (igroup .eq. 3) GTYPE='CXZ ' ! B type 44 if (igroup .eq. 4) GTYPE='CXY ' ! C type 45 if (igroup .eq. 5) GTYPE='F ' 46 if (igroup .eq. 6) GTYPE='B ' 47 if (igroup .eq. 7) GTYPE='R ' 48 if (igroup .eq. 8) GTYPE='H ' 49 if (igroup .eq. 9) GTYPE='H ' 50 if (igroup .gt. 10) GTYPE='P ' !just in case 51 52c this should works for H,P; but for F, B, C R, CONVVEC/PRIMCOORD 53c should probably be used !!!!!!!!!!!!!1 54 Call ReadF1('PRIMVEC', 0, idim) !dvec & idvec are assigned 55 Call ReadF1('PRIMCOORD', 0, idim) 56 57c get the inverse matrix 58 call Invert3x3(dvec,iRdvec) 59 60 if ((igroup .gt. 1) .and. (igroup .lt. 8)) then 61 Call ReadF1('CONVVEC', 0, idim) !dvec & idvec are assigned 62 if (igroup .eq. 7) then 63c for rhombohedral hexagonal vectors and rhobohedral coordinates 64c must be specified 65 do i=1,3 66 do j=1,3 67 idvec(j,i)=iRdvec(j,i) 68 enddo 69 enddo 70 endif 71 endif 72 73 Close(11) 74 75 Call SaveStructFile(gtype, dvec, idvec, natr, nat, x, y, z) 76 END 77 78 79 subroutine SaveStructFile(gtype, dvec, idvec, natr, nat, x, y, z) 80 IMPLICIT REAL*8 (A-H,O-Z) 81 include 'param.inc' 82 INTEGER NAT(NAC), INN(NAC) 83 INTEGER MAXNAT ! maximum atomic number in structure 84 INTEGER NPT, ISPLIT ! used for Wien97 85 REAL*8 X(NAC),Y(NAC),Z(NAC),DVEC(3,3),IDVEC(3,3),TMPVEC(3,3) 86 REAL*8 A(NAC),B(NAC),C(NAC), MIN_NNDIST(NAC), RMT(NAC), R0(NAC) 87 REAL*8 SYMMOP(1,3,3), SYMMTR(1,3) 88 REAL*8 V1(3), V2(3), V3(3) 89 CHARACTER*10 ATOM(100), ATOM_UPPER(100) 90 CHARACTER*4 gtype 91 REAL*8 FRMT(NAC) 92 93 PARAMETER (RAD2DEG = 57.295779513084, TOLMIN=1.d-7) 94 95 include 'atoms.inc' 96 97C CONVERSION FROM BOHRS TO ANGS. 98 BOHR=0.529177d0 99 100c **************************************************************** 101c now get the fractional coordinates of atoms and 102c translate them if necessary ( in WIEN97 fractional coorinates 103c are in range [0,1) ) 104 maxnat = 1 105 do i=1,natr 106 nat(i) = nat(i) - int( nat(i) / 100 ) * 100 ! nat from 0 to 99 107 if(nat(i).gt.maxnat) maxnat = nat(i) 108 a(i) = x(i)*idvec(1,1) + y(i)*idvec(2,1) + z(i)*idvec(3,1) 109 b(i) = x(i)*idvec(1,2) + y(i)*idvec(2,2) + z(i)*idvec(3,2) 110 c(i) = x(i)*idvec(1,3) + y(i)*idvec(2,3) + z(i)*idvec(3,3) 111 if(a(i).lt.0.0d0) a(i) = a(i) + 1.0d0 112 if(b(i).lt.0.0d0) b(i) = b(i) + 1.0d0 113 if(c(i).lt.0.0d0) c(i) = c(i) + 1.0d0 114c Write(6,21) nat(i), a(i), b(i), c(i) 115 enddo 116 117c get the nearest neighbor distances 118 do i=1,natr 119 min_nndist(i) = Dist(dvec(1,1), dvec(1,2), dvec(1,3), 120 $ dvec(2,1), dvec(2,2), dvec(2,3)) !this always greater then nndist 121 do j=1,natr 122 if (i.ne.j) then 123 a2 = a(j) + nint(a(i)-a(j)) 124 b2 = b(j) + nint(b(i)-b(j)) 125 c2 = c(j) + nint(c(i)-c(j)) 126 127 x1 = a(i)*dvec(1,1) + b(i)*dvec(2,1) + c(i)*dvec(3,1) 128 y1 = a(i)*dvec(1,2) + b(i)*dvec(2,2) + c(i)*dvec(3,2) 129 z1 = a(i)*dvec(1,3) + b(i)*dvec(2,3) + c(i)*dvec(3,3) 130 131 x2 = a2*dvec(1,1) + b2*dvec(2,1) + c2*dvec(3,1) 132 y2 = a2*dvec(1,2) + b2*dvec(2,2) + c2*dvec(3,2) 133 z2 = a2*dvec(1,3) + b2*dvec(2,3) + c2*dvec(3,3) 134 135 d = Dist(x1,y1,z1, x2,y2,z2) 136 if (d.lt.min_nndist(i)) then 137 min_nndist(i) = d 138 inn(i) = j 139 endif 140 endif 141 enddo 142 enddo 143 144c ************************* 145c get a,b,c,alpha,beta,gamma 146 ap = VecSize3(dvec,1) / BOHR 147 bp = VecSize3(dvec,2) / BOHR 148 cp = VecSize3(dvec,3) / BOHR 149 do i=1,3 150 v1(i) = dvec(1,i) 151 v2(i) = dvec(2,i) 152 v3(i) = dvec(3,i) 153 enddo 154 alpha = dacos( ScalarProduct(v2, v3) / 155 $ (VecSize(v2) * VecSize(v3)) ) * RAD2DEG 156 beta = dacos( ScalarProduct(v3,v1) / 157 $ (VecSize(v3) * VecSize(v1)) ) * RAD2DEG 158 gamma = dacos( ScalarProduct(v1,v2) / 159 $ (VecSize(v1) * VecSize(v2)) ) * RAD2DEG 160 161c ******************************************* 162c get some aproximation for muffin tin radius 163c ******************************************* 164c let the smallest offset between RMT spheres be 0.1 a.u. 165 OFFSET = 0.1d0/2.d0 ! each sphre gets half of offset 166 167 do i=1,natr 168 if (nat(i).eq.1) frmt(i) = 0.50d0 169 if (nat(i).gt.1 .and. nat(i).lt.21) frmt(i) = 1.00d0 170 if (nat(i).ge.21 .and. nat(i).lt.39) frmt(i) = 1.10d0 171 if (nat(i).ge.39 .and. nat(i).lt.57) frmt(i) = 1.15d0 172 if (nat(i).ge.57) frmt(i) = 1.20d0 173 174 r0(i) = 0.000005d0 175 if (nat(i).le.71) r0(i) = 0.00001d0 176 if (nat(i).le.36) r0(i) = 0.00005d0 177 if (nat(i).le.18) r0(i) = 0.0001d0 178 enddo 179 180 do i=1,natr 181 rmt(i) = frmt(i) * min_nndist(i) / 182 $ ( BOHR * (frmt(i) + frmt(inn(i))) ) - OFFSET 183 enddo 184 isplit = 8 185 npt = 781 186 187c **************************** 188c now write WIEN97 struct file 189c **************************** 190 write(*,'(a)') 191 $ 'Wien97 struct file generated by XCrySDen program' 192 write(*,'(a4,a23,i3)') gtype, 'LATTICE.NONEQUIV. ATOMS',natr 193 if (maxnat.lt.40) then 194 write(*,'(a13,a4)') 'MODE OF CALC=', 'NREL' 195 else 196 write(*,'(a13,a4)') 'MODE OF CALC=', 'RELA' 197 endif 198 write(*,'(6f10.6)') ap, bp, cp, alpha, beta, gamma 199 200 do i=1,natr 201c i should be positive for cubic symmetry and negative for non-cubic 202c *** ** * CHANGE THAT * ** *** 203 write(*,'(a4,i4,a4,f10.8,a3,f10.8,a3,f10.8)') 204 $ 'ATOM', -i,': X=', a(i), ' Y=', b(i), ' Z=', c(i) 205 write(*,'(a15,i2,a17,i2)') 'MULT=', 1, 'ISPLIT=', isplit 206 write(*,'(a10,a5,i5,a5,f10.8,a5,f10.4,a5,f5.1)') 207 $ atom(nat(i)),'NPT=',npt,'R0=',r0(i), 208 $ 'RMT=',rmt(i),'Z:', real(nat(i)) 209 write(*,'(a20,3f10.7)') 'LOCAL ROT MATRIX: ', 1.0,0.0,0.0 210 write(*,'(20x,3f10.7)') 0.0, 1.0, 0.0 211 write(*,'(20x,3f10.7)') 0.0, 0.0, 1.0 212 enddo 213 write(*,'(i4,A)') 1,' NUMBER OF SYMMETRY OPERATIONS' 214 write(*,'(3i2,f10.7)') 1, 0, 0, 0.0 215 write(*,'(3i2,f10.7)') 0, 1, 0, 0.0 216 write(*,'(3i2,f10.7)') 0, 0, 1, 0.0 217 write(*,'(i8)') 1 218 21 FORMAT(I5,3F16.10) 219 220 END 221 222 223 224