1 /*------------------------------------------------------------------- 2 Copyright 2011 Ravishankar Sundararaman, Deniz Gunceler 3 4 This file is part of JDFTx. 5 6 JDFTx is free software: you can redistribute it and/or modify 7 it under the terms of the GNU General Public License as published by 8 the Free Software Foundation, either version 3 of the License, or 9 (at your option) any later version. 10 11 JDFTx is distributed in the hope that it will be useful, 12 but WITHOUT ANY WARRANTY; without even the implied warranty of 13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14 GNU General Public License for more details. 15 16 You should have received a copy of the GNU General Public License 17 along with JDFTx. If not, see <http://www.gnu.org/licenses/>. 18 -------------------------------------------------------------------*/ 19 20 #include <commands/command.h> 21 #include <electronic/Everything.h> 22 #include <fluid/Euler.h> 23 24 struct CommandIon : public Command 25 { CommandIonCommandIon26 CommandIon() : Command("ion", "jdftx/Ionic/Geometry") 27 { 28 format = "<species-id> <x0> <x1> <x2> [v <vx0> <vx1> <vx2>] <moveScale> [<constraint type>=" 29 + constraintTypeMap.optionList() + " <d0> <d1> <d2> [<group> [HyperPlane <d0> ...]]]"; 30 comments = 31 "Add an atom of species <species-id> at coordinates (<x0>,<x1>,<x2>).\n" 32 "\n" 33 "Optionally, for dynamics, specify ion velocity <v0>,<v1>,<v2> after keyword 'v'.\n" 34 "\n" 35 "<moveScale> preconditions the motion of this ion (set 0 to hold fixed)\n" 36 "\n" 37 "In addition, the ion may be constrained to a line or a plane with line\n" 38 "direction or plane normal equal to (<d0>,<d1>,<d2>) in the coordinate\n" 39 "system selected by command coords-type. Note that the constraints must\n" 40 "be consistent with respect to symmetries (if enabled).\n" 41 "\n" 42 "The HyperPlane constraint allows constraining collective motion of many\n" 43 "ions by restricting their motion to a hyperplane with normal specified\n" 44 "by (<d0>,<d1>,<d2>) for all ions specifying a hyperplane constraint.\n" 45 "By default, all hyperplane-constrained ions are included in a single\n" 46 "group; use optional <group> label to specify multiple hyper-planes.\n" 47 "Multiple hyperplane constraints may also be added to each atom,\n" 48 "but this requires an explicit group label for each hyperplane.\n" 49 "\n" 50 "Note that when coords-type is lattice, the constraints are in covariant\n" 51 "lattice coordinates (like direction indices) for line constraints, but\n" 52 "contravariant coordinates (like plane indices) for plane constraints."; 53 54 allowMultiple = true; 55 56 require("ion-species"); 57 //Dependencies due to coordinate system option: 58 require("latt-scale"); 59 require("coords-type"); 60 } 61 processCommandIon62 void process(ParamList& pl, Everything& e) 63 { //Find species: 64 string id; pl.get(id, string(), "species-id", true); 65 auto sp = findSpecies(id, e); 66 if(!sp) throw string("Species "+id+" has not been defined"); 67 //Read coordinates: 68 vector3<> pos; 69 for(int k=0; k<3; k++) 70 { ostringstream oss; oss << "x" << k; 71 pl.get(pos[k], 0., oss.str(), true); 72 } 73 //Transform coordinates if necessary 74 if(e.iInfo.coordsType == CoordsCartesian) 75 pos = inv(e.gInfo.R) * pos; 76 //Add position to list: 77 sp->atpos.push_back(pos); 78 sp->velocities.push_back(vector3<>(NAN,NAN,NAN)); 79 80 //Look for optional velocity and get moveScale: 81 string key; 82 pl.get(key, string(), "moveScale|v", true); 83 SpeciesInfo::Constraint constraint; 84 if(key == "v") 85 { //Read velocity (update the NAN added above): 86 vector3<>& vel = sp->velocities.back(); 87 for(int k=0; k<3; k++) 88 { ostringstream oss; oss << "v" << k; 89 pl.get(vel[k], 0., oss.str(), true); 90 } 91 //Transform coordinates if necessary 92 if(e.iInfo.coordsType == CoordsCartesian) 93 vel = inv(e.gInfo.R) * vel; 94 //Get moveScale from command line beyond the velocity: 95 pl.get(constraint.moveScale, 0., "moveScale", true); 96 } 97 else 98 { //No velocity: get moveScale from key 99 ParamList(key+" ").get(constraint.moveScale, 0., "moveScale", true); 100 } 101 102 //Parse constraints 103 if(constraint.moveScale < 0.) // Check for negative moveScales 104 throw string("moveScale cannot be negative"); 105 pl.get(constraint.type, SpeciesInfo::Constraint::None, constraintTypeMap, "Type"); 106 if(constraint.type != SpeciesInfo::Constraint::None) 107 { if(!constraint.moveScale) 108 throw string("Constraint specified after moveScale = 0"); 109 pl.get(constraint.d[0], 0.0, "d0", true); 110 pl.get(constraint.d[1], 0.0, "d1", true); 111 pl.get(constraint.d[2], 0.0, "d2", true); 112 if(e.iInfo.coordsType == CoordsLattice) //Transform to Cartesian (taking care of covariant/contravariant for line/plane directions) 113 switch(constraint.type) 114 { case SpeciesInfo::Constraint::Linear: constraint.d = e.gInfo.R * constraint.d; break; 115 case SpeciesInfo::Constraint::Planar: constraint.d = ~inv(e.gInfo.R) * constraint.d; break; 116 case SpeciesInfo::Constraint::HyperPlane: constraint.d = ~inv(e.gInfo.R) * constraint.d; break; 117 default: break; 118 } 119 if(constraint.type == SpeciesInfo::Constraint::HyperPlane) 120 { string groupLabel; 121 pl.get(groupLabel, string(), "group"); //optional group label for hyperplane constraint 122 constraint.hyperplane.assign(1, std::make_pair(constraint.d, groupLabel)); 123 constraint.d = vector3<>(); //zero out the constraint d to prevent incorrect symmetry checks 124 //Look for additional hyperplane constraints: 125 while(true) 126 { string key; pl.get(key, string(), "Type"); 127 if(not key.length()) break; //End of constraint list 128 if(key != "HyperPlane") throw string("Additional constraints must be of type HyperPlane"); 129 vector3<> d; 130 pl.get(d[0], 0.0, "d0", true); 131 pl.get(d[1], 0.0, "d1", true); 132 pl.get(d[2], 0.0, "d2", true); 133 pl.get(groupLabel, string(), "group", true); //group label required for additional hyperplanes 134 if(e.iInfo.coordsType == CoordsLattice) d = ~inv(e.gInfo.R) * d; //transform to Cartesian 135 constraint.hyperplane.push_back(std::make_pair(d, groupLabel)); 136 } 137 } 138 else //Note d can be null for a guiven atom with HyperPlane 139 { if(not constraint.d.length_squared()) 140 throw string("Constraint vector must be non-null"); 141 } 142 } 143 sp->constraints.push_back(constraint); 144 } 145 printStatusCommandIon146 void printStatus(Everything& e, int iRep) 147 { int iIon=0; 148 for(auto sp: e.iInfo.species) 149 for(unsigned at=0; at<sp->atpos.size(); at++) 150 { if(iIon==iRep) 151 { //Species and position: 152 vector3<> pos = sp->atpos[at]; 153 if(e.iInfo.coordsType == CoordsCartesian) 154 pos = e.gInfo.R * pos; //report cartesian positions 155 logPrintf("%s %19.15lf %19.15lf %19.15lf", sp->name.c_str(), pos[0], pos[1], pos[2]); 156 //Optional velocity: 157 vector3<> vel = sp->velocities[at]; 158 if(not std::isnan(vel.length_squared())) 159 { if(e.iInfo.coordsType == CoordsCartesian) 160 vel = e.gInfo.R * vel; //report cartesian velocities 161 logPrintf(" v %19.15lf %19.15lf %19.15lf", vel[0], vel[1], vel[2]); 162 } 163 logPrintf(" %lg", sp->constraints[at].moveScale); 164 if(sp->constraints[at].type != SpeciesInfo::Constraint::None) 165 sp->constraints[at].print(globalLog, e); 166 } 167 iIon++; 168 } 169 } 170 } 171 commandIon; 172 173 //------------------------------------------------------------------------------------------------- 174 175 struct CommandInitialMagneticMoments : public Command 176 { CommandInitialMagneticMomentsCommandInitialMagneticMoments177 CommandInitialMagneticMoments() : Command("initial-magnetic-moments", "jdftx/Initialization") 178 { 179 format = "<species> <M1> <M2> ... <Mn> [<species2> ...]\n" 180 " | <species> <M1> <theta1> <phi1> ... <Mn> <thetan> <phin> [<species2> ...]"; 181 comments = 182 "Specify initial magnetic moments, defined as the difference between\n" 183 "up and down electron counts, on each atom of one or more species.\n" 184 "\n" 185 "The second syntax with polar angles (in degrees) must be used\n" 186 "for noncollinear magnetism calculations.\n" 187 "\n" 188 "For each species, the initial magnetic moments are applied\n" 189 "to the atoms in the order of ion commands for that species.\n" 190 "This may be used to construct a spin-polarized reference density\n" 191 "for LCAO initialization of the Kohn-Sham orbitals."; 192 193 require("ion"); 194 require("spintype"); 195 } 196 processCommandInitialMagneticMoments197 void process(ParamList& pl, Everything& e) 198 { if(e.eInfo.spinType==SpinNone || e.eInfo.spinType==SpinOrbit) 199 throw string("Cannot specify magnetic moments in an unpolarized calculation"); 200 string id; 201 pl.get(id, string(), "species", true); 202 while(id.length()) 203 { auto sp = findSpecies(id, e); 204 if(!sp) throw string("Species "+id+" has not been defined"); 205 sp->initialMagneticMoments.resize(sp->atpos.size()); 206 for(unsigned a=0; a<sp->atpos.size(); a++) 207 { double M=0., theta=0., phi=0.; 208 pl.get(M, 0., "M", true); 209 if(e.eInfo.spinType == SpinVector) 210 { pl.get(theta, 0., "theta", true); 211 pl.get(phi, 0., "phi", true); 212 } 213 //Store as Cartesian: 214 sp->initialMagneticMoments[a] = M * polarUnitVector(phi*M_PI/180, theta*M_PI/180); 215 } 216 //Check for additional species: 217 pl.get(id, string(), "species"); 218 } 219 } 220 printStatusCommandInitialMagneticMoments221 void printStatus(Everything& e, int iRep) 222 { for(auto sp: e.iInfo.species) 223 if(sp->initialMagneticMoments.size()) 224 { logPrintf(" \\\n\t%s", sp->name.c_str()); 225 for(const vector3<>& M: sp->initialMagneticMoments) 226 { if(e.eInfo.spinType == SpinVector) 227 { vector3<> euler; 228 if(M.length()) getEulerAxis(M, euler); 229 euler *= 180./M_PI; //convert to degrees 230 logPrintf(" %lg %lg %lg ", M.length(), euler[1], euler[0]); 231 } 232 else logPrintf(" %lg", M[2]); //SpinZ 233 } 234 } 235 } 236 } 237 commandInitialMagneticMoments; 238 239 //------------------------------------------------------------------------------------------------- 240 241 struct CommandInitialOxidationState : public Command 242 { CommandInitialOxidationStateCommandInitialOxidationState243 CommandInitialOxidationState() : Command("initial-oxidation-state", "jdftx/Initialization") 244 { 245 format = "<species> <oxState> [<species2> ...]"; 246 comments = 247 "Specify initial oxidation state assumed for each species in LCAO.\n" 248 "This may significantly improve LCAO convergence in some cases."; 249 250 require("ion-species"); 251 } 252 processCommandInitialOxidationState253 void process(ParamList& pl, Everything& e) 254 { string id; 255 pl.get(id, string(), "species", true); 256 while(id.length()) 257 { auto sp = findSpecies(id, e); 258 if(!sp) throw string("Species "+id+" has not been defined"); 259 pl.get(sp->initialOxidationState, 0., "oxState", true); 260 //Check for additional species: 261 pl.get(id, string(), "species"); 262 } 263 } 264 printStatusCommandInitialOxidationState265 void printStatus(Everything& e, int iRep) 266 { for(auto sp: e.iInfo.species) 267 if(sp->initialOxidationState) 268 logPrintf(" \\\n\t%s %lg", sp->name.c_str(), sp->initialOxidationState); 269 } 270 } 271 commandInitialOxidationState; 272 273 EnumStringMap<coreOverlapCheck> overlapCheckMap 274 ( additive, "additive", 275 vector, "vector", 276 none, "none" 277 ); 278 279 280 struct CommandCoreOverlapCheck : public Command 281 { CommandCoreOverlapCheckCommandCoreOverlapCheck282 CommandCoreOverlapCheck() : Command("core-overlap-check", "jdftx/Ionic/Optimization") 283 { 284 format = "<condition>"; 285 comments = "Checks for core overlaps between ionic pseudopotentials based on <condition>:\n" 286 "+ additive: checks for interatomic distance < (R1 + R2)\n" 287 "+ vector: checks for interatomic distance < sqrt(R1^2 + R2^2) (default)\n" 288 "+ none"; 289 hasDefault = true; 290 } 291 processCommandCoreOverlapCheck292 void process(ParamList& pl, Everything& e) 293 { pl.get(e.iInfo.coreOverlapCondition, vector, overlapCheckMap, "overlap check"); 294 } 295 printStatusCommandCoreOverlapCheck296 void printStatus(Everything& e, int iRep) 297 { logPrintf("%s", overlapCheckMap.getString(e.iInfo.coreOverlapCondition)); 298 } 299 } 300 commandCoreOverlapCheck; 301