1! 2! Copyright (C) 1996-2016 The SIESTA group 3! This file is distributed under the terms of the 4! GNU General Public License: see COPYING in the top directory 5! or http://www.gnu.org/copyleft/gpl.txt. 6! See Docs/Contributors.txt for a list of contributors. 7! 8 9 program fcbuild 10 11c ********************************************************************* 12c Build supercell coordinates for Force Constant Matrix calculation, 13c for clusters, linear chains, slabs and 3D xtals. 14c Compatible with vibra.f 15c 16c Uses the FDF (Flexible Data Format) package (version 0.66.1.5) 17c of J.M.Soler and A. Garcia, 18c for input and output 19c 20c Produces an output file in FDF format, to be used by SIESTA 21c 22c Written by P.Ordejon, August'98 23c 24c ********************************************************************** 25 26 use fdf 27 implicit none 28 29 include 'vibra.h' 30 31c Internal variables ... 32 33 logical overflow 34 35 character 36 . filein*20, fileout*20 37 38 character 39 . slabel*20, sname*150, slabel_defect*150, sname_defect*20 40 41 42 integer 43 . i, i1, i2, iatom, imass(maxa), imasssc(maxasc), iunit, ix, j, 44 . lx, ly, lz, lxmax, lymax, lzmax, 45 . lx_defect, ly_defect, lz_defect, 46 . na_defect, natoms, ncells, nnat 47 48 real*8 49 . dx, alat, alp, blp, clp, alplp, betlp, gamlp, pi, xxx 50 51 real*8 52 . b(3,maxa), cell(3,3), r(3), scell(3,3), xa(3,maxasc), 53 . xmass(maxa) 54 55 data pi / 3.1415926d0 / 56 data overflow /.false./ 57 58 logical :: has_constr 59 60 type(block_fdf) :: bfdf 61 type(parsed_line), pointer :: pline 62 63c ... 64 65 66C ****************** READ DATA FROM FDF FILE ********************* 67 68c Set up fdf ... 69 filein = 'stdin' 70 fileout = 'out.fdf' 71 call fdf_init(filein,fileout) 72c ... 73 74c Defile Name of the system ... 75 sname_defect = ' ' 76 sname = fdf_string('SystemName',sname_defect) 77 write(6,'(a,a)') 78 . 'redata: System Name = ',sname 79c ... 80 81c Defile System Label (short name to label files) ... 82 slabel_defect = 'vibra' 83 slabel = fdf_string('SystemLabel',slabel_defect) 84 write(6,'(a,a)') 85 . 'redata: System Label = ',slabel 86c ... 87 88c Read Number of Atoms in Unit cell ... 89 na_defect = 0 90 natoms = fdf_integer('NumberOfAtoms',na_defect) 91 if (natoms .le. 0) then 92 write(6,'(a)') 93 . 'ERROR: Number of atoms must be larger than zero.' 94 write(6,'(a)') 95 . ' You MUST specify NumberOfatoms in fdf file.' 96 stop 97 endif 98 write(6,'(a,i5)') 99 . 'Number of Atoms = ',natoms 100c check if dimension of number of atomos is large enough 101 call chkdim('fcbuild', 'maxa', maxa, natoms, 1) 102c ... 103 104c Lattice constant of unit cell... 105 alat = fdf_physical('LatticeConstant',0.d0,'Bohr') 106 if (alat .eq. 0.d0) then 107 write(6,'(a)') 108 . 'ERROR: No valid lattice constant specified.' 109 write(6,'(a)') 110 . ' You MUST specify LatticeConstant in fdf file.' 111 stop 112 endif 113 write(6,'(a,f10.5,a)') 'Lattice Constant = ',alat,' Bohr' 114c ... 115 116c Whether the constrainst block exists 117 has_constr = fdf_block('GeometryConstraints',bfdf) 118 has_constr = has_constr .or. 119 & fdf_block('Geometry.Constraints',bfdf) 120c ... 121 122c Lattice vectors of unit cell... 123 if ( fdf_block('LatticeParameters',bfdf) .and. 124 . fdf_block('LatticeVectors',bfdf) ) then 125 write(6,'(2a)')'ERROR: Lattice vectors doubly ', 126 . 'specified: by LatticeVectors and by LatticeParameters.' 127 stop 128 endif 129 130 if ( fdf_block('LatticeParameters',bfdf) ) then 131 if (.not. fdf_bline(bfdf, pline)) 132 $ call die("No LatticeParameters") 133 if (.not. (fdf_bmatch(pline, 'vvvvvv') )) then 134 call die ('LatticeParameters: Error in syntax') 135 endif 136 alp = fdf_bvalues(pline,1) 137 blp = fdf_bvalues(pline,2) 138 clp = fdf_bvalues(pline,3) 139 alplp = fdf_bvalues(pline,4) 140 betlp = fdf_bvalues(pline,5) 141 gamlp = fdf_bvalues(pline,6) 142 write(6,'(a)') 143 . 'Lattice Parameters (units of Lattice Constant) =' 144 write(6,'(a,3f10.5,3f9.3)') 145 . ' ',alp,blp,clp,alplp,betlp,gamlp 146 alplp = alplp * pi/180.d0 147 betlp = betlp * pi/180.d0 148 gamlp = gamlp * pi/180.d0 149 cell(1,1) = alp 150 cell(2,1) = 0.d0 151 cell(3,1) = 0.d0 152 cell(1,2) = blp * cos(gamlp) 153 cell(2,2) = blp * sin(gamlp) 154 cell(3,2) = 0.d0 155 cell(1,3) = clp * cos(betlp) 156 xxx = (cos(alplp) - cos(betlp)*cos(gamlp))/sin(gamlp) 157 cell(2,3) = clp * xxx 158 cell(3,3) = clp * sqrt(sin(betlp)*sin(betlp) - xxx*xxx) 159 elseif ( fdf_block('LatticeVectors',bfdf) ) then 160 do i = 1,3 161 if (.not. fdf_bline(bfdf, pline)) 162 . call die('redcel: ERROR in LatticeVectors block') 163 cell(1,i) = fdf_bvalues(pline,1) 164 cell(2,i) = fdf_bvalues(pline,2) 165 cell(3,i) = fdf_bvalues(pline,3) 166 enddo 167 else 168 do i = 1,3 169 do j = 1,3 170 cell(i,j) = 0.d0 171 enddo 172 cell(i,i) = 1.d0 173 enddo 174 endif 175 write(6,'(a)') 176 . 'Lattice vectors (in units of Lattice Constant) =' 177 do i = 1,3 178 write(6,'(a,3f10.5)') 179 . ' ',(cell(j,i), j=1,3) 180 enddo 181c ... 182 183c Multiply cell vectors by by lattice constant ... 184 do i = 1,3 185 do ix = 1,3 186 cell(ix,i) = alat * cell(ix,i) 187 enddo 188 enddo 189 write(6,'(a)') 190 . 'Lattice vectors (in Bohr) =' 191 do i = 1,3 192 write(6,'(a,3f10.5)') ' ',(cell(j,i), j=1,3) 193 enddo 194c ... 195 196c Read atomic coordinates and species of unit cell... 197 call recoor(overflow,cell,alat,b,imass,xmass,natoms) 198c ... 199 200c Define number of unit cells in the supercell ... 201 lx_defect = 0 202 ly_defect = 0 203 lz_defect = 0 204 lxmax = fdf_integer('SuperCell_1',lx_defect) 205 lymax = fdf_integer('SuperCell_2',ly_defect) 206 lzmax = fdf_integer('SuperCell_3',lz_defect) 207 call chkdim('fcbuild', 'maxx', maxx, lxmax, 1) 208 call chkdim('fcbuild', 'maxy', maxy, lymax, 1) 209 call chkdim('fcbuild', 'maxz', maxz, lzmax, 1) 210 ncells = (2*lxmax+1)*(2*lymax+1)*(2*lzmax+1) 211 write(6,'(a,i5)') 'lxmax = ',lxmax 212 write(6,'(a,i5)') 'lymax = ',lymax 213 write(6,'(a,i5)') 'lzmax = ',lzmax 214 write(6,'(a,i5)') 'Number of unit cells in Supercell = ',ncells 215c ... 216 217c Read the atomic displacement for Force Constant calculation ... 218 dx = fdf_physical('AtomicDispl',0.04d0,'Bohr') 219 write(6,'(a,f10.5,a)') 'Atomic displacements = ',dx,' Bohr' 220c ... 221 222C *************** END READ DATA FROM FDF FILE ******************** 223 224c Build lattice vector of the supercell ... 225 do ix=1,3 226 scell(ix,1) = (2*lxmax+1)*cell(ix,1)/alat 227 scell(ix,2) = (2*lymax+1)*cell(ix,2)/alat 228 scell(ix,3) = (2*lzmax+1)*cell(ix,3)/alat 229 enddo 230c ... 231 232C Build atomic coordinates in the supercell ... 233c loop over unit cells within supercell 234 iatom=0 235 do lx=-lxmax,lxmax 236 do ly=-lymax,lymax 237 do lz=-lzmax,lzmax 238 r(1) = lx*cell(1,1) + ly*cell(1,2) + lz*cell(1,3) 239 r(2) = lx*cell(2,1) + ly*cell(2,2) + lz*cell(2,3) 240 r(3) = lx*cell(3,1) + ly*cell(3,2) + lz*cell(3,3) 241c loop over atoms in unit cell 242 do i=1,natoms 243 iatom=iatom+1 244 imasssc(iatom)=imass(i) 245 do ix=1,3 246 xa(ix,iatom) = b(ix,i) + r(ix) 247 enddo 248 enddo 249 enddo 250 enddo 251 enddo 252 253 nnat = natoms * (2*lxmax+1) * (2*lymax+1) * (2*lzmax+1) 254 if (iatom .ne. nnat) stop 'Error computing number of atoms' 255c ... 256 257c Determine the indices of the atoms in the central (lx=ly=lz=0) cell 258c (those that need to be displaced to calculate the force constants) ... 259 i1 = natoms * (4*lxmax*lymax*lzmax + 260 . 2*(lxmax*lymax + lxmax*lzmax + lymax*lzmax) + 261 . lxmax + lymax + lzmax) + 1 262 i2 = i1 + natoms - 1 263c ... 264 265C ************** WRITE FDF FORMAT FILE WITH SUPERCELL DATA *********** 266 267c Open file to write ... 268 269 call io_assign(iunit) 270 open(unit=iunit,file='FC.fdf',status='new',err=100) 271 272c Write new constrained file if the GeometryConstraints 273c block exists 274 if ( has_constr ) then 275 write(iunit,*) 276 write(iunit,'(a)') '# GeometryConstraints block found' 277 write(iunit,'(a)') '# defaulting to constrained FC' 278 write(iunit,'(a,a)') 'Vibra.FC ',trim(slabel)//'.FCC' 279 write(iunit,*) 280 end if 281c ... 282 283c Write Number of atoms in Supercell ... 284 write(iunit,'(a,i5)') 'NumberOfAtoms ',nnat 285 write(iunit,*) 286c ... 287 288c Write Lattice Constant... 289 write(iunit,'(a,f15.10,a)') 'LatticeConstant ',alat,' Bohr' 290 write(iunit,*) 291c ... 292 293c Write Lattice vectors of supercell... 294 write(iunit,'(a)') '%block LatticeVectors' 295 do i=1,3 296 write(iunit,10) (scell(ix,i),ix=1,3) 297 enddo 298 write(iunit,'(a)') '%endblock LatticeVectors' 299 write(iunit,*) 300c ... 301 302c Write atomic coordinates format ... 303 write(iunit,'(a)') 304 . 'AtomicCoordinatesFormat NotScaledCartesianBohr' 305 write(iunit,*) 306c ... 307 308c Write atomic coordinates and species in supercell ... 309 write(iunit,'(a)') 310 . '%block AtomicCoordinatesAndAtomicSpecies' 311 do i=1,nnat 312 write(iunit,11) (xa(ix,i),ix=1,3),imasssc(i) 313 enddo 314 write(iunit,'(a)') 315 . '%endblock AtomicCoordinatesAndAtomicSpecies' 316 write(iunit,*) 317c ... 318 319c Write MD options: Force Constant Calculation ... 320 write(iunit,'(2a)') 'MD.TypeOfRun ',' FC' 321 write(iunit,'(a,i5)') 'MD.FCfirst ',i1 322 write(iunit,'(a,i5)') 'MD.FClast ',i2 323 write(iunit,'(a,f15.10,a)') 'MD.FCdispl ',dx,' Bohr' 324c ... 325 326 write(6,'(/,3(/,a,/),/)') 327 . '************************************************************', 328 . ' Output (in FDF format) has been written in FC.fdf', 329 . '************************************************************' 330 33110 format(3(f15.10,2x)) 33211 format(3(f15.10,2x),i5) 333 334 stop 335 336100 continue 337 write(6,'(/,3(/,a,/),/)') 338 . '********************************************', 339 . ' ERROR: File FC.fdf already exists!', 340 . '********************************************' 341 stop 342 343 CONTAINS 344 345 subroutine die(str) 346 character(len=*), intent(in), optional:: str 347 if (present(str)) then 348 write(6,"(a)") str 349 write(0,"(a)") str 350 endif 351 STOP 352 end subroutine die 353 354 end program fcbuild 355 356