diff --git a/dependencies/matid b/dependencies/matid index 015f6985283260c2da173d651fe3a1cccc4968c5..ac8f7a8e9f16f6a84808ebde5c1f0cd317a686f4 160000 --- a/dependencies/matid +++ b/dependencies/matid @@ -1 +1 @@ -Subproject commit 015f6985283260c2da173d651fe3a1cccc4968c5 +Subproject commit ac8f7a8e9f16f6a84808ebde5c1f0cd317a686f4 diff --git a/gui/public/env.js b/gui/public/env.js index ac192bb1524cd4d87956abaf931bef9955850a4a..5b7b471502d7218fb530d199067fe910c34b0978 100644 --- a/gui/public/env.js +++ b/gui/public/env.js @@ -107,6 +107,14 @@ window.nomadEnv = { "label": "Entry type", "align": "left" }, + "upload_name": { + "label": "Upload name", + "align": "left" + }, + "upload_id": { + "label": "Upload id", + "align": "left" + }, "upload_create_time": { "label": "Upload time", "align": "left" diff --git a/gui/src/components/entry/properties/MaterialCardTopology.js b/gui/src/components/entry/properties/MaterialCardTopology.js index 995f210677ca81eb754bee67845e26bec18f5d8a..f3f7840ee5c758ed6eeb825d5172f89bad201290 100644 --- a/gui/src/components/entry/properties/MaterialCardTopology.js +++ b/gui/src/components/entry/properties/MaterialCardTopology.js @@ -241,12 +241,12 @@ const useTopologyItemStyles = makeStyles({ justifyContent: 'center', flexDirection: 'column', padding: props.theme.spacing(0.25, 0), - flexGrow: 1 + flexGrow: 1, + height: props.theme.spacing(5) }), nodeLabelPrimary: (props) => ({ textTransform: 'uppercase', - fontSize: 14, - marginBottom: -props.theme.spacing(0.3) + fontSize: 14 }), nodeLabelSecondary: (props) => ({ fontSize: 12 @@ -281,9 +281,11 @@ const TopologyItem = React.memo(({node, level, selected}) => { color: isSelected ? 'white' : theme.palette.text.secondary }} > - {`${node.structural_type}` + (!isEmpty(node?.child_systems) && node.structural_type === 'group' - ? ` (${[...new Set(node.child_systems.map(x => x.structural_type))].join(', ')})` - : '' + {node.label === 'original' + ? undefined + : `${node.structural_type}` + (!isEmpty(node?.child_systems) && node.structural_type === 'group' + ? ` (${[...new Set(node.child_systems.map(x => x.structural_type))].join(', ')})` + : '' )} </Typography> </div> diff --git a/gui/src/components/visualization/StructureNGL.js b/gui/src/components/visualization/StructureNGL.js index 6c6ccd0101d55bf5f42722b2a38ae4ff16f58858..8cdf3768ce17435830b2797237e31bcdf5776871 100644 --- a/gui/src/components/visualization/StructureNGL.js +++ b/gui/src/components/visualization/StructureNGL.js @@ -115,7 +115,7 @@ const StructureNGL = React.memo(({ const color = elementColors[symbolUpperCase] return { label: symbol, - color: color ? `#${color.toString(16)}` : '#ff1493', + color: color ? `#${color.toString(16).padStart(6, 0)}` : '#ff1493', radius: vdwRadii[atomicNumber] || 0.1, atomicNumber: atomicNumber || 0 } @@ -334,20 +334,40 @@ const StructureNGL = React.memo(({ // The objects are added to a dummy 'unitcell' representation that is // supported by NGL. const nglCell = component?.object?.unitcell - const collapsed = [false, false, false] - const pbc = [true, true, true] if (nglCell) { + const pbc = topologyMap[root]?.cell?.pbc || [true, true, true] component.addRepresentation('unitcell', {opacity: 0}) const fracToCart = nglCell.fracToCart + fracToCart.elements = nglCell.fracToCart.elements.map((e) => isNaN(e) ? 0 : e) const a = new THREE.Vector3(1, 0, 0).applyMatrix4(fracToCart) const b = new THREE.Vector3(0, 1, 0).applyMatrix4(fracToCart) const c = new THREE.Vector3(0, 0, 1).applyMatrix4(fracToCart) + + // If some of the basis vectors are collapsed, we need to create + // artificial ones for the wrapping + const collapsed = [nglCell.a === 0, nglCell.b === 0, nglCell.c === 0] const basis = [a, b, c] + const validBases = [] + const invalidBases = [] + for (let i = 0; i < collapsed.length; ++i) { + const isCollapsed = collapsed[i] + if (!isCollapsed) validBases.push(basis[i]) + else invalidBases.push(i) + } + if (validBases.length === 2) { + const newBasis = new THREE.Vector3().crossVectors( + validBases[0], + validBases[1] + ).normalize() + basis[invalidBases[0]] = newBasis + } + const cartToFrac = new THREE.Matrix4().makeBasis(basis[0], basis[1], basis[2]).clone().invert() + nglCell.cartToFrac = cartToFrac + const cell = createCell( new THREE.Vector3(), basis, collapsed, - pbc, '#000', 1, 0, @@ -372,6 +392,8 @@ const StructureNGL = React.memo(({ ) addObject3DToStage(latticeConstants, stageRef.current) component.unitCell = cell + component.cartToFrac = cartToFrac + component.fracToCart = fracToCart component.basis = basis component.pbc = pbc component.latticeConstants = latticeConstants @@ -420,26 +442,22 @@ const StructureNGL = React.memo(({ const topParent = topologyMap[selection].parent_system const structuralType = topSelection.structural_type const isMonomer = structuralType === 'monomer' - const isGroup = structuralType === 'group' - const isSubSystem = structuralType === 'subsystem' + const isMolecule = structuralType === 'molecule' + const isGroup = structuralType === 'group' || topSelection.label === 'subsystem' + const independent = topSelection.atoms || topSelection.atoms_ref || isMolecule || isMonomer const child_types = topSelection.child_systems ? new Set(topSelection.child_systems.map(x => x.structural_type)) : new Set() const isMonomerGroup = isGroup && isEqual(child_types, new Set(['monomer'])) - const isMolecule = structuralType === 'molecule' - const showParent = isGroup || isSubSystem // Determine the selection to center on. - selectionRef.current = representationMap.current[(isMonomerGroup - ? topParent - : selection - )]?.sele + selectionRef.current = representationMap.current[independent ? selection : topParent]?.sele // Determine the selections to show opaque const opaque = new Set([selection]) // Determine the selections to show transparent - const transparent = new Set(showParent ? [topParent] : []) + const transparent = new Set(isGroup ? [topParent] : []) // Determine whether to show cell const cellVisible = !(isMolecule || isMonomer || isMonomerGroup) @@ -607,7 +625,6 @@ function getAlignment(alignments, directions) { * @param origin - The origin position for the cell * @param basis - The cell basis vectors * @param collapsed - Whether a basis vector is collapsed (length = 0) - * @param periodicity - Periodicity of the cell as three boolean, one for each basis vector * @param color - Color for the cell wireframe * @param linewidth - Cell wireframe line width * @param dashSize - Cell wireframe dash size. Defaults to 0. Provide a value > @@ -615,7 +632,7 @@ function getAlignment(alignments, directions) { * @param gapSize - Cell wireframe dash size. Defaults to 0. Provide a value > * 0 for a dashed line. */ -function createCell(origin, basis, collapsed, periodicity, color, linewidth, dashSize, gapSize) { +function createCell(origin, basis, collapsed, color, linewidth, dashSize, gapSize) { const cell = new THREE.Object3D() let lineMaterial if (!(dashSize === 0 && gapSize === 0)) { @@ -632,42 +649,45 @@ function createCell(origin, basis, collapsed, periodicity, color, linewidth, das }) } - function addEdge(points, isDim, isCollapsed) { - if (!(isDim && isCollapsed)) { - const lineGeometry = new THREE.BufferGeometry().setFromPoints(points) - const line = new THREE.Line(lineGeometry, lineMaterial) - cell.add(line) - line.computeLineDistances() - } + function addEdge(points) { + const lineGeometry = new THREE.BufferGeometry().setFromPoints(points) + const line = new THREE.Line(lineGeometry, lineMaterial) + cell.add(line) + line.computeLineDistances() } for (let len = basis.length, i = 0; i < len; ++i) { + const isFirstCollapsed = collapsed[i] + if (isFirstCollapsed) continue const basisVector = basis[i] - const isCollapsed = collapsed[i] // First edge - const isDim1 = !periodicity[i] const points1 = [origin, basisVector.clone().add(origin)] - addEdge(points1, isDim1, isCollapsed) + addEdge(points1) // Second edge const secondIndex = (i + 1) % len + const isSecondCollapsed = collapsed[secondIndex] const secondAddition = basis[secondIndex].clone() - const isDim2 = !periodicity[i] || !periodicity[(i + 1) % 3] - const points2 = [secondAddition.clone().add(origin), basisVector.clone().add(secondAddition).add(origin)] - addEdge(points2, isDim2, isCollapsed) + if (!isSecondCollapsed) { + const points2 = [secondAddition.clone().add(origin), basisVector.clone().add(secondAddition).add(origin)] + addEdge(points2) + } // Third edge const thirdIndex = (i + 2) % len + const isThirdCollapsed = collapsed[thirdIndex] const thirdAddition = basis[thirdIndex].clone() - const isDim3 = !periodicity[i] || !periodicity[(i + 2) % 3] - const points3 = [thirdAddition.clone().add(origin), basisVector.clone().add(thirdAddition).add(origin)] - addEdge(points3, isDim3, isCollapsed) + if (!isThirdCollapsed) { + const points3 = [thirdAddition.clone().add(origin), basisVector.clone().add(thirdAddition).add(origin)] + addEdge(points3) + } // Fourth edge - const isDim4 = !periodicity[i] || !periodicity[(i + 2) % 3] || !periodicity[(i + 1) % 3] - const points4 = [secondAddition.clone().add(thirdAddition).add(origin), basisVector.clone().add(secondAddition).add(thirdAddition).add(origin)] - addEdge(points4, isDim4, isCollapsed) + if (!isSecondCollapsed && !isThirdCollapsed) { + const points4 = [secondAddition.clone().add(thirdAddition).add(origin), basisVector.clone().add(secondAddition).add(thirdAddition).add(origin)] + addEdge(points4) + } } return cell } @@ -1099,7 +1119,8 @@ function addObject3DToStage(object, stage) { function wrapPositions(component, wrap) { let posNew const basis = component.basis - const unitCell = component.object.unitcell + const cartToFrac = component.cartToFrac + const fracToCart = component.fracToCart const pbc = component.pbc if (!basis) return @@ -1120,10 +1141,10 @@ function wrapPositions(component, wrap) { posNew = component.posWrap } else { posNew = component.posCart.map(pos => { - return pos.clone().applyMatrix4(unitCell.cartToFrac) + return pos.clone().applyMatrix4(cartToFrac) }) const eps = 1e-2 - const epsArray = new THREE.Vector3(eps, eps, eps).applyMatrix4(unitCell.cartToFrac).toArray() + const epsArray = new THREE.Vector3(eps, eps, eps).applyMatrix4(cartToFrac).toArray() const center = [0.5, 0.5, 0.5] // Positions will be nearest to this location const shift = center.map((center, i) => pbc[i] ? center - 0.5 - epsArray[i] : 0) for (let len = posNew.length, i = 0; i < len; ++i) { @@ -1137,7 +1158,7 @@ function wrapPositions(component, wrap) { iFracPos.setComponent(i, remainder) } } - posNew = posNew.map(pos => pos.applyMatrix4(unitCell.fracToCart)) + posNew = posNew.map(pos => pos.applyMatrix4(fracToCart)) component.posWrap = posNew } // Use original cartesian positions diff --git a/gui/src/metainfo.json b/gui/src/metainfo.json index 5168d8421a44e770320e559c41d71b243212c9d6..7400400fe9f35fe64d2a7d3a6022da92b1800d0c 100644 --- a/gui/src/metainfo.json +++ b/gui/src/metainfo.json @@ -8554,6 +8554,25 @@ "type_data": "float64" }, "unit": "kilogram / meter ** 3" + }, + { + "m_def": "nomad.metainfo.metainfo.Quantity", + "m_parent_index": 9, + "m_parent_sub_section": "quantities", + "m_annotations": { + "elasticsearch": [ + "results.material.topology.cell.pbc" + ] + }, + "name": "pbc", + "description": "Periodic boundary conditions of the cell lattice vectors, given in order a, b, c.", + "type": { + "type_kind": "python", + "type_data": "bool" + }, + "shape": [ + 3 + ] } ] }, diff --git a/gui/src/searchQuantities.json b/gui/src/searchQuantities.json index 50f91d18eda746252f92a186963ea75bcf570a9d..94c5e467c92f9d1499fe0a0b6187b27a8d7ad3c3 100644 --- a/gui/src/searchQuantities.json +++ b/gui/src/searchQuantities.json @@ -1848,6 +1848,18 @@ "unit": "kilogram / meter ** 3", "aggregatable": false }, + "results.material.topology.cell.pbc": { + "name": "pbc", + "description": "Periodic boundary conditions of the cell lattice vectors, given in order a, b, c.", + "type": { + "type_kind": "python", + "type_data": "bool" + }, + "shape": [ + 3 + ], + "aggregatable": true + }, "results.material.topology.symmetry.bravais_lattice": { "name": "bravais_lattice", "description": "Identifier for the Bravais lattice in Pearson notation. The first lowercase letter\nidentifies the crystal family and can be one of the following: a (triclinic), b\n(monoclinic), o (orthorhombic), t (tetragonal), h (hexagonal) or c (cubic). The\nsecond uppercase letter identifies the centring and can be one of the following: P\n(primitive), S (face centred), I (body centred), R (rhombohedral centring) or F\n(all faces centred).", diff --git a/nomad/app/v1/routers/systems.py b/nomad/app/v1/routers/systems.py index 521e9e7086369be2155e63e5fe85af5502fac704..10158f9bdd6755d612d55259ee9e0a92a61287bb 100644 --- a/nomad/app/v1/routers/systems.py +++ b/nomad/app/v1/routers/systems.py @@ -24,9 +24,10 @@ from fastapi import APIRouter, Depends, Path, Query, status, HTTPException from fastapi.responses import StreamingResponse import ase.io import ase.build +from MDAnalysis.lib.util import NamedStream from nomad.utils import strip, deep_get, query_list_to_dict -from nomad.normalizing.common import ase_atoms_from_nomad_atoms +from nomad.normalizing.common import ase_atoms_from_nomad_atoms, mda_universe_from_nomad_atoms from nomad.datamodel.metainfo.simulation.system import Atoms as NOMADAtoms from .entries import answer_entry_archive_request, _bad_id_response @@ -42,7 +43,7 @@ format_list: List[dict] = [ 'description': 'Protein Data Bank file', 'extension': 'pdb', 'mime_type': 'chemical/x-pdb', - 'reader': 'ase', + 'reader': 'mdanalysis', } ] format_map: Dict[str, dict] = OrderedDict() @@ -54,6 +55,9 @@ format_description = "\n".join([f' - `{format["label"]}`: {format["description"] ase_format_map = { 'pdb': 'proteindatabank' } +mda_format_map = { + 'pdb': 'pdb' +} class TempFormatEnum(str, Enum): @@ -128,13 +132,26 @@ async def get_entry_raw_file( # Transform the system into ase atoms, catch any errors format_info = format_map[format] stringio = StringIO() + try: + atoms = NOMADAtoms.m_from_dict(result_dict) if format_info['reader'] == 'ase': - atoms = ase_atoms_from_nomad_atoms(NOMADAtoms.m_from_dict(result_dict)) + atoms = ase_atoms_from_nomad_atoms(atoms) ase_format = ase_format_map[format] ase.io.write(stringio, atoms, ase_format) elif format_info['reader'] == 'mdanalysis': - pass + universe = mda_universe_from_nomad_atoms(atoms) + mda_format = mda_format_map[format] + # For some reason NamedStream tries to close the stream: here we + # disable closing. The memory will be freed when stringio goes out + # of scope. + stringio.close = lambda: None # type: ignore[assignment] + namedstream = NamedStream(stringio, f'temp.{mda_format}', close=False, reset=False) + universe = mda_universe_from_nomad_atoms(atoms) + selection = universe.select_atoms("all") + selection.write(namedstream) + else: + raise ValueError("No reader for filetype.") except Exception as e: raise HTTPException( status_code=status.HTTP_500_INTERNAL_SERVER_ERROR, diff --git a/nomad/datamodel/results.py b/nomad/datamodel/results.py index a38792b18e35dcc3d1ead01c0dd1b762fe514b59..14ac370a9555fcc57a295359961e8d9f00f01cd2 100644 --- a/nomad/datamodel/results.py +++ b/nomad/datamodel/results.py @@ -245,14 +245,14 @@ class BandGap(MSection): ''' ) index = Quantity( - type=np.dtype(np.int32), + type=np.int32, description=''' Index of the data, e.g. spin channel index. ''', a_elasticsearch=Elasticsearch(material_entry_type), ) value = Quantity( - type=np.dtype(np.float64), + type=np.float64, shape=[], unit='joule', description=''' @@ -281,7 +281,7 @@ class BandGapElectronic(BandGap): ) energy_highest_occupied = Quantity( - type=np.dtype(np.float64), + type=np.float64, unit='joule', shape=[], description=''' @@ -289,7 +289,7 @@ class BandGapElectronic(BandGap): ''', ) energy_lowest_unoccupied = Quantity( - type=np.dtype(np.float64), + type=np.float64, unit='joule', shape=[], description=''' @@ -305,7 +305,7 @@ class LatticeParameters(MSection): ''', ) a = Quantity( - type=np.dtype(np.float64), + type=np.float64, unit='m', description=''' Length of the first basis vector. @@ -313,7 +313,7 @@ class LatticeParameters(MSection): a_elasticsearch=Elasticsearch(material_entry_type), ) b = Quantity( - type=np.dtype(np.float64), + type=np.float64, unit='m', description=''' Length of the second basis vector. @@ -321,7 +321,7 @@ class LatticeParameters(MSection): a_elasticsearch=Elasticsearch(material_entry_type), ) c = Quantity( - type=np.dtype(np.float64), + type=np.float64, unit='m', description=''' Length of the third basis vector. @@ -329,7 +329,7 @@ class LatticeParameters(MSection): a_elasticsearch=Elasticsearch(material_entry_type), ) alpha = Quantity( - type=np.dtype(np.float64), + type=np.float64, unit='radian', description=''' Angle between second and third basis vector. @@ -337,7 +337,7 @@ class LatticeParameters(MSection): a_elasticsearch=Elasticsearch(material_entry_type), ) beta = Quantity( - type=np.dtype(np.float64), + type=np.float64, unit='radian', description=''' Angle between first and third basis vector. @@ -345,7 +345,7 @@ class LatticeParameters(MSection): a_elasticsearch=Elasticsearch(material_entry_type), ) gamma = Quantity( - type=np.dtype(np.float64), + type=np.float64, unit='radian', description=''' Angle between first and second basis vector. @@ -381,19 +381,19 @@ class WyckoffSet(MSection): ''' ) x = Quantity( - type=np.dtype(np.float64), + type=np.float64, description=''' The free parameter x if present. ''' ) y = Quantity( - type=np.dtype(np.float64), + type=np.float64, description=''' The free parameter y if present. ''' ) z = Quantity( - type=np.dtype(np.float64), + type=np.float64, description=''' The free parameter z if present. ''' @@ -430,7 +430,7 @@ class Structure(MSection): a_elasticsearch=Elasticsearch(material_entry_type), ) lattice_vectors = Quantity( - type=np.dtype('float64'), + type=np.float64, shape=[3, 3], unit='m', description=''' @@ -438,7 +438,7 @@ class Structure(MSection): ''' ) cartesian_site_positions = Quantity( - type=np.dtype('float64'), + type=np.float64, shape=['n_sites', 3], unit='m', description=''' @@ -466,7 +466,7 @@ class Structure(MSection): ''' ) cell_volume = Quantity( - type=np.dtype(np.float64), + type=np.float64, unit='m ** 3', description=''' Volume of the cell. @@ -474,14 +474,14 @@ class Structure(MSection): a_elasticsearch=Elasticsearch(material_entry_type), ) atomic_density = Quantity( - type=np.dtype(np.float64), + type=np.float64, unit='1 / m ** 3', description=''' Atomic density of the material (atoms/volume).' ''' ) mass_density = Quantity( - type=np.dtype(np.float64), + type=np.float64, unit='kg / m ** 3', description=''' Mass density of the material. @@ -565,7 +565,7 @@ class Symmetry(MSection): ], ) hall_number = Quantity( - type=np.dtype(np.int32), + type=np.int32, shape=[], description=''' The Hall number for this system. @@ -595,7 +595,7 @@ class Symmetry(MSection): ], ) space_group_number = Quantity( - type=np.dtype(np.int32), + type=np.int32, shape=[], description=''' Specifies the International Union of Crystallography (IUC) number of the 3D space @@ -666,7 +666,7 @@ class Cell(MSection): ''', ) a = Quantity( - type=np.dtype(np.float64), + type=np.float64, unit='m', description=''' Length of the first basis vector. @@ -674,7 +674,7 @@ class Cell(MSection): a_elasticsearch=Elasticsearch(material_entry_type), ) b = Quantity( - type=np.dtype(np.float64), + type=np.float64, unit='m', description=''' Length of the second basis vector. @@ -682,7 +682,7 @@ class Cell(MSection): a_elasticsearch=Elasticsearch(material_entry_type), ) c = Quantity( - type=np.dtype(np.float64), + type=np.float64, unit='m', description=''' Length of the third basis vector. @@ -690,7 +690,7 @@ class Cell(MSection): a_elasticsearch=Elasticsearch(material_entry_type), ) alpha = Quantity( - type=np.dtype(np.float64), + type=np.float64, unit='radian', description=''' Angle between second and third basis vector. @@ -698,7 +698,7 @@ class Cell(MSection): a_elasticsearch=Elasticsearch(material_entry_type), ) beta = Quantity( - type=np.dtype(np.float64), + type=np.float64, unit='radian', description=''' Angle between first and third basis vector. @@ -706,7 +706,7 @@ class Cell(MSection): a_elasticsearch=Elasticsearch(material_entry_type), ) gamma = Quantity( - type=np.dtype(np.float64), + type=np.float64, unit='radian', description=''' Angle between first and second basis vector. @@ -714,7 +714,7 @@ class Cell(MSection): a_elasticsearch=Elasticsearch(material_entry_type), ) volume = Quantity( - type=np.dtype(np.float64), + type=np.float64, unit='m ** 3', description=''' Volume of the cell. @@ -722,7 +722,7 @@ class Cell(MSection): a_elasticsearch=Elasticsearch(material_entry_type), ) atomic_density = Quantity( - type=np.dtype(np.float64), + type=np.float64, unit='1 / m ** 3', description=''' Atomic density of the material (atoms/volume).' @@ -730,13 +730,21 @@ class Cell(MSection): a_elasticsearch=Elasticsearch(material_entry_type), ) mass_density = Quantity( - type=np.dtype(np.float64), + type=np.float64, unit='kg / m ** 3', description=''' Mass density of the material. ''', a_elasticsearch=Elasticsearch(material_entry_type), ) + pbc = Quantity( + type=bool, + shape=[3], + description=''' + Periodic boundary conditions of the cell lattice vectors, given in order a, b, c. + ''', + a_elasticsearch=Elasticsearch(material_entry_type), + ) class Prototype(MSection): @@ -828,7 +836,7 @@ class SymmetryNew(MSection): ], ) hall_number = Quantity( - type=np.dtype(np.int32), + type=np.int32, shape=[], description=''' The Hall number for this system. @@ -858,7 +866,7 @@ class SymmetryNew(MSection): ], ) space_group_number = Quantity( - type=np.dtype(np.int32), + type=np.int32, shape=[], description=''' Specifies the International Union of Crystallography (IUC) number of the 3D space @@ -907,7 +915,7 @@ class SymmetryNew(MSection): ''' ) origin_shift = Quantity( - type=np.dtype(np.float64), + type=np.float64, shape=[3], description=''' Vector $\\mathbf{p}$ from the origin of the standardized system to the origin of @@ -918,7 +926,7 @@ class SymmetryNew(MSection): ''' ) transformation_matrix = Quantity( - type=np.dtype(np.float64), + type=np.float64, shape=[3, 3], description=''' Matrix $\\mathbf{P}$ that is used to transform the standardized coordinates to the @@ -953,8 +961,8 @@ class Relation(MSection): class System(MSection): - '''Describes a physical structure as identified in an entry. Can be a part - of a larger structural hierarchy, i.e. the topology. + '''Describes a physical system as identified in an entry. Can be a part of a + larger structural hierarchy. ''' m_def = Section( description=''' @@ -1185,7 +1193,7 @@ class System(MSection): a_elasticsearch=Elasticsearch(material_type) ) indices = Quantity( - type=np.dtype(np.int64), + type=np.int64, shape=['*', '*'], description=''' Indices of the atoms belonging to this group. These indices refer to @@ -1466,7 +1474,7 @@ class DFT(MSection): ] ) exact_exchange_mixing_factor = Quantity( - type=np.dtype(np.float64), + type=np.float64, description='Amount of exact exchange mixed in with the XC functional (value range = [0,1]).', a_elasticsearch=Elasticsearch(material_entry_type) ) @@ -1605,7 +1613,7 @@ class Precision(MSection): ''' ) k_line_density = Quantity( - type=np.dtype(np.float64), + type=np.float64, shape=[], unit='m', description=''' @@ -1790,7 +1798,7 @@ class DOSElectronic(DOS): a_elasticsearch=Elasticsearch(material_entry_type, nested=True) ) energy_fermi = Quantity( - type=np.dtype(np.float64), + type=np.float64, unit='joule', shape=[], description=''' @@ -1866,7 +1874,7 @@ class BandStructureElectronic(BandStructure): a_elasticsearch=Elasticsearch(material_entry_type, nested=True) ) energy_fermi = Quantity( - type=np.dtype(np.float64), + type=np.float64, unit='joule', shape=[], description=''' @@ -2015,7 +2023,7 @@ class BulkModulus(MSection): ], ) value = Quantity( - type=np.dtype(np.float64), + type=np.float64, description='Bulk modulus value.', unit='pascal', a_elasticsearch=Elasticsearch(material_entry_type), @@ -2041,7 +2049,7 @@ class ShearModulus(MSection): ], ) value = Quantity( - type=np.dtype(np.float64), + type=np.float64, description='Shear modulus value.', unit='pascal', a_elasticsearch=Elasticsearch(material_entry_type), @@ -2126,7 +2134,7 @@ class QuantityDynamic(MSection): """ ) time = Quantity( - type=np.dtype(np.float64), + type=np.float64, shape=[], unit='second', description=""" @@ -2135,7 +2143,7 @@ class QuantityDynamic(MSection): """, ) time_step = Quantity( - type=np.dtype(np.float64), + type=np.float64, unit='second', description=""" The time step between successive evaluations. Provide either @@ -2143,7 +2151,7 @@ class QuantityDynamic(MSection): """, ) time_start = Quantity( - type=np.dtype(np.float64), + type=np.float64, unit='second', description=""" The time at which the evaluation started. Provide either this and @@ -2159,7 +2167,7 @@ class VolumeDynamic(QuantityDynamic): """ ) value = Quantity( - type=np.dtype(np.float64), + type=np.float64, shape=[], unit="m ** 3", description=""" @@ -2175,7 +2183,7 @@ class PressureDynamic(QuantityDynamic): """ ) value = Quantity( - type=np.dtype(np.float64), + type=np.float64, shape=[], unit="pascal", description=""" @@ -2191,7 +2199,7 @@ class TemperatureDynamic(QuantityDynamic): """ ) value = Quantity( - type=np.dtype(np.float64), + type=np.float64, shape=[], unit="kelvin", description=""" @@ -2207,7 +2215,7 @@ class EnergyDynamic(QuantityDynamic): """ ) value = Quantity( - type=np.dtype(np.float64), + type=np.float64, shape=[], unit="joule", description=""" @@ -2361,7 +2369,7 @@ class SolarCell(MSection): ''' ) efficiency = Quantity( - type=np.dtype(np.float64), + type=np.float64, shape=[], description=''' Power conversion effciency of a solar cell in percentage %. @@ -2369,7 +2377,7 @@ class SolarCell(MSection): a_elasticsearch=Elasticsearch(material_entry_type) ) fill_factor = Quantity( - type=np.dtype(np.float64), + type=np.float64, shape=[], description=''' Fill factor of a solar cell in absolute values (from 0 to 1). @@ -2377,7 +2385,7 @@ class SolarCell(MSection): a_elasticsearch=Elasticsearch(material_entry_type) ) open_circuit_voltage = Quantity( - type=np.dtype(np.float64), + type=np.float64, unit='V', shape=[], description=''' @@ -2386,7 +2394,7 @@ class SolarCell(MSection): a_elasticsearch=Elasticsearch(material_entry_type) ) short_circuit_current_density = Quantity( - type=np.dtype(np.float64), + type=np.float64, unit='A / m**2', shape=[], description=''' @@ -2395,7 +2403,7 @@ class SolarCell(MSection): a_elasticsearch=Elasticsearch(material_entry_type) ) illumination_intensity = Quantity( - type=np.dtype(np.float64), + type=np.float64, unit=('W/m**2'), shape=[], description=''' @@ -2404,7 +2412,7 @@ class SolarCell(MSection): a_elasticsearch=Elasticsearch(material_entry_type) ) device_area = Quantity( - type=np.dtype(np.float64), + type=np.float64, unit=('m**2'), shape=[], description=''' diff --git a/nomad/normalizing/common.py b/nomad/normalizing/common.py index 26785c6acecd60f39be96d51f594a88905fe1f2f..ffc60fcc1f29f841d191a056a326509ad5529cbd 100644 --- a/nomad/normalizing/common.py +++ b/nomad/normalizing/common.py @@ -19,6 +19,7 @@ import numpy as np from ase import Atoms from typing import List, Set, Any, Optional from nptyping import NDArray +import MDAnalysis as mda from matid import SymmetryAnalyzer # pylint: disable=import-error from matid.symmetry.wyckoffset import WyckoffSet as WyckoffSetMatID # pylint: disable=import-error import matid.geometry # pylint: disable=import-error @@ -146,13 +147,12 @@ def cell_from_ase_atoms(atoms: Atoms) -> Cell: cell.gamma = None if np.isnan(gamma) else gamma * ureg.radian volume = atoms.cell.volume * ureg.angstrom ** 3 + cell.volume = volume mass = atomutils.get_summed_atomic_mass(atoms.get_atomic_numbers()) * ureg.kg - mass_density = mass / volume + cell.mass_density = None if volume == 0 else mass / volume number_of_atoms = atoms.get_number_of_atoms() - atomic_density = number_of_atoms / volume - cell.volume = volume - cell.atomic_density = atomic_density - cell.mass_density = mass_density + cell.atomic_density = None if volume == 0 else number_of_atoms / volume + cell.pbc = np.zeros(3, bool)[:] = atoms.get_pbc() return cell @@ -288,47 +288,91 @@ def ase_atoms_from_structure(system: Structure) -> Atoms: ) +def mda_universe_from_nomad_atoms(system: Atoms, logger=None) -> mda.Universe: + '''Returns an instance of mda.Universe from a NOMAD Atoms-section. + + Args: + system: The atoms to transform + + Returns: + A new mda.Universe created from the given data. + ''' + n_atoms = len(system.positions) + n_residues = 1 + atom_resindex = [0] * n_atoms + residue_segindex = [0] + + universe = mda.Universe.empty( + n_atoms, + n_residues=n_residues, + atom_resindex=atom_resindex, + residue_segindex=residue_segindex, + trajectory=True + ) + + # Add positions + universe.atoms.positions = system.positions.to(ureg.angstrom).magnitude + + # Add atom attributes + atom_names = system.labels + universe.add_TopologyAttr('name', atom_names) + universe.add_TopologyAttr('type', atom_names) + universe.add_TopologyAttr('element', atom_names) + + # Add the box dimensions + universe.atoms.dimensions = atomutils.cell_to_cellpar( + system.lattice_vectors.to(ureg.angstrom).magnitude, + degrees=True + ) + + return universe + + def structures_2d(original_atoms, logger=None): conv_atoms = None prim_atoms = None wyckoff_sets = None spg_number = None + symm_system = original_atoms.copy() + n_pbc = sum(original_atoms.get_pbc()) + try: - # Get dimension of system by also taking into account the covalent radii - dimensions = matid.geometry.get_dimensions(original_atoms, [True, True, True]) - basis_dimensions = np.linalg.norm(original_atoms.get_cell(), axis=1) - gaps = basis_dimensions - dimensions - periodicity = gaps <= config.normalize.cluster_threshold - - # If two axis are not periodic, return. This only happens if the vacuum - # gap is not aligned with a cell vector or if the linear gap search is - # unsufficient (the structure is "wavy" making also the gap highly - # nonlinear). - if sum(periodicity) != 2: - # TODO: @ Lauri: I'm not sure if this is smart/wise, but the periodicity calculation doesn't seem to work for conventional cell systems. The calculated periodicity is always (True, True, True). At least that's my interpretation of what's happening. This is my quick and dirty and at the moment only solution for this problem. - periodicity = original_atoms.pbc - # if logger: - # logger.error("could not detect the periodic dimensions in a 2D system") - # return conv_atoms, prim_atoms, wyckoff_sets, spg_number - - # Center the system in the non-periodic direction, also taking - # periodicity into account. The get_center_of_mass()-function in MatID - # takes into account periodicity and can produce the correct CM unlike - # the similar function in ASE. - pbc_cm = matid.geometry.get_center_of_mass(original_atoms) - cell_center = 0.5 * np.sum(original_atoms.get_cell(), axis=0) - translation = cell_center - pbc_cm - symm_system = original_atoms.copy() - symm_system.translate(translation) - symm_system.wrap() - - # Set the periodicity according to detected periodicity in order - # for SymmetryAnalyzer to use the symmetry analysis designed for 2D - # systems. - symm_system.set_pbc(periodicity) + # If the given system if fully periodic, try to extract a 2D system by + # checking the presence of vacuum gaps. + if n_pbc == 3: + # Get dimension of system by also taking into account the covalent radii + dimensions = matid.geometry.get_dimensions(original_atoms, [True, True, True]) + basis_dimensions = np.linalg.norm(original_atoms.get_cell(), axis=1) + gaps = basis_dimensions - dimensions + periodicity = gaps <= config.normalize.cluster_threshold + + # If two axis are not periodic, return. This only happens if the vacuum + # gap is not aligned with a cell vector or if the linear gap search is + # unsufficient (the structure is "wavy" making also the gap highly + # nonlinear). + if sum(periodicity) != 2: + if logger: + logger.error("could not detect the periodic dimensions in a 2D system") + return conv_atoms, prim_atoms, wyckoff_sets, spg_number + + # Center the system in the non-periodic direction, also taking + # periodicity into account. The get_center_of_mass()-function in MatID + # takes into account periodicity and can produce the correct CM unlike + # the similar function in ASE. + pbc_cm = matid.geometry.get_center_of_mass(original_atoms) + cell_center = 0.5 * np.sum(original_atoms.get_cell(), axis=0) + translation = cell_center - pbc_cm + symm_system = original_atoms.copy() + symm_system.translate(translation) + symm_system.wrap() + + # Set the periodicity according to detected periodicity in order + # for SymmetryAnalyzer to use the symmetry analysis designed for 2D + # systems. + symm_system.set_pbc(periodicity) symmetry_analyzer = SymmetryAnalyzer( symm_system, - config.normalize.symmetry_tolerance, + 0.4, # The value is increased here to better match 2D materials. config.normalize.flat_dim_threshold ) @@ -341,17 +385,6 @@ def structures_2d(original_atoms, logger=None): # dimensions. conv_atoms = atomutils.get_minimized_structure(conv_atoms) prim_atoms = atomutils.get_minimized_structure(prim_atoms) - - # Swap the cell axes so that the non-periodic one is always the - # last basis (=c) - swap_dim = 2 - for i, periodic in enumerate(conv_atoms.get_pbc()): - if not periodic: - non_periodic_dim = i - break - if non_periodic_dim != swap_dim: - atomutils.swap_basis(conv_atoms, non_periodic_dim, swap_dim) - atomutils.swap_basis(prim_atoms, non_periodic_dim, swap_dim) except Exception as e: if logger: logger.error( diff --git a/nomad/normalizing/topology.py b/nomad/normalizing/topology.py index 65d00f861f4faf72bb53c7a9a22c33865ba215fd..4fd7ddead8c5e9d2503e171e2ca298c8b4b87918 100644 --- a/nomad/normalizing/topology.py +++ b/nomad/normalizing/topology.py @@ -16,21 +16,16 @@ # limitations under the License. # -from typing import Dict, List, Tuple, Optional +from typing import Dict, List, Optional from collections import defaultdict from ase import Atoms from ase.data import chemical_symbols -from ase.geometry.cell import complete_cell import numpy as np -import matid.geometry # pylint: disable=import-error -from matid.classification.structureclusterer import StructureClusterer # pylint: disable=import-error -from matid import Classifier # pylint: disable=import-error -from matid.classifications import Class0D, Atom, Class1D, Material2D, Surface, Class3D, Class2D, Unknown # pylint: disable=import-error +from matid.clustering import Clusterer, Cluster, Classification # pylint: disable=import-error from matid.symmetry.symmetryanalyzer import SymmetryAnalyzer # pylint: disable=import-error from nomad import atomutils -from nomad.units import ureg from nomad.atomutils import Formula from nomad.datamodel.results import Symmetry, Material, System, Relation, Prototype from nomad.datamodel.metainfo.simulation.system import Atoms as NOMADAtoms @@ -65,11 +60,6 @@ def get_topology_original(material: Material, atoms: NOMADAtoms = None) -> Syste method='parser', label='original', description='A representative system chosen from the original simulation.', - material_id=material.material_id, - material_name=material.material_name, - structural_type=material.structural_type, - functional_type=material.functional_type, - compound_type=material.compound_type, chemical_formula_hill=material.chemical_formula_hill, chemical_formula_iupac=material.chemical_formula_iupac, chemical_formula_anonymous=material.chemical_formula_anonymous, @@ -250,11 +240,13 @@ class TopologyNormalizer(): atoms = ase_atoms_from_nomad_atoms(nomad_atoms) except Exception: return None - # TODO: MatID does not currently support non-periodic structures - if not any(atoms.pbc): + # TODO: MatID does not currently support non-periodic structures or + # structures with zero-sized cell + cell = atoms.get_cell() + if not any(atoms.pbc) or cell is None or cell.volume == 0: return None - # TODO: Currently the topology creation is skipped completely. Needs to - # be re-enabled once we work out performance problems. + # TODO: The matid processing is currently completely skipped. Should be + # enabled after more testing. n_atoms = len(atoms) if n_atoms > 0: return None @@ -266,22 +258,25 @@ class TopologyNormalizer(): # Add all meaningful clusters to the topology clusters = self._perform_matid_clustering(atoms) - cluster_indices_list, cluster_symmetries = self._filter_clusters(clusters) - for indices, symm in zip(cluster_indices_list, cluster_symmetries): - subsystem_atoms = atoms[indices] - subsystem = self._create_subsystem(subsystem_atoms, indices) + filtered_clusters = self._filter_clusters(clusters) + for cluster in filtered_clusters: + subsystem = self._create_subsystem(cluster) structural_type = subsystem.structural_type - # TODO: Currently only 2D and surface subsystems are included + # Currently only 2D and surface subsystems are included if structural_type in {'surface', '2D'}: add_system(subsystem, topology, original) add_system_info(subsystem, topology) - conventional_cell = self._create_conv_cell_system(symm, structural_type) - add_system(conventional_cell, topology, subsystem) - add_system_info(conventional_cell, topology) + try: + conventional_cell = self._create_conv_cell_system(cluster, structural_type) + except Exception as e: + raise ValueError("conventional cell infomation could not be created") from e + else: + add_system(conventional_cell, topology, subsystem) + add_system_info(conventional_cell, topology) return list(topology.values()) - def _create_subsystem(self, atoms: Atoms, indices: List[int]) -> System: + def _create_subsystem(self, cluster: Cluster) -> System: ''' Creates a new subsystem as detected by MatID. ''' @@ -290,13 +285,31 @@ class TopologyNormalizer(): label='subsystem', description='Automatically detected subsystem.', system_relation=Relation(type='subsystem'), - indices=[indices] + indices=[list(cluster.indices)] ) - subsystem.structural_type = self._system_type_analysis(atoms) + + classification = 'unavailable' + try: + classification = cluster.classification() + except Exception as e: + self.logger.error( + 'matid project system classification failed', exc_info=e, error=str(e) + ) + type_map = { + Classification.Class3D: 'bulk', + Classification.Atom: 'atom', + Classification.Class0D: 'molecule / cluster', + Classification.Class1D: '1D', + Classification.Class2D: '2D', + Classification.Surface: 'surface', + Classification.Material2D: '2D', + Classification.Unknown: 'unavailable' + } + subsystem.structural_type = type_map.get(classification, 'unavailable') return subsystem - def _create_conv_cell_system(self, symm, structural_type: str): + def _create_conv_cell_system(self, cluster: Cluster, structural_type: str): ''' Creates a new topology item for a conventional cell. ''' @@ -305,27 +318,18 @@ class TopologyNormalizer(): label='conventional cell', system_relation=Relation(type='subsystem'), ) - conv_system = symm.get_conventional_system() if structural_type == 'surface': - symmsystem.description = 'The conventional cell of the bulk material from which the surface is constructed from.' - symmsystem.structural_type = 'bulk' - symmsystem.atoms = nomad_atoms_from_ase_atoms(conv_system) + self._add_conventional_surface(cluster, symmsystem) elif structural_type == '2D': - symmsystem.description = 'The conventional cell of the 2D material.' - symmsystem.structural_type = '2D' - symmsystem.atoms = nomad_atoms_from_ase_atoms(conv_system) + self._add_conventional_2d(cluster, symmsystem) - if structural_type == 'surface': - symmsystem = self._create_symmsystem_surface(symm, symmsystem) - elif structural_type == '2D': - symmsystem = self._create_symmsystem_2D(symm, symmsystem) return symmsystem def _perform_matid_clustering(self, atoms: Atoms) -> List: ''' Performs the clustering with MatID ''' - clusterer = StructureClusterer() + clusterer = Clusterer() clusters = clusterer.get_clusters( atoms, max_cell_size=6, @@ -336,10 +340,9 @@ class TopologyNormalizer(): ) return clusters - def _filter_clusters(self, clusters: StructureClusterer) -> Tuple[List[List[int]], List[Symmetry]]: + def _filter_clusters(self, clusters: List[Cluster]) -> List[Cluster]: ''' - Filters all clusters < 2 atoms and creates a cluster indices list of the remaining - clusters + Filter out clusters that do not have a meaningful representation in NOMAD. ''' # TODO: Add proper decision mechanism for identifying when the detected # cluster are OK. For now, we filter out cluster that have very small @@ -350,63 +353,18 @@ class TopologyNormalizer(): # only on the number of atoms in the cluster. Maybe the cluster given # out by MatID could already be grouped a bit better, e.g. all outliers # would be grouped into one cluster? - filtered_cluster = filter(lambda x: len(x.indices) > 1, clusters) - cluster_indices_list: List[List[int]] = [] - cluster_symmetries: List[Symmetry] = [] - for cluster in filtered_cluster: - indices = list(cluster.indices) - cluster_indices_list += [indices] - regions = cluster.regions - number_of_atoms: List[int] = [] - for region in regions: - if region: - number_of_atoms.append(region.cell.get_number_of_atoms()) - else: - number_of_atoms.append(-1) - - # If there are 2 regions with the same size, the one with the smaller index is selected - largest_region_index = number_of_atoms.index(max(number_of_atoms)) - largest_region_system = regions[largest_region_index].cell - - # TODO: only SymmetryAnalyzer for 2D and surface - symm = SymmetryAnalyzer(largest_region_system) - cluster_symmetries += [symm] - return cluster_indices_list, cluster_symmetries + return list(filter(lambda x: len(x.indices) > 1, clusters)) - def _system_type_analysis(self, atoms: Atoms) -> matid.classifications: - ''' - Classifies ase.atoms and returns the MatID system type as a string. - ''' - system_types = {Class3D: 'bulk', - Atom: 'atom', - Class0D: 'molecule / cluster', - Class1D: '1D', - Class2D: 'unavailable', - Surface: 'surface', - Material2D: '2D', - Unknown: 'unavailable'} - try: - classifier = Classifier(radii='covalent', cluster_threshold=1) - cls = classifier.classify(atoms) - except Exception as e: - self.logger.error( - 'matid project system classification failed', exc_info=e, error=str(e)) - return 'unavailable' - else: - classification = type(cls) - try: - system_type = system_types[classification] - except Exception as e: - self.logger.error( - 'matid project system classification unavailable', exc_info=e, error=str(e)) - system_type = 'unavailable' - return system_type - - def _create_symmsystem_surface(self, symm: SymmetryAnalyzer, subsystem: System) -> System: + def _add_conventional_surface(self, cluster: Cluster, subsystem: System) -> None: ''' Creates the subsystem with the symmetry information of the conventional cell ''' + cell = cluster.cell() + symm = SymmetryAnalyzer(cell) + subsystem.description = 'The conventional cell of the bulk material from which the surface is constructed from.' + subsystem.structural_type = 'bulk' conv_system = symm.get_conventional_system() + subsystem.atoms = nomad_atoms_from_ase_atoms(conv_system) prototype = self._create_prototype(symm, conv_system) spg_number = symm.get_space_group_number() subsystem.prototype = prototype @@ -416,36 +374,35 @@ class TopologyNormalizer(): material_id = material_id_bulk(spg_number, wyckoff_sets) subsystem.material_id = material_id subsystem.symmetry = symmetry - return subsystem - def _create_symmsystem_2D(self, symm: SymmetryAnalyzer, subsystem: System) -> System: + def _add_conventional_2d(self, cluster: Cluster, subsystem: System) -> None: ''' Creates the subsystem with the symmetry information of the conventional cell ''' - subsystem_atoms = Atoms( - symbols=subsystem.atoms.labels, - positions=subsystem.atoms.positions.to(ureg.angstrom), - cell=complete_cell(subsystem.atoms.lattice_vectors.to(ureg.angstrom)), - pbc=np.array(subsystem.atoms.periodic, dtype=bool) - ) - conv_atoms, __, wyckoff_sets, spg_number = structures_2d(subsystem_atoms) + cell = cluster.cell() + subsystem.description = 'The conventional cell of the 2D material.' + subsystem.structural_type = '2D' + conv_atoms, _, wyckoff_sets, spg_number = structures_2d(cell) + subsystem.cell = cell_from_ase_atoms(conv_atoms) + subsystem.atoms = nomad_atoms_from_ase_atoms(conv_atoms) # Here we zero out the irrelevant lattice parameters to correctly handle # 2D systems with nonzero thickness (e.g. MoS2). - if subsystem.cell.c: - subsystem.cell.c = 0.0 - if subsystem.cell.alpha: + if subsystem.cell.c is not None: + subsystem.cell.c = None + if subsystem.cell.alpha is not None: subsystem.cell.alpha = None - if subsystem.cell.beta: + if subsystem.cell.beta is not None: subsystem.cell.beta = None + if subsystem.cell.atomic_density is not None: + subsystem.cell.atomic_density = None + if subsystem.cell.mass_density is not None: + subsystem.cell.mass_density = None + if subsystem.cell.volume is not None: + subsystem.cell.volume = None - prototype = self._create_prototype(symm, conv_atoms) - subsystem.prototype = prototype subsystem.material_id = material_id_2d(spg_number, wyckoff_sets) - symmetry = None - subsystem.symmetry = symmetry - return subsystem def _create_symmetry(self, symm: SymmetryAnalyzer) -> Symmetry: international_short = symm.get_space_group_international_short() diff --git a/tests/normalizing/conftest.py b/tests/normalizing/conftest.py index 5d4fcef1312caa23c5a1c8e4aed5258e0b0a77d0..99178e70016c066b2c164026e5c200cc289ef300 100644 --- a/tests/normalizing/conftest.py +++ b/tests/normalizing/conftest.py @@ -37,6 +37,7 @@ from nomad.datamodel.results import ( System as ResultSystem ) from nomad.datamodel.optimade import Species +from nomad.normalizing.common import cell_from_ase_atoms, nomad_atoms_from_ase_atoms from nomad.datamodel.metainfo.simulation.run import Run, Program from nomad.datamodel.metainfo.simulation.method import ( Method, BasisSet, Electronic, DFT, XCFunctional, Functional, @@ -1215,16 +1216,33 @@ def crystal_structure_properties(crystal_structure): return crystal_structure_property -def single_Cu_surface_atoms(): - # create simple Cu surface - Cu_fcc_100 = ase.build.fcc100('Cu', (3, 5, 5), vacuum=10, periodic=True) - Cu_fcc_100.rattle(stdev=0.001, seed=None, rng=None) - Cu_fcc_110 = ase.build.fcc110('Cu', (3, 5, 5), vacuum=10, periodic=True) - Cu_fcc_110.rattle(stdev=0.001, seed=None, rng=None) - return Cu_fcc_100, Cu_fcc_110 +def conv_fcc(symbols): + return ase.build.bulk(symbols, crystalstructure='fcc', a=3.61, cubic=True) -def single_Cu_surface_prototype(): +def conv_bcc(symbols): + return ase.build.bulk(symbols, crystalstructure='bcc', cubic=True) + + +def surf(conv_cell, indices, layers=[3, 3, 2], vacuum=10): + surface = ase.build.surface(conv_cell, indices, layers[2], vacuum=vacuum, periodic=True) + surface *= [layers[0], layers[1], 1] + return surface + + +def stack(a, b): + stacked = ase.build.stack(a, b, axis=2, distance=3, maxstrain=6.7) + stacked.set_pbc([True, True, False]) + ase.build.add_vacuum(stacked, 10) + return stacked + + +def rattle(atoms): + atoms.rattle(stdev=0.001, seed=7, rng=None) + return atoms + + +def single_cu_surface_prototype(): prototype_Cu_fcc = Prototype() prototype_Cu_fcc.aflow_id = "A_cF4_225_a" prototype_Cu_fcc.assignment_method = "normalized-wyckoff" @@ -1233,30 +1251,24 @@ def single_Cu_surface_prototype(): return prototype_Cu_fcc -def single_Cu_surface_topology() -> List[ResultSystem]: +def single_cu_surface_topology() -> List[ResultSystem]: '''Copper surface topology''' - prototype_Cu_fcc = single_Cu_surface_prototype() - cell = Cell( - a=3.610000000000001 * ureg.angstrom, - b=3.610000000000001 * ureg.angstrom, - c=3.610000000000001 * ureg.angstrom, - alpha=1.5707963267948966 * ureg.rad, - beta=1.5707963267948966 * ureg.rad, - gamma=1.5707963267948966 * ureg.rad - ) + prototype_Cu_fcc = single_cu_surface_prototype() + conv_cell = conv_fcc('Cu') + surface = surf(conv_cell, (1, 0, 0)) - # create Cu topology system label_Cu = 'subsystem' + n_atoms = len(surface) structural_type_Cu = 'surface' elements_Cu = ['Cu'] - formula_hill_Cu = 'Cu75' - formula_reduced_Cu = 'Cu75' - formula_anonymous_Cu = 'A75' + formula_hill_Cu = f'Cu{n_atoms}' + formula_reduced_Cu = f'Cu{n_atoms}' + formula_anonymous_Cu = f'A{n_atoms}' system_relation = Relation() system_relation.type = "subsystem" - indices_Cu = [i for i in range(75)] + indices = list(range(n_atoms)) - subsystem_Cu = create_system(label_Cu, structural_type_Cu, elements_Cu, formula_hill_Cu, formula_reduced_Cu, formula_anonymous_Cu, system_relation, indices=indices_Cu) + subsystem_Cu = create_system(label_Cu, structural_type_Cu, elements_Cu, formula_hill_Cu, formula_reduced_Cu, formula_anonymous_Cu, system_relation, indices=indices) topologies_Cu = [subsystem_Cu] label_Cu_conv = 'conventional cell' @@ -1266,11 +1278,6 @@ def single_Cu_surface_topology() -> List[ResultSystem]: formula_anonymous_Cu_conv = 'A4' material_id_Cu_conv = "3M6onRRrQbutydx916-Y15I79Z_X" - atoms_Cu_conv = NOMADAtoms() - atoms_Cu_conv.periodic = [False, False, False] - atoms_Cu_conv.lattice_vectors = [[3.609999999999999, 0.0, 0.0], [0.0, 3.609999999999999e-10, 0.0], [0.0, 0.0, 3.609999999999999]] * ureg.angstrom - atoms_Cu_conv.positions = [[0.0, 0.0, 0.0], [0.0, 1.8049999999999995, 1.8049999999999995], [1.8049999999999995, 0.0, 1.8049999999999995], [1.8049999999999995, 1.8049999999999995, 0.0]] * ureg.angstrom - atoms_Cu_conv.labels = ["Cu", "Cu", "Cu", "Cu"] species = Species() species.name = "Cu" species.chemical_symbols = ["Cu"] @@ -1280,22 +1287,15 @@ def single_Cu_surface_topology() -> List[ResultSystem]: wyckoff_sets.indices = [0, 1, 2, 3] wyckoff_sets.element = "Cu" + cell = cell_from_ase_atoms(conv_cell) + atoms = nomad_atoms_from_ase_atoms(conv_cell) Symmetry_fcc = crystal_structure_properties('fcc') - convsystem_Cu = create_system(label_Cu_conv, structural_type_Cu_conv, elements_Cu, formula_hill_Cu_conv, formula_reduced_Cu_conv, formula_anonymous_Cu_conv, system_relation, material_id=material_id_Cu_conv, atoms=atoms_Cu_conv, cell=cell, symmetry=Symmetry_fcc, prototype=prototype_Cu_fcc) + convsystem_Cu = create_system(label_Cu_conv, structural_type_Cu_conv, elements_Cu, formula_hill_Cu_conv, formula_reduced_Cu_conv, formula_anonymous_Cu_conv, system_relation, material_id=material_id_Cu_conv, atoms=atoms, cell=cell, symmetry=Symmetry_fcc, prototype=prototype_Cu_fcc) topologies_Cu.append(convsystem_Cu) return topologies_Cu -def single_Cr_surface_atoms() -> Atoms: - '''Cr surface system''' - Cr_bcc_100 = ase.build.bcc100('Cr', (5, 5, 5), vacuum=10, periodic=True) - Cr_bcc_100.rattle(stdev=0.001, seed=None, rng=None) - Cr_bcc_110 = ase.build.bcc110('Cr', (5, 5, 5), vacuum=10, periodic=True) - Cr_bcc_110.rattle(stdev=0.001, seed=None, rng=None) - return Cr_bcc_100, Cr_bcc_110 - - -def single_Cr_surface_topology() -> List[ResultSystem]: +def single_cr_surface_topology() -> List[ResultSystem]: '''Cr surface topology''' prototype_Cr_bcc = Prototype() prototype_Cr_bcc.aflow_id = "A_cI2_229_a" @@ -1303,15 +1303,18 @@ def single_Cr_surface_topology() -> List[ResultSystem]: prototype_Cr_bcc.label = "229-W-cI2" prototype_Cr_bcc.formula = "W2" + conv_cell = conv_bcc('Cr') + surface = surf(conv_cell, (1, 0, 0)) + n_atoms = len(surface) label = 'subsystem' structural_type = 'surface' elements = ['Cr'] - formula_hill = 'Cr125' - formula_reduced = 'Cr125' - formula_anonymous = 'A125' + formula_hill = f'Cr{n_atoms}' + formula_reduced = f'Cr{n_atoms}' + formula_anonymous = f'A{n_atoms}' system_relation = Relation() system_relation.type = "subsystem" - indices = [i for i in range(75)] + indices = list(range(n_atoms)) subsystem_Cr = create_system(label, structural_type, elements, formula_hill, formula_reduced, formula_anonymous, system_relation, indices=indices) topologies_Cr = [subsystem_Cr] @@ -1323,11 +1326,6 @@ def single_Cr_surface_topology() -> List[ResultSystem]: formula_reduced = 'Cr2' formula_anonymous = 'A2' material_id = "MDlo8h4C2Ppy-kLY9fHRovgnTN9T" - atoms = NOMADAtoms() - atoms.periodic = [False, False, False] - atoms.lattice_vectors = [[2.8800000000000014, 0.0, 0.0], [0.0, 2.8800000000000014e-10, 0.0], [0.0, 0.0, 2.8800000000000014]] * ureg.angstrom - atoms.positions = [[0.0, 0.0, 0.0], [1.4400000000000007, 1.4400000000000007, 1.4400000000000007]] * ureg.angstrom - atoms.labels = ["Cr", "Cr"] species = Species() species.name = "Cr" species.chemical_symbols = ["Cr"] @@ -1336,31 +1334,30 @@ def single_Cr_surface_topology() -> List[ResultSystem]: wyckoff_sets.wyckoff_letter = "a" wyckoff_sets.indices = [0, 1] wyckoff_sets.element = "Cr" - cell = Cell( - a=2.8800000000000014 * ureg.angstrom, - b=2.8800000000000014 * ureg.angstrom, - c=2.8800000000000014 * ureg.angstrom, - alpha=1.5707963267948966 * ureg.rad, - beta=1.5707963267948966 * ureg.rad, - gamma=1.5707963267948966 * ureg.rad, - ) + + atoms = nomad_atoms_from_ase_atoms(conv_cell) + cell = cell_from_ase_atoms(conv_cell) Symmetry_bcc = crystal_structure_properties('bcc') convsystem_Cr = create_system(label, structural_type, elements, formula_hill, formula_reduced, formula_anonymous, system_relation, material_id=material_id, atoms=atoms, cell=cell, symmetry=Symmetry_bcc, prototype=prototype_Cr_bcc) topologies_Cr.append(convsystem_Cr) return topologies_Cr -def single_Ni_surface_topology() -> List[ResultSystem]: +def single_ni_surface_topology() -> List[ResultSystem]: '''Ni surface topology''' + conv_cell = conv_fcc('Ni') + surface = surf(conv_cell, (1, 0, 0)) + n_atoms = len(surface) + label = 'subsystem' structural_type = 'surface' elements = ['Ni'] - formula_hill = 'Ni75' - formula_reduced = 'Ni75' - formula_anonymous = 'A75' + formula_hill = f'Ni{n_atoms}' + formula_reduced = f'Ni{n_atoms}' + formula_anonymous = f'A{n_atoms}' system_relation = Relation() system_relation.type = "subsystem" - indices = [i for i in range(75, 150)] + indices = list(range(n_atoms)) subsystem_Ni = create_system(label, structural_type, elements, formula_hill, formula_reduced, formula_anonymous, system_relation, indices=indices) topologies_Ni = [subsystem_Ni] @@ -1372,20 +1369,6 @@ def single_Ni_surface_topology() -> List[ResultSystem]: formula_reduced = 'Ni4' formula_anonymous = 'A4' material_id = "NdIWxnQzlp-aeP1IM2d8YJ04h6T0" - atoms = NOMADAtoms() - atoms.periodic = [False, False, False] - atoms.lattice_vectors = [ - [3.52, 0.0, 0.0], - [0.0, 3.52, 0.0], - [0.0, 0.0, 3.52] - ] * ureg.angstrom - atoms.positions = [ - [0.0, 0.0, 0.0], - [0.0, 1.76, 1.76], - [1.76, 0.0, 1.76], - [1.76, 1.76, 0.0] - ] * ureg.angstrom - atoms.labels = ["Ni", "Ni", "Ni", "Ni"] species = Species() species.name = "Ni" species.chemical_symbols = ["Ni"] @@ -1394,38 +1377,30 @@ def single_Ni_surface_topology() -> List[ResultSystem]: wyckoff_sets.wyckoff_letter = "a" wyckoff_sets.indices = [0, 1, 2, 3] wyckoff_sets.element = "Ni" - cell = Cell( - a=3.52 * ureg.angstrom, - b=3.52 * ureg.angstrom, - c=3.52 * ureg.angstrom, - alpha=1.5707963267948966 * ureg.rad, - beta=1.5707963267948966 * ureg.rad, - gamma=1.5707963267948966 * ureg.rad - ) + cell = cell_from_ase_atoms(conv_cell) + atoms = nomad_atoms_from_ase_atoms(conv_cell) Symmetry_fcc = crystal_structure_properties('fcc') - prototype_Cu_fcc = single_Cu_surface_prototype() + prototype_Cu_fcc = single_cu_surface_prototype() convsystem_Ni = create_system(label, structural_type, elements, formula_hill, formula_reduced, formula_anonymous, system_relation, material_id=material_id, atoms=atoms, cell=cell, symmetry=Symmetry_fcc, prototype=prototype_Cu_fcc) topologies_Ni.append(convsystem_Ni) return topologies_Ni -def stacked_Cu_Ni_surface() -> Atoms: - '''Stacked Cu and Ni surfaces''' - Cu_fcc_111 = ase.build.fcc111('Cu', (3, 5, 5), vacuum=None, periodic=False) - topologies_Cu = single_Cu_surface_topology() - - Ni_fcc_111 = ase.build.fcc111('Ni', (3, 5, 5), vacuum=None, periodic=False) - topologies_Ni = single_Ni_surface_topology() +def stacked_cu_ni_surface_topology() -> List[ResultSystem]: + topologies_cu = single_cu_surface_topology() + topologies_ni = single_ni_surface_topology() - CuNi_fcc_111 = ase.build.stack(Cu_fcc_111, Ni_fcc_111, axis=2, distance=2, maxstrain=2.4) - CuNi_fcc_111.rattle(stdev=0.001, seed=None, rng=None) + # Indices are modified + n_atoms_cu = len(topologies_cu[0].indices) + n_atoms_ni = len(topologies_ni[0].indices) + topologies_cu[0].indices = list(range(0, n_atoms_cu)) + topologies_ni[0].indices = list(range(n_atoms_cu, n_atoms_cu + n_atoms_ni)) - # create stacked Cu and Ni surface topology - topologies_Cu_Ni = topologies_Cu + topologies_Ni - return CuNi_fcc_111, topologies_Cu_Ni + topologies_cu_ni = topologies_cu + topologies_ni + return topologies_cu_ni -def single_2D_graphene_layer_atoms() -> Atoms: +def graphene() -> Atoms: '''Graphene system''' symbols_C = ['C', 'C'] positions_C = [ @@ -1443,11 +1418,10 @@ def single_2D_graphene_layer_atoms() -> Atoms: cell=cell_C, pbc=True ) * [4, 4, 1] - system_C.rattle(stdev=0.001, seed=None, rng=None) return system_C -def single_2D_graphene_layer_topology() -> List[ResultSystem]: +def graphene_topology() -> List[ResultSystem]: '''Graphene topology''' label_C = 'subsystem' structural_type_C = '2D' @@ -1486,17 +1460,17 @@ def single_2D_graphene_layer_topology() -> List[ResultSystem]: species.chemical_symbols = ["C"] species.concentration = [1.0] cell = Cell( - a=2.4677241413662866 * ureg.angstrom, - b=2.4677241413662866 * ureg.angstrom, - c=0, - gamma=2.0943951023931957 * ureg.rad + a=2.470 * ureg.angstrom, + b=2.470 * ureg.angstrom, + gamma=2.0943951023931957 * ureg.rad, + pbc=[True, True, False] ) convsystem_C = create_system(label_C_conv, structural_type_C_conv, elements_C_conv, formula_hill_C_conv, formula_reduced_C_conv, formula_anonymous_C_conv, system_relation, material_id=material_id_C_conv, atoms=atoms_C_conv, cell=cell, symmetry=None, prototype=None) topologies_C.append(convsystem_C) return topologies_C -def single_2D_BN_layer_atoms() -> Atoms: +def boron_nitride() -> Atoms: '''Boron nitride system''' symbols_BN = ['B', 'N'] positions_BN = [ @@ -1519,11 +1493,10 @@ def single_2D_BN_layer_atoms() -> Atoms: BN_4 = ase.build.stack(BN_2, BN_2, axis=1) BN_8 = ase.build.stack(BN_4, BN_4, axis=0) BN_16 = ase.build.stack(BN_8, BN_8, axis=1) - BN_16.rattle(stdev=0.001, seed=None, rng=None) return BN_16 -def single_2D_BN_layer_topology() -> List[ResultSystem]: +def boron_nitride_topology() -> List[ResultSystem]: '''Boron nitride topology''' label = 'subsystem' structural_type = '2D' @@ -1563,18 +1536,17 @@ def single_2D_BN_layer_topology() -> List[ResultSystem]: species_N.chemical_symbols = ["N"] species_N.concentration = [1.0] cell = Cell( - a=2.510266994011973 * ureg.angstrom, - b=2.510266994011973 * ureg.angstrom, - c=0, - gamma=2.0943951023931957 * ureg.rad + a=2.513 * ureg.angstrom, + b=2.513 * ureg.angstrom, + gamma=2.0943951023931957 * ureg.rad, + pbc=[True, True, False] ) convsystem_BN = create_system(label, structural_type, elements, formula_hill, formula_reduced, formula_anonymous, system_relation, material_id=material_id, atoms=atoms, cell=cell, symmetry=None, prototype=None) topologies_BN.append(convsystem_BN) return topologies_BN -def single_2D_MoS2_layer_atoms() -> Atoms: - '''MoS2 system''' +def mos2() -> Atoms: symbols_MoS2 = ['Mo', 'S', 'S'] positions_MoS2 = [ [0.0, 0.0, 9.063556323175761], @@ -1595,12 +1567,10 @@ def single_2D_MoS2_layer_atoms() -> Atoms: MoS2_2D = ase.build.surface(system_MoS2, (1, 1, 0), layers=4, vacuum=None, periodic=True) stacked_2D_MoS2 = ase.build.stack(MoS2_2D, MoS2_2D, axis=2, distance=2.5) stacked_2D_MoS2_2 = ase.build.stack(stacked_2D_MoS2, stacked_2D_MoS2, axis=2) - stacked_2D_MoS2_2.rattle(stdev=0.001, seed=None, rng=None) return stacked_2D_MoS2_2 -def single_2D_MoS2_layer_topology() -> List[ResultSystem]: - '''MoS2 topology''' +def mos2_topology() -> List[ResultSystem]: label = 'subsystem' structural_type = '2D' elements = ['Mo', 'S'] @@ -1643,25 +1613,30 @@ def single_2D_MoS2_layer_topology() -> List[ResultSystem]: species_S.chemical_symbols = ["S"] species_S.concentration = [1.0] cell = Cell( - a=3.253646631826119 * ureg.angstrom, - b=3.253646631826119 * ureg.angstrom, - c=0, - gamma=2.0943951023931957 * ureg.rad + a=3.22 * ureg.angstrom, + b=3.22 * ureg.angstrom, + gamma=2.0943951023931957 * ureg.rad, + pbc=[True, True, False] ) convsystem_MoS2 = create_system(label, structural_type, elements, formula_hill, formula_reduced, formula_anonymous, system_relation, material_id=material_id, atoms=atoms, cell=cell, symmetry=None, prototype=None) topologies_MoS2.append(convsystem_MoS2) return topologies_MoS2 -def stacked_C_BN_2D_layers() -> Atoms: - '''System with stacked graphene and boron nitride''' - C_32 = single_2D_graphene_layer_atoms() - topologies_C = single_2D_graphene_layer_topology() - topologies_C[0].indices = [i for i in range(32, 64)] - BN_16 = single_2D_BN_layer_atoms() - stacked_C_BN = ase.build.stack(BN_16, C_32, axis=2, maxstrain=6.7) - stacked_C_BN = ase.build.surface(stacked_C_BN, (0, 0, 1), layers=1, vacuum=10, periodic=True) - stacked_C_BN.rattle(stdev=0.001, seed=None, rng=None) - topologies_BN = single_2D_BN_layer_topology() - topologies_C_BN = topologies_C + topologies_BN - return stacked_C_BN, topologies_C_BN +def stacked_graphene_boron_nitride_topology() -> Atoms: + topologies_c = graphene_topology() + topologies_bn = boron_nitride_topology() + + # Indices are modified + n_atoms_c = len(topologies_c[0].indices) + n_atoms_bn = len(topologies_bn[0].indices) + topologies_c[0].indices = list(range(0, n_atoms_c)) + topologies_bn[0].indices = list(range(n_atoms_c, n_atoms_c + n_atoms_bn)) + + # Lattice parameters are modified as the cells are strained + topologies_c[1].cell.a = 2.16 * ureg.angstrom + topologies_c[1].cell.b = 2.16 * ureg.angstrom + topologies_bn[1].cell.a = 2.16 * ureg.angstrom + topologies_bn[1].cell.b = 2.16 * ureg.angstrom + + return topologies_c + topologies_bn diff --git a/tests/normalizing/test_material.py b/tests/normalizing/test_material.py index f2325f397c9ad667beedbffc0fec16db5ca64dfb..0b6ee9a4c7d41570b9481a68ffde249542c944db 100644 --- a/tests/normalizing/test_material.py +++ b/tests/normalizing/test_material.py @@ -25,6 +25,7 @@ from matid.symmetry.wyckoffset import WyckoffSet # pylint: disable=import-error from nomad.units import ureg from nomad import atomutils from nomad.utils import hash +from nomad.normalizing.common import ase_atoms_from_nomad_atoms, ase_atoms_from_structure from tests.normalizing.conftest import get_template_for_structure @@ -177,25 +178,6 @@ def test_material_2d(two_d): assert material.n_elements == 1 assert material.symmetry is None - # Conventional structure - conv = two_d.results.properties.structures.structure_conventional - assert_structure(conv) - assert conv.n_sites == 2 - assert conv.species_at_sites == ['C', 'C'] - assert np.array_equal(conv.dimension_types, [1, 1, 0]) - assert conv.lattice_parameters.a.to(ureg.angstrom).magnitude == pytest.approx(2.461, abs=1e-3) - assert conv.lattice_parameters.b.to(ureg.angstrom).magnitude == pytest.approx(2.461, abs=1e-3) - assert conv.lattice_parameters.c.to(ureg.angstrom).magnitude == 0 - assert conv.lattice_parameters.alpha is None - assert conv.lattice_parameters.beta is None - assert conv.lattice_parameters.gamma.magnitude == pytest.approx(120 / 180 * np.pi) - - # Original structure - assert_structure(two_d.results.properties.structures.structure_original) - - # Primitive structure - assert_structure(two_d.results.properties.structures.structure_primitive) - def test_material_surface(surface): material = surface.results.material @@ -207,9 +189,6 @@ def test_material_surface(surface): assert material.material_name is None assert material.symmetry is None - properties = surface.results.properties - assert_structure(properties.structures.structure_original) - def test_material_bulk(bulk): # Material @@ -540,25 +519,24 @@ two_d_swap_expected = Atoms( @pytest.mark.parametrize( 'atoms, expected', [ - # 1D with cell boundary in the middle of the structure - (one_d_split, one_d_split_expected), - # 2D with cell boundary in the middle of the structure - (two_d_split, two_d_split_expected), - # 2D with cell where the nonperiodic axis is not last by default in the - # conventional cell. - (two_d_swap, two_d_swap_expected), + pytest.param(one_d_split, one_d_split_expected, id='1D with cell boundary in the middle of the structure'), + pytest.param(two_d_split, two_d_split_expected, id='2D with cell boundary in the middle of the structure'), + pytest.param(two_d_swap, two_d_swap_expected, id='2D cell where the nonperiodic axis is not last by default in the conventional cell.'), ] ) def test_conventional_structure(atoms, expected): '''Tests that the conventional structure has the correct form. ''' entry = get_template_for_structure(atoms) - structure_conventional = entry.results.properties.structures.structure_conventional - pos = structure_conventional.cartesian_site_positions.to(ureg.angstrom).magnitude - cell = structure_conventional.lattice_vectors.to(ureg.angstrom).magnitude - pbc = np.array(structure_conventional.dimension_types, dtype=bool) - - assert np.array_equal(pbc, expected.get_pbc()) - assert np.allclose(pos, expected.get_positions()) - assert np.array_equal(structure_conventional.species_at_sites, expected.get_chemical_symbols()) - assert np.allclose(cell, expected.get_cell()) + topology = entry.results.material.topology + if topology: + for top in topology: + if top.label == 'conventional cell': + conv = ase_atoms_from_nomad_atoms(top.atoms) + else: + conv = ase_atoms_from_structure(entry.results.properties.structures.structure_conventional) + + assert np.array_equal(conv.get_pbc(), expected.get_pbc()) + assert np.allclose(conv.get_positions(), expected.get_positions()) + assert np.array_equal(conv.get_chemical_symbols(), expected.get_chemical_symbols()) + assert np.allclose(conv.get_cell(), expected.get_cell()) diff --git a/tests/normalizing/test_topology.py b/tests/normalizing/test_topology.py index 44c9ade53d7cabccd983f3d21fbfa6f4c95a5caf..b939b21c28104f32ec0136ea4b82220bcddd466b 100644 --- a/tests/normalizing/test_topology.py +++ b/tests/normalizing/test_topology.py @@ -20,23 +20,27 @@ import numpy as np from collections import defaultdict import pytest -from tests.normalizing.conftest import ( +from nomad.units import ureg +from tests.normalizing.conftest import ( # pylint: disable=unused-import get_template_for_structure, get_template_topology, - single_Cu_surface_atoms, - single_Cu_surface_topology, - single_Cr_surface_atoms, - single_Cr_surface_topology, - stacked_Cu_Ni_surface, - single_2D_graphene_layer_atoms, - single_2D_graphene_layer_topology, - single_2D_BN_layer_atoms, - single_2D_BN_layer_topology, - single_2D_MoS2_layer_atoms, - single_2D_MoS2_layer_topology, - stacked_C_BN_2D_layers + conv_bcc, + conv_fcc, + rattle, + stack, + surf, + single_cu_surface_topology, + single_cr_surface_topology, + stacked_cu_ni_surface_topology, + graphene, + graphene_topology, + boron_nitride, + boron_nitride_topology, + mos2, + mos2_topology, + stacked_graphene_boron_nitride_topology, + projection ) -from tests.normalizing.conftest import projection # pylint: disable=unused-import def assert_topology(topology): @@ -81,7 +85,6 @@ def test_topology_calculation(pbc): # Test the original structure original = topology[0] - assert original.structural_type == 'unavailable' assert original.atoms_ref.positions.shape == (6, 3) assert len(original.atoms_ref.labels) == 6 assert original.atoms_ref.lattice_vectors.shape == (3, 3) @@ -165,26 +168,61 @@ def test_no_topology(fixture, request): @pytest.mark.skip @pytest.mark.parametrize('surface, ref_topologies', [ - pytest.param(single_Cu_surface_atoms()[0], single_Cu_surface_topology(), + pytest.param(rattle(surf(conv_fcc('Cu'), [1, 0, 0])), single_cu_surface_topology(), id='single surface Cu FCC 100'), - pytest.param(single_Cu_surface_atoms()[1], single_Cu_surface_topology(), + pytest.param(rattle(surf(conv_fcc('Cu'), [1, 1, 0])), single_cu_surface_topology(), id='single surface Cu FCC 110'), - pytest.param(single_Cr_surface_atoms()[0], single_Cr_surface_topology(), + pytest.param(rattle(surf(conv_bcc('Cr'), [1, 0, 0])), single_cr_surface_topology(), id='single surface Cr BCC 100'), - pytest.param(single_Cr_surface_atoms()[1], single_Cr_surface_topology(), + pytest.param(rattle(surf(conv_bcc('Cr'), [1, 1, 0])), single_cr_surface_topology(), id='single surface Cr BCC 110'), - pytest.param(stacked_Cu_Ni_surface()[0], stacked_Cu_Ni_surface()[1], - id='stacked surfaces of Cu and Ni'), - pytest.param(single_2D_graphene_layer_atoms(), single_2D_graphene_layer_topology(), + pytest.param(rattle( + stack( + surf(conv_fcc('Cu'), [1, 0, 0], vacuum=0), + surf(conv_fcc('Ni'), [1, 0, 0], vacuum=0) + )), + stacked_cu_ni_surface_topology(), + id='stacked surfaces of Cu and Ni'), + pytest.param(rattle(graphene()), graphene_topology(), id='single 2D layer of graphene'), - pytest.param(single_2D_BN_layer_atoms(), single_2D_BN_layer_topology(), + pytest.param(rattle(boron_nitride()), boron_nitride_topology(), id='single 2D layer of BN'), - pytest.param(single_2D_MoS2_layer_atoms(), single_2D_MoS2_layer_topology(), + pytest.param(rattle(mos2()), mos2_topology(), id='single 2D layer of MoS2'), - pytest.param(stacked_C_BN_2D_layers()[0], stacked_C_BN_2D_layers()[1], - id='stacked layer of BN and C') + pytest.param(rattle( + stack( + graphene(), + boron_nitride() + )), + stacked_graphene_boron_nitride_topology(), + id='stacked layer of BN and C' + ) ]) -def test_surface_2D_topology(surface, ref_topologies): +def test_2D_topology(surface, ref_topologies): + + def compare_section(real, ref): + """Used to compare two metainfo sections.""" + for name in ref.m_def.all_quantities.keys(): + compare_quantity(name, real, ref) + + def compare_quantity(name, value, ref): + """Used to compare two metainfo quantities.""" + real_value = getattr(value, name) + ref_value = getattr(ref, name) + if ref_value is None: + assert real_value is None + elif isinstance(ref_value, ureg.Quantity): + assert real_value.magnitude == pytest.approx(ref_value.magnitude, rel=0.01, abs=0) + elif isinstance(ref_value, (np.ndarray, list)): + real_array = np.array(real_value) + ref_array = np.array(ref_value) + if ref_array.dtype == bool: + assert np.array_equal(real_array, ref_array) + else: + assert np.isclose(real_array, ref_array, rtol=0.01, atol=0) + else: + raise ValueError("Could not compare values.") + entry_archive = get_template_for_structure(surface) topology = entry_archive.results.material.topology assert_topology(topology) @@ -209,12 +247,7 @@ def test_surface_2D_topology(surface, ref_topologies): if subsystem_topology.label == 'conventional cell': # Cell - ref_cell = ref_topologies[ref_index].cell - cell = subsystem_topology.cell - if ref_structural_type == '2D': - assert np.allclose(list(cell.values())[:4], list(ref_cell.values()), rtol=1e-05, atol=1e-9) - else: - assert np.allclose(list(cell.values())[:6], list(ref_cell.values()), rtol=1e-05, atol=1e-9) + compare_section(subsystem_topology.cell, ref_topologies[ref_index].cell) # Symmetry if ref_topologies[ref_index].symmetry: @@ -264,7 +297,7 @@ def test_surface_2D_topology(surface, ref_topologies): assert ref_atoms_property == atoms_property elif subsystem_topology.label == 'subsystem': - # Indices: passes if the index overlapp is large enough + # Indices: passes if the index overlap is large enough ref_indices = ref_topologies[ref_index].indices indices = subsystem_topology.indices[0] indices_overlap = set(ref_indices).intersection(set(indices)) @@ -290,7 +323,6 @@ def test_topology_projection(projection): assert len(topology) == 2 assert topology[0].label == 'original' assert topology[1].label == system.atoms_group[-1].label - assert topology[0].structural_type == 'bulk' assert topology[1].structural_type == 'group' assert len(topology[0].child_systems) == 1 assert topology[0].child_systems[-1] == topology[1].system_id