1#!/usr/bin/env ruby
2
3# Time-stamp: <2015-11-19 17:23:47 togo>
4#
5#   Copyright (C) 2005 Atsushi Togo
6#   togo.atsushi@gmail.com
7#
8#   This program is free software; you can redistribute it and/or
9#   modify it under the terms of the GNU General Public License
10#   as published by the Free Software Foundation; either version 2
11#   of the License, or (at your option) any later version.
12#
13#   This program is distributed in the hope that it will be useful,
14#   but WITHOUT ANY WARRANTY; without even the implied warranty of
15#   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16#   GNU General Public License for more details.
17#
18#   You should have received a copy of the GNU General Public License
19#   along with this program; if not, write to
20#   the Free Software Foundation, Inc., 51 Franklin Street,
21#   Fifth Floor, Boston, MA 02110-1301, USA, or see
22#   http://www.gnu.org/copyleft/gpl.html
23
24require "optparse"
25
26module Crystal
27  # constant
28  AtomicNumber = {"H"=>1, "He"=>2, "Li"=>3, "Be"=>4, "B"=>5, "C"=>6, "N"=>7, "O"=>8,
29    "F"=>9, "Ne"=>10, "Na"=>11, "Mg"=>12, "Al"=>13, "Si"=>14, "P"=>15, "S"=>16,
30    "Cl"=>17, "Ar"=>18, "K"=>19, "Ca"=>20, "Sc"=>21, "Ti"=>22, "V"=>23, "Cr"=>24,
31    "Mn"=>25, "Fe"=>26, "Co"=>27, "Ni"=>28, "Cu"=>29, "Zn"=>30, "Ga"=>31, "Ge"=>32,
32    "As"=>33, "Se"=>34, "Br"=>35, "Kr"=>36, "Rb"=>37, "Sr"=>38, "Y"=>39, "Zr"=>40,
33    "Nb"=>41, "Mo"=>42, "Tc"=>43, "Ru"=>44, "Rh"=>45, "Pd"=>46, "Ag"=>47, "Cd"=>48,
34    "In"=>49, "Sn"=>50, "Sb"=>51, "Te"=>52, "I"=>53, "Xe"=>54, "Cs"=>55, "Ba"=>56,
35    "La"=>57, "Ce"=>58, "Pr"=>59, "Nd"=>60, "Pm"=>61, "Sm"=>62, "Eu"=>63, "Gd"=>64,
36    "Tb"=>65, "Dy"=>66, "Ho"=>67, "Er"=>68, "Tm"=>69, "Yb"=>70, "Lu"=>71, "Hf"=>72,
37    "Ta"=>73, "W"=>74, "Re"=>75, "Os"=>76, "Ir"=>77, "Pt"=>78, "Au"=>79, "Hg"=>80,
38    "Tl"=>81, "Pb"=>82, "Bi"=>83, "Po"=>84, "At"=>85, "Rn"=>86, "Fr"=>87, "Ra"=>88,
39    "Ac"=>89, "Th"=>90, "Pa"=>91, "U"=>92, "Np"=>93, "Pu"=>94, "Am"=>95, "Cm"=>96,
40    "Bk"=>97, "Cf"=>98, "Es"=>99, "Fm"=>100, "Md"=>101, "No"=>102, "Lr"=>103, "Rf"=>104,
41    "Db"=>105, "Sg"=>106, "Bh"=>107, "Hs"=>108, "Mt"=>109, "Ds"=>110, "Rg"=>111}
42
43  # method
44  def transform2axis(alpha,beta,gamma,cellA,cellB,cellC)
45    alpha = alpha / 180 * Math::PI
46    beta  = beta  / 180 * Math::PI
47    gamma = gamma / 180 * Math::PI
48#     cz = cellC
49#     bz = Math.cos(alpha) * cellB
50#     by = Math.sin(alpha) * cellB
51#     az = Math.cos(beta)  * cellA
52#     ay = ( Math.cos(gamma) - Math.cos(beta) * Math.cos(alpha) ) / Math.sin(alpha) * cellA
53#     ax = Math.sqrt(cellA**2 - ay**2 - az**2)
54#     return [ax,ay,az], [0,by,bz], [0,0,cz]
55    ax = cellA
56    bx = Math.cos(gamma) * cellB
57    by = Math.sin(gamma) * cellB
58    cx = Math.cos(beta)  * cellC
59    cy = ( Math.cos(alpha) - Math.cos(beta) * Math.cos(gamma) ) / Math.sin(gamma) * cellC
60    cz = Math.sqrt(cellC**2 - cy**2 - cx**2)
61    return [ax,0,0], [bx,by,0], [cx,cy,cz]
62  end
63
64  def axisAngle(inputAxis)
65    axis = []
66    3.times {|i| axis << Vector.elements(inputAxis[i])}
67    alpha = Math.acos(axis[1].inner_product(axis[2])/axis[1].r/axis[2].r)/Math::PI*180
68    beta  = Math.acos(axis[2].inner_product(axis[0])/axis[2].r/axis[0].r)/Math::PI*180
69    gamma = Math.acos(axis[0].inner_product(axis[1])/axis[0].r/axis[1].r)/Math::PI*180
70    return alpha, beta, gamma
71  end
72
73  def axisLength(axis)
74    length = []
75    axis.each do |vec|
76      length << Vector.elements(vec).r
77    end
78    length
79  end
80
81  module_function :transform2axis, :axisAngle, :axisLength
82
83  # class
84  class Atom
85    attr_accessor :position
86    def initialize(position)
87      @position = position
88    end
89  end
90
91  # This class is mainly considered in Cartesian coordinate.
92  class Cell
93    require 'matrix'
94
95    attr_accessor :axis, :atoms, :comment
96    def initialize(axis, atoms, comment = nil)
97      @axis = axis
98      @atoms = atoms
99      @comment = comment
100    end
101
102    def axisAngle
103      Crystal.axisAngle(@axis)
104    end
105
106    def axisLength
107      Crystal.axisLength(@axis)
108    end
109
110    def volume
111      a = @axis
112      a[0][0]*(a[1][1]*a[2][2]-a[1][2]*a[2][1])+a[0][1]*(a[1][2]*a[2][1]-a[1][0]*a[2][2])+a[0][2]*(a[1][0]*a[2][1]-a[1][1]*a[2][0])
113    end
114
115    def distance(num1, num2)
116       (Matrix.rows(@axis).transpose * (Vector.elements(@atoms[num1].position) - Vector.elements(@atoms[num2].position))).r
117    end
118
119    def minDistance(num1, num2)
120      min = distance(num1, num2)
121      (-1..1).each do |a|
122        (-1..1).each do |b|
123          (-1..1).each do |c|
124            position = Vector.elements(@atoms[num2].position) + Vector[a,b,c]
125            distance = (Matrix.rows(@axis).transpose * (Vector.elements(@atoms[num1].position) - position )).r
126            min = distance if min > distance
127          end
128        end
129      end
130      min
131    end
132
133    def positionReal(position)
134      (Matrix.rows(@axis).transpose * Vector.elements(position)).to_a
135    end
136  end
137
138end
139
140module Vasp
141
142  class CellToPoscar
143
144    include Crystal
145
146    def initialize(cell)
147      @cell = cell
148    end
149
150    def print
151      puts createPoscar
152    end
153
154    def write(filename)
155      File.open(filename, "w") {|file| file.puts createPoscar}
156    end
157
158    private
159
160    def createPoscar
161      @cell.comment = "POSCAR generated by cell class" if @cell.comment == nil
162      poscar = sprintf("#{@cell.comment.strip}\n")       # comment
163      poscar << (sprintf "1.0\n")                        # scale
164      @cell.axis.each {|row| poscar << sprintf("%20.16f  %20.16f  %20.16f\n", row[0], row[1], row[2])} # axis
165      # number of atoms & sort by atom name
166      names = []
167      @cell.atoms.each {|atom| names << atom.name}
168      atomHash = Hash.new
169      names.uniq.each {|name| atomHash[name] = []}
170      @cell.atoms.each {|atom| atomHash[atom.name] << atom}
171      names.uniq.each {|name| poscar << sprintf(" #{atomHash[name].size}")}
172      poscar << sprintf("\n")
173      poscar << sprintf("Direct\n")                      # Direct only
174      # position
175      names.uniq.each do |name|
176        counter = 0
177        atomHash[name].each do |atom|
178          counter += 1
179          3.times do |i|
180            if sprintf("%20.16f", atom.position[i]).to_f.abs < 1e-16
181              atom.position[i] = 0.0
182            end
183          end
184          poscar << sprintf("%20.16f  %20.16f  %20.16f # #{name}#{counter}\n", atom.position[0], atom.position[1], atom.position[2])
185        end
186      end
187      poscar
188    end
189  end
190
191  class Poscar
192
193    def initialize(filename = "POSCAR", potcarName = "POTCAR")
194      parse(filename, potcarName)
195    end
196
197    def cell
198      Crystal::Cell.new(@axis, @atoms, @comment)
199    end
200
201    private
202
203    def atomNameByPotcar(filename)  # pick up atom names from POTCAR
204      name = []
205      open(filename, "r").each do |line|
206        if line =~ /VRHFIN\s*=\s*([A-Za-z]*)/
207          name << $1
208        end
209      end
210      name
211    end
212
213    def axisPoscar(scale)
214      axis = []
215      3.times {|i| axis[i] = @input.readline.strip.split(/\s+/)} # a,b,c axis
216      # multiply universal scaling factor
217      3.times do |i|
218        3.times do |j|
219          axis[i][j] = axis[i][j].to_f * scale
220        end
221      end
222      axis
223    end
224
225    def namePoscar(numAtoms, potcarName)
226      if File.exist?(potcarName)
227        atomName = atomNameByPotcar(potcarName)
228      else
229        atomName = "ABCDEFGHIJKLMNOPQRSTUVWXYZ".split(//)
230      end
231      name = []
232      numAtoms.size.times do |i|
233        numAtoms[i].times {|j| name << "#{atomName[i]}"}
234      end
235      name
236    end
237
238    def numAtomsPoscar
239      num = @input.readline.strip.split(/\s+/)
240      num.size.times {|i| num[i] = num[i].to_i}
241      num
242    end
243
244    def parse(filename, potcarName)
245      @input = open(filename, "r")
246      @comment = @input.readline.chomp # line 1: comment (string)
247      scale = scalePoscar(filename)    # line 2: universal scaling factor (float)
248      @axis = axisPoscar(scale)        # line 3-5: axis (3x3: float)
249      numAtoms = numAtomsPoscar        # line 6: number of atoms ([] integer), atom name ([name, name, ...])  example [Sn, Sn, O, O]
250      if numAtoms[0] == 0
251        numAtoms = numAtomsPoscar
252      end
253      name = namePoscar(numAtoms, potcarName)
254      #
255      # line 7-(8): 'Selective dynamics' or not (bool)
256      #
257      @selectiveDynamics = false
258      temp = @input.readline
259      if temp =~ /^\s*s/i then           # check 's' or 'S' of the first word
260        @selectiveDynamics = true
261        temp = @input.readline           # when this situation, reading one more line is nessesarry
262      end
263      if (! temp =~ /^\s*d/i)
264        # allow only 'Direct' now
265        warn "#{filename}: line 7 or 8"
266        warn "poscar.rb can handle only direct coordinates."
267        exit
268      end
269
270      @atoms = positionPoscar(numAtoms, name)    # line 9(8): [Atom, ...]
271      # initial states of MD is ignored.
272      @input.close
273    end
274
275    def positionPoscar(numEachAtom, name)
276      atoms = []
277      numEachAtom.each do |num|
278        num.times do
279          lineArr = @input.gets.strip.split(/\s+/)
280          position = lineArr[0..2].collect! {|x| x.to_f}
281          if lineArr.size >= 5 then
282            movable = lineArr[3..5].collect! {|x| x = true if x =~ /^t/i}
283          else
284            movable = [false, false, false]
285          end
286          atoms << Atom.new(position, name.shift, movable)
287        end
288      end
289      atoms
290    end
291
292    def scalePoscar(filename)
293      scale = @input.readline.to_f
294      if scale < 0 # minus value is neglected
295        print "#{filename} line 2\n"
296        print "poscar.rb cannot use negative scaling factor.\n"
297        exit
298      end
299      scale
300    end
301  end
302
303  class Atom < Crystal::Atom
304    attr_accessor :name, :movable
305    def initialize(position, name = nil, movable = nil)
306      @position = position
307      @name = name
308      @movable = movable
309    end
310  end
311end
312
313
314def printFindsym(cell, precision)
315  puts "poscar2findsym"
316  puts "#{precision}"
317  puts "1"
318  cell.axis.each do |vec|
319    printf("%20.16f %20.16f %20.16f\n", vec[0], vec[1], vec[2])
320  end
321  puts "1"
322  puts "1 0 0"
323  puts "0 1 0"
324  puts "0 0 1"
325  puts "#{cell.atoms.size}"
326  # number of atoms & sort by atom name
327  names = []
328  cell.atoms.each {|atom| names << atom.name}
329  names.uniq.each_with_index do |name, i|
330    cell.atoms.each do |atom|
331      if atom.name == name
332        print "#{i+1} "
333      end
334    end
335  end
336  puts
337  names.uniq.each_with_index do |name, i|
338    cell.atoms.each do |atom|
339      if atom.name == name
340        vec = atom.position
341        printf("%20.16f %20.16f %20.16f\n", vec[0], vec[1], vec[2])
342      end
343    end
344  end
345
346  printf("\n")
347end
348
349precision = 1e-5
350opt = OptionParser.new
351opt.on('--precision=', '-p', 'precision') {|tmp| precision = tmp.to_f}
352opt.parse!(ARGV)
353
354printFindsym(Vasp::Poscar.new(ARGV.shift).cell, precision)
355