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