Source code for scm.plams.interfaces.thirdparty.crystal

from os.path import join as opj
import numpy as np

from ...core.basejob  import SingleJob
from ...core.settings import Settings
from ...core.errors import PlamsError
from ...mol.molecule  import Molecule

__all__ = ['CrystalJob','mol2CrystalConf']



[docs]class CrystalJob(SingleJob): """ A class representing a single computational job with `CRYSTAL <https://www.crystal.unito.it/>` Use Crystal version >= 14, lower versions have an even stricter and more complicated input. """ _command = 'crystal' _filenames = {'inp':'INPUT', 'run':'$JN.run', 'out':'$JN.out', 'err': '$JN.err'}
[docs] def get_input(self): """ Transform all contents of ``input`` branch of ``settings`` into string. """ #these are the geometry keys from the manual, first mandatory block _geom_keys = ['CRYSTAL','SLAB','POLYMER','HELIX','MOLECULE','EXTERNAL','DLVINPUT'] #second mandatory block: basis set #we use the option to end the geometry section with BASISSET instead of END, #because that is much more like PLAMS and way nicer to code. #Use CUSTOM keyword to give custom basis sets, see Crystal manual _basis_keys = ['BASISSET'] #third mandatory block: Hamiltonian and SCF control #this needs to be after geometry and basis and this is also the last block #so everything else goes in here _option_keys = ['OPTIONS'] def parse(key, value): ret = '' key = key.upper() if isinstance(value, Settings): if key in _geom_keys: #geometry block should be available as a list of lines if not isinstance(value, list): raise PlamsError('Geometry block does not support subblocks') #add geometry lines if '_h' in value: ret += '{}\n'.format(value['_h'].upper()) i = 1 while ('_'+str(i)) in value: ret += parse('',value['_'+str(i)]) i += 1 for el in value: if not el.startswith('_'): ret += parse(el,value[el]) elif key in _basis_keys: for subkey, item in value: #we have a CUSTOM basis set if not subkey.upper() == 'CUSTOM': raise PlamsError('BasisSet block only supports subblock CUSTOM') if '_h' in item: ret += '{}\n'.format(item['_h'].upper()) i = 1 while ('_'+str(i)) in item: ret += parse('',item['_'+str(i)]) i += 1 for el in item: if not el.startswith('_'): ret += parse('',item[el]) ret += 'END\n' #option block has no start key elif key in _option_keys: if '_h' in value: ret += '{}\n'.format(value['_h'].upper()) i = 1 while ('_'+str(i)) in value: ret += parse('',value['_'+str(i)]) i += 1 for el in value: if not el.startswith('_'): ret += parse(el,value[el]) else: #one line per item ret += '{}\n'.format(key) if '_h' in value: ret += '{}\n'.format(value['_h'].upper()) i = 1 while ('_'+str(i)) in value: ret += parse('',value['_'+str(i)]) i += 1 for el in value: if not el.startswith('_'): ret += parse(el,value[el]) ret += 'END\n' elif isinstance(value, str) and key in [*_geom_keys,*_basis_keys,*_option_keys]: ret += '{}\n'.format(value.upper()) elif isinstance(value, list): if key != '': ret += '{}\n'.format(key) for el in value: ret += '{}\n'.format(str(el).upper()) elif key == '': ret += '{}\n'.format(str(value).upper()) elif value == '' or value is True: ret += '{}\n'.format(key) elif value is False: pass else: ret += '{}\n{}\n'.format(key, str(value).upper()) return ret inp = '' use_molecule = ('ignore_molecule' not in self.settings) or (self.settings.ignore_molecule == False) if use_molecule: self._parsemol() #we need a certain ordering, so make a copy of the settings instance #and convert all first-level keys to uppercase tmp = Settings() for key in self.settings.input: tmp[key.upper()] = self.settings.input[key] #check for three blocks if not any([ x in _geom_keys for x in tmp ]): raise PlamsError('One geometry block is necessary for a Crystal Job') if not any([ x in _basis_keys for x in tmp ]): raise PlamsError('BasisSet block is necessary for a Crystal Job') #first title inp += '{}\n'.format(self.name) #geometry block next for item in _geom_keys: if item in tmp: inp += parse(item, tmp[item]) #do not close block, it is closed by the BASISSET keyword following del tmp[item] #basis set block next for item in _basis_keys: if item in tmp: inp += 'BASISSET\n' inp += parse(item, tmp[item]) del tmp[item] #everything else now for item in tmp: inp += parse(item, tmp[item]) inp += 'END' return inp
def _parsemol(self): if 'ase' in Molecule._writeformat: #ASE has a write function for Crystal coordinate files, use that if possible filename = opj(self.path, 'fort.34') #The Crystal IO of ASE only supports filenames, not descirptors right now self.molecule.writease(filename) self.settings.input.external = True else: raise PlamsError('Crystal Interface has no builtin Molecule support, install ASE or use function crystalMol2Conf() and set self.settings.ignore_molecule. See Doc for details.')
[docs] def get_runscript(self): """ Run Crystal. """ ret = self._command ret += ' < ' + self._filename('inp') if self.settings.runscript.stdout_redirect: ret += ' >' + self._filename('out') ret += '\n\n' return ret
[docs] def check(self): """ Look for the normal termination signal in output. Note, that does not mean your calculation was successful! """ termination = self.results.grep_output('TERMINATION') return len(termination) == 1
[docs]def mol2CrystalConf(molecule): """ Call this function to create a CRYSTAL-type input of your structure: Returns a given |Molecule| object as a geomkey and a list of strings that can be used to create a Settings instance for Crystal. >>> print(crystalMol2Conf(mol)) 'GEOMKEY', ['0 0 0', '1', 'lattice', 'nAtoms', 'ElementNumber1 X1 Y1 Z1','ElementNumber2 X2 Y2 Z2', ...] - IFLAG,IFHR and IFSO are returned as 0,0,0 by default with Symmetry group P1 (number 1). This should allow most calculations to run. The user needs to change them if he wants to take advantage of symmetry. - The geometry key is guessed from the number of lattice vectors. For special stuff change it by hand. - The number of lattice vectors in the given molecule should correspond to the dimensionality of the system. Do not fill them with zeros or unit vectors, this will result in a 3D-Periodic system with wrong fractional coordinates. So stick with the standard PLAMS way of doing things. """ geomList = [] #add lattice keyword if len(molecule.lattice) == 0: geomKey = 'MOLECULE' elif len(molecule.lattice) == 1: geomKey = 'POLYMER' elif len(molecule.lattice) == 2: geomKey = 'SLAB' elif len(molecule.lattice) == 3: geomKey = 'CRYSTAL' #add line for IFLAG,IFHR,IFSO: 0 0 0 always works with P1, only for CRYSTAL if geomKey == 'CRYSTAL': geomList.append('0 0 0') #add a line for space group, assume P1 Symmetry because this always works geomList.append('1') #add lattice information: vector lengths and angles (if there is any) lengths = [] angles = [] lattice = molecule.lattice[:] for vec in lattice: lengths.append(np.linalg.norm(vec)) if len(lattice) == 2: v1_u = lattice[0] / np.linalg.norm(lattice[0]) v2_u = lattice[1] / np.linalg.norm(lattice[1]) angle = np.arccos(np.clip(np.dot(v1_u, v2_u), -1.0, 1.0)) / np.pi*180.0 angles.append(angle) elif len(lattice) == 3: for first in range(0,3): second = first - 1 third = first - 2 v1_u = lattice[second] / np.linalg.norm(lattice[second]) v2_u = lattice[third] / np.linalg.norm(lattice[third]) angle = np.arccos(np.clip(np.dot(v1_u, v2_u), -1.0, 1.0)) / np.pi*180.0 angles.append(angle) if len(lattice) > 0: geomList.append((("{} ")*(len(lengths)+len(angles))).format(*lengths,*angles)) #add number of atoms geomList.append(str(len(molecule))) #now add atoms in fractional coordinates or real coords if no lattice #transpose lattice nDim = len(lattice) if nDim > 0: if nDim == 1: lattice.append(np.array([0.0,1.0,0.0])) lattice.append(np.array([0.0,0.0,1.0])) elif nDim == 2: lattice.append(np.array([0.0,0.0,1.0])) #transpose lattice, since PLAMS saves vectors as rows, we need it as columns lattice = np.transpose(lattice) latticeMatInv = np.linalg.inv(lattice) for atom in molecule: atomVec = np.dot(latticeMatInv, np.array(atom.coords)).tolist() #make the crystal coordinates go from 0 to 1 for i in range(0,nDim): if atomVec[i] < 0.0: while atomVec[i] < 0.0: atomVec[i] += 1.0 if atomVec[i] > 1.0: while atomVec[i] > 1.0: atomVec[i] -= 1.0 geomList.append('{:<2} {:>14f} {:>14f} {:>14f}'.format(atom.atnum,*atomVec)) else: for atom in molecule: geomList.append('{:<2} {:}'.format(atom.atnum,atom.str(symbol=False))) return geomKey, geomList