1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38 package ffx.potential.commands;
39
40 import ffx.crystal.Crystal;
41 import ffx.crystal.LatticeSystem;
42 import ffx.potential.bonded.Atom;
43 import ffx.potential.cli.PotentialCommand;
44 import ffx.utilities.FFXBinding;
45 import picocli.CommandLine.Command;
46 import picocli.CommandLine.Option;
47 import picocli.CommandLine.Parameters;
48
49 import java.io.BufferedWriter;
50 import java.io.File;
51 import java.io.FileWriter;
52 import java.io.IOException;
53 import java.util.HashMap;
54 import java.util.Iterator;
55 import java.util.Map;
56
57 import static java.lang.String.format;
58 import static org.apache.commons.io.FilenameUtils.getName;
59 import static org.apache.commons.io.FilenameUtils.removeExtension;
60 import static org.apache.commons.math3.util.FastMath.cos;
61
62
63
64
65
66
67
68 @Command(name = "XYZtoQE", description = "Generate QE input from a XYZ file.")
69 public class ExportQE extends PotentialCommand {
70
71
72 @Option(names = {"--ns", "--nstep"}, paramLabel = "500", defaultValue = "500",
73 description = "Number of structural optimization steps performed in this run.")
74 private int nStep;
75
76
77 @Option(names = {"--ec", "--etot_conv_thr"}, paramLabel = "1.0e-6", defaultValue = "1.0e-6",
78 description = "Convergence threshold on total energy (a.u) for ionic minimization.")
79 private double etotConvThr;
80
81
82 @Option(names = {"--ef", "--forc_conv_thr"}, paramLabel = "1.0e-4", defaultValue = "1.0e-4",
83 description = "Convergence threshold on forces (a.u) for ionic minimization.")
84 private double forcConvThr;
85
86
87 @Option(names = {"--ke", "--ecutwfc"}, paramLabel = "50.0", defaultValue = "50.0",
88 description = "Kinetic energy cutoff (Ry) for wavefunctions.")
89 private double ecutwfc;
90
91
92 @Option(names = {"--rho", "--ecutrho"}, paramLabel = "500.0", defaultValue = "500.0",
93 description = "Kinetic energy cutoff (Ry) for charge density and potential.")
94 private double ecutrho;
95
96
97 @Option(names = {"--em", "--electron_maxstep"}, paramLabel = "1500", defaultValue = "1500",
98 description = "Maximum number of iterations in a scf step.")
99 private int electronMaxstep;
100
101
102 @Option(names = {"--ct", "--conv_thr"}, paramLabel = "1.0e-8", defaultValue = "1.0e-8",
103 description = "Convergence threshold for self consistency.")
104 private double convThr;
105
106
107 @Option(names = {"--mb", "--mixing_beta"}, paramLabel = "0.5", defaultValue = "0.5",
108 description = "Mixing factor for self-consistency.")
109 private double mixingBeta;
110
111
112 @Option(names = {"--hx", "--hexagonal"}, paramLabel = "true", defaultValue = "true",
113 description = "Perform QE on hexagonal system.")
114 private boolean hexagonal;
115
116
117 @Parameters(arity = "1", paramLabel = "file",
118 description = "XYZ file to be converted.")
119 private String filename = null;
120
121 public ExportQE() { super(); }
122 public ExportQE(FFXBinding binding) { super(binding); }
123 public ExportQE(String[] args) { super(args); }
124
125 @Override
126 public ExportQE run() {
127
128 if (!init()) {
129 return this;
130 }
131
132
133 activeAssembly = getActiveAssembly(filename);
134 if (activeAssembly == null) {
135 logger.info(helpString());
136 return this;
137 }
138
139
140 filename = activeAssembly.getFile().getAbsolutePath();
141 logger.info(format("\n Converting %s to QE format\n", filename));
142
143 String dirString = getBaseDirString(filename);
144 String name = removeExtension(getName(filename));
145 File modelFile = new File(dirString + name + ".in");
146
147 Crystal crystal = activeAssembly.getCrystal().getUnitCell();
148 double xtalA = crystal.a;
149 double xtalB = crystal.b;
150 double xtalC = crystal.c;
151
152 HashMap<String, Double> atomTypes = new HashMap<>();
153 StringBuilder atomicPositions = new StringBuilder();
154
155 Atom[] atoms = activeAssembly.getAtomArray();
156 for (Atom atom : atoms) {
157 if (!atomTypes.containsKey(atom.getName())) {
158 atomTypes.put(atom.getName(), atom.getAtomType().atomicWeight);
159 }
160 double[] xyz = atom.getXYZ(null);
161 crystal.toFractionalCoordinates(xyz, xyz);
162 atomicPositions.append(format("%2s %16.12f %16.12f %16.12f%n", atom.getName(), xyz[0], xyz[1], xyz[2]));
163 }
164
165 try (BufferedWriter bwQE = new BufferedWriter(new FileWriter(modelFile))) {
166
167 bwQE.write(format("&CONTROL%n" +
168 "\tcalculation = 'vc-relax',%n" +
169 "\trestart_mode = 'from_scratch',%n" +
170 "\tprefix = '%s',%n" +
171 "\tetot_conv_thr = %6.4E,%n" +
172 "\tforc_conv_thr = %6.4E,%n" +
173 "\tnstep = %d,%n" +
174 "/%n", name, etotConvThr, forcConvThr, nStep));
175
176
177 bwQE.write(format("&SYSTEM%n" +
178 "\tspace_group = %d,%n" +
179 "\tnat = %d,%n" +
180 "\tntyp = %d,%n" +
181 "\ta = %16.12f%n" +
182 "\tb = %16.12f%n" +
183 "\tc = %16.12f%n" +
184 "\tcosAB = %16.12f%n" +
185 "\tcosAC = %16.12f%n" +
186 "\tcosBC = %16.12f%n" +
187 "\tecutwfc = %6.4f,%n" +
188 "\tecutrho = %6.4f,%n" +
189 "\tvdw_corr = 'XDM',%n",
190 crystal.spaceGroup.number,
191 activeAssembly.getAtomList().size(),
192 atomTypes.size(),
193 xtalA, xtalB, xtalC,
194 cos(crystal.gamma), cos(crystal.beta), cos(crystal.alpha),
195 ecutwfc, ecutrho));
196 if (crystal.spaceGroup.latticeSystem == LatticeSystem.HEXAGONAL_LATTICE) {
197 bwQE.write("\trhombohedral = .FALSE.,\n");
198 }
199 bwQE.write("/\n");
200
201
202 bwQE.write(format("&ELECTRONS%n" +
203 "\telectron_maxstep = %d,%n" +
204 "\tconv_thr = %6.4E,%n" +
205 "\tscf_must_converge = .TRUE.,%n" +
206 "\tmixing_beta = %5.3f,%n" +
207 "/%n", electronMaxstep, convThr, mixingBeta));
208
209
210 bwQE.write("&IONS\n" +
211 "\tion_dynamics = 'bfgs',\n" +
212 "/\n");
213
214
215 bwQE.write("&CELL\n" +
216 "\tcell_dynamics = 'bfgs',\n" +
217 "/\n");
218
219
220 StringBuilder line = new StringBuilder();
221 Iterator<Map.Entry<String, Double>> it = atomTypes.entrySet().iterator();
222 while (it.hasNext()) {
223 Map.Entry<String, Double> pair = it.next();
224 line.append(" ").append(pair.getKey()).append(" ").append(pair.getValue())
225 .append(" ").append(pair.getKey()).append(".b86bpbe.UPF\n");
226 }
227 bwQE.write("ATOMIC_SPECIES\n" + line + "\n");
228
229
230 bwQE.write("ATOMIC_POSITIONS crystal_sg\n" + atomicPositions + "\n");
231
232
233 int k1; int k2; int k3;
234 if (xtalA < 5) { k1 = 8; } else if (xtalA <= 7) { k1 = 6; } else if (xtalA <= 12) { k1 = 4; } else { k1 = 2; }
235 if (xtalB < 5) { k2 = 8; } else if (xtalB <= 7) { k2 = 6; } else if (xtalB <= 12) { k2 = 4; } else { k2 = 2; }
236 if (xtalC < 5) { k3 = 8; } else if (xtalC <= 7) { k3 = 6; } else if (xtalC <= 12) { k3 = 4; } else { k3 = 2; }
237 bwQE.write("K_POINTS automatic\n" + k1 + " " + k2 + " " + k3 + " 1 1 1\n");
238 } catch (IOException e) {
239 logger.severe(" Error writing QE input file: " + e.getMessage());
240 }
241
242 logger.info(format(" Saved QE file: %s", modelFile));
243 return this;
244 }
245 }