diff --git a/CHANGELOG.md b/CHANGELOG.md index 8c29e6a..e0ddd54 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,2 +1,2 @@ ### 1.0.0 -First Ambuild public release under GPL v3.0 +First AMBUILD public release under GPL v3.0 diff --git a/README.md b/README.md index 74eb7ac..18da778 100644 --- a/README.md +++ b/README.md @@ -2,11 +2,7 @@ Ambuild is a python program for creating polymeric molecular structures. Please feel free to follow Ambuild on [twitter](https://twitter.com/Ambuild2). -The code is developed by [Abbie Trewin's](https://twitter.com/AbbieTrewin) group at the [University of Lancaster](https://www.lancaster.ac.uk/sci-tech/about-us/people/abbie-trewin). - -Detailed instructions on how to run Ambuild can be found on our [wiki](https://github.com/linucks/ambuild/wiki/01_Introduction) - -### Visit our YouTube page to watch our [installation video](https://www.youtube.com/watch?v=cHBXM7fytNM), which talks you through how to install Ambuild on a Ubuntu or Debian machine. Alternatively, read on! +The code is developed by [Abbie Trewin's](https://twitter.com/AbbieTrewin) group at the [University of Lancaster](https://www.lancaster.ac.uk/sci-tech/about-us/people/abbie-trewin). Ambuild has previously been tested on Ubuntu/Debian machines, but should work for other Linux environments too. The instructions given below relate to installation on a Ubuntu/Debian Linux environment. We explain here how to install each component required to run Ambuild onto a new machine. Please feel free to skip any steps describing how to install any components which you have already installed. @@ -25,7 +21,7 @@ With numpy installed Ambuild can be used to create molecular structures, but can #### 2. Install Docker -These instructions are derived from the [Docker Website](https://docs.docker.com/engine/install/ubuntu/). If you already have Docker installed, please skip ahead to [section 4](#4-install-nvidia-drivers). +Instructions from: https://docs.docker.com/engine/install/ubuntu/ 1. Firstly, remove any old versions with the command: ``` @@ -101,7 +97,7 @@ For more examples and ideas, visit: https://docs.docker.com/get-started/ ``` #### 4. Install NVIDIA Drivers -In order for applications to take advantage of GPU acceleration, you will need to install the NVIDIA GPU drivers for your card - the drivers are the piece of software that allow different programmes to communicate with the GPU card. There are instructions for how to do this on the [NVIDIA website](https://docs.nvidia.com/cuda/cuda-installation-guide-linux/index.html) +In order for applications within the Docker container to take advantage of GPU acceleration, you will need to install the NVIDIA GPU drivers for your card - the drivers are the piece of software that allow different programmes to communicate with the GPU card. There are instructions for how to do this on the [NVIDIA website](https://docs.nvidia.com/cuda/cuda-installation-guide-linux/index.html) On Ubuntu, the easiest way to do this seems to be with the command: ``` @@ -109,29 +105,31 @@ sudo ubuntu-drivers autoinstall ``` #### 5. Install NVIDIA Docker Runtime -This is the layer of software that enables the Docker container to talk to the GPU card via the NVIDIA Drivers. - -These instructions are derived from: https://github.com/NVIDIA/nvidia-docker +Instructions from: https://github.com/NVIDIA/nvidia-docker 1. Run the following commands to install the nvidia-container-toolkit: -``` + ``` distribution=$(. /etc/os-release;echo $ID$VERSION_ID) -``` -``` + ``` + ``` curl -s -L https://nvidia.github.io/nvidia-docker/gpgkey | sudo apt-key add - -``` -``` + ``` + ``` curl -s -L https://nvidia.github.io/nvidia-docker/$distribution/nvidia-docker.list | sudo tee /etc/apt/sources.list.d/nvidia-docker.list -``` -``` + ``` + ``` sudo apt-get update && sudo apt-get install -y nvidia-container-toolkit -``` -``` + ``` + ``` sudo systemctl restart docker -``` + ``` + +2. Test nvidia-smi with the latest official CUDA image -2. Install the nvidia-container-runtime + ```docker run --gpus all nvidia/cuda:10.0-base nvidia-smi``` + +3. Install the nvidia-container-runtime ``` sudo apt install nvidia-container-runtime ``` @@ -139,7 +137,7 @@ There is currently a bug with the nvidia docker container runtime as detailed he To work around the bug, carry out the following additional step: -3. Create a file called /etc/docker/daemon.json with the following content by typing the following (cut and paste the following command from 'sudo' to the second 'EOF' into the terminal): +4. Create a file called /etc/docker/daemon.json with the following content by typing the following (cut and paste the following command from 'sudo' to the second 'EOF' into the terminal): ``` sudo tee -a /etc/docker/daemon.json << EOF @@ -153,49 +151,11 @@ sudo tee -a /etc/docker/daemon.json << EOF } EOF ``` -4. Restart docker: +5. Restart docker: ``` sudo systemctl restart docker ``` -5. Test nvidia-smi with the latest official CUDA image - - ```docker run --gpus all nvidia/cuda:10.0-base nvidia-smi``` - -This should generate output similar to the following: -``` -Tue Jul 21 08:54:21 2020 -+-----------------------------------------------------------------------------+ -| NVIDIA-SMI 435.21 Driver Version: 435.21 CUDA Version: 10.1 | -|-------------------------------+----------------------+----------------------+ -| GPU Name Persistence-M| Bus-Id Disp.A | Volatile Uncorr. ECC | -| Fan Temp Perf Pwr:Usage/Cap| Memory-Usage | GPU-Util Compute M. | -|===============================+======================+======================| -| 0 Quadro M4000 Off | 00000000:01:00.0 On | N/A | -| 46% 34C P8 10W / 120W | 43MiB / 8123MiB | 0% Prohibited | -+-------------------------------+----------------------+----------------------+ -| 1 Tesla K20c Off | 00000000:02:00.0 Off | 0 | -| 30% 30C P8 16W / 225W | 0MiB / 4743MiB | 0% Default | -+-------------------------------+----------------------+----------------------+ -``` - -If you do not see an output like the one given above, then there may be an issue with the NVIDIA installation. This can often be caused by conflicts with an existing NVIDIA/CUDA installation. This can often be resolved by uninstalling all previous NVIDIA/CUDA installations and just installing the NVIDA drivers. This can be done with the following commands: - -``` -sudo apt-get --purge remove "*cublas*" "cuda" "nsight" -``` -``` -sudo apt-get --purge remove $(dpkg-query -l | grep nvidia | awk '{print $2}') -``` -``` -sudo apt-get autoremove -``` -``` -sudo ubuntu-drivers autoinstall -``` - -If you still do not see an output as in the example given with step 5, we suggest contacting your local Linux specialist. - #### 6. Disable secondary video GPU card If you have more than one GPU card (e.g. you have a card specifically for running jobs), then you may need to disable your video GPU card for running jobs so that any GPU jobs are placed on the specialised card rather than using the video card. This will not disable the video card for viewing your screen - it will just prevent it being used to run computational simulation jobs. @@ -214,34 +174,41 @@ sudo nvidia-smi -c 2 -i GPU-4030396e-e7b4-aa4d-e035-22758536dba5 ``` E.g. in the example output above, the UUID string will be: ```GPU-66dc2593-494d-4b44-4574-2b92976db56b```, making the command in step 2 read: ```sudo nvidia-smi -c 2 -i GPU-66dc2593-494d-4b44-4574-2b92976db56b``` -#### 7. Get Ambuild - -A tar.gz file of the latest release of Ambuild can be downloaded from the [releases page](https://github.com/linucks/ambuild/releases). +If you do not see an output like the one given as an example in step 1, we advise you to firstly try running the three commands immediately below. If you still do not see an output as in the example given with step 1, we suggest contacting your local Linux specialist. -You can also download the ambuild-1.0.0.tar.gz file directly with the following command: ``` -curl -OL https://github.com/linucks/ambuild/archive/1.0.0.tar.gz +sudo apt-get --purge remove "*cublas*" "cuda" "nsight" +``` +``` +sudo apt-get --purge remove "*nvidia*" ``` - -Once downloaded, the file can be unpacked with the command: ``` -tar -xzf 1.0.0.tar.gz +sudo apt-get autoremove +``` +``` +sudo ubuntu-drivers autoinstall ``` -This will create a directory called ```ambuild-1.0.0``` containing the Ambuild source code. +#### 7. Get Ambuild -Now you are ready to run Ambuild! Please see our [wiki](https://github.com/linucks/ambuild/wiki/08_Running_Ambuild) for how to get started. +1. Firstly, install git using: +``` +sudo apt-get install git +``` + +2. Checkout ambuild from the GitHub repository. We recommend that you cd into /opt, and then run this command. +``` +git clone https://github.com/linucks/ambuild.git +``` ## Installation of optional dependencies ### Poreblazer To install [poreblazer](https://github.com/richardjgowers/poreblazer) for use by Ambuild, the following steps are required. 1. Checkout or download poreblazer from GitHub: -```git clone https://github.com/SarkisovGroup/PoreBlazer``` +```git clone https://github.com/richardjgowers/poreblazer.git``` 2. Install the [gfortran](https://gcc.gnu.org/wiki/GFortran) compiler. On Ubuntu/Debian, this should just be a case of running: ```sudo apt-get install gfortran``` -3. Compile the poreblazer executable. This is done in the ```src``` directory of the poreblazer directory, so cd into this directory and open the ```Makefile```. Then, edit ```FORTRAN_COMPILER``` line so it reads: ```FORTRAN_COMPILER= gfortran``` - -4. Run the command: ```make``` This should create the ```poreblazer.exe``` executable in this directory. +3. Compile the poreblazer executable. This is done in the ```src``` directory of the poreblazer directory, so cd into this directory and then run the command: ```make``` This should create the ```poreblazer.exe``` executable in this directory. diff --git a/ambuild/ab_analyse.py b/ambuild/ab_analyse.py index 01bc3a4..3a8f439 100644 --- a/ambuild/ab_analyse.py +++ b/ambuild/ab_analyse.py @@ -48,7 +48,7 @@ def __init__(self, cell, logfile="ambuild.csv"): def start(self): """Called whenever we start a step""" - assert self._stepTime is None + assert self._stepTime == None assert self.last self.step += 1 self._stepTime = time.time() diff --git a/ambuild/ab_block.py b/ambuild/ab_block.py index 35a65c8..5a0b77c 100644 --- a/ambuild/ab_block.py +++ b/ambuild/ab_block.py @@ -28,7 +28,6 @@ from ambuild.ab_fragment import Fragment from ambuild import xyz_core from ambuild import xyz_util -from ambuild import ab_fragment logger = logging.getLogger(__name__) @@ -43,36 +42,17 @@ class Block(object): * bonds just point to two endGroups that are bonded to each other """ - def __init__( - self, - filePath=None, - fragmentType=None, - initFragment=None, - solvent=None, - markBonded=None, - catalyst=None, - ): + def __init__(self, filePath=None, fragmentType=None, initFragment=None): """ Constructor """ # Need to change so cannot create block withough fragmentType - self.fragments = [] if filePath: - if initFragment: - raise RuntimeError("Cannot include initFragment with filePath") - if not os.path.isfile(filePath): - raise RuntimeError(f"Missing fragment file: {filePath}") - if not fragmentType: - raise RuntimeError( - f"Missing fragmentType definition for file: {filePath}" - ) - initFragment = ab_fragment.fragmentFactory( - fragmentType, - filePath, - solvent=solvent, - markBonded=markBonded, - catalyst=catalyst, + assert os.path.isfile(filePath) and fragmentType, "Missing file: {}".format( + filePath ) + initFragment = Fragment(filePath, fragmentType) + self.fragments = [] if initFragment: self.fragments.append(initFragment) self._blockBonds = [] # List of bond objects between blocks @@ -85,7 +65,6 @@ def __init__( self._bodies = ( [] ) # List of which body in the block each atom belongs too - frags can contain multiple bodies - self._fragmentTypeDict = {} # The list of atoms that are endGroups and their corresponding angleAtoms self._freeEndGroups = {} self._numFreeEndGroups = 0 @@ -286,7 +265,8 @@ def blockBonds(self): ] def bondBlock(self, bond): - """Add newBlock to this one""" + """ Add newBlock to this one + """ assert bond.endGroup1.block() == self # Tried optimising this by passing in the bond to update and only updating those fragments/ # endGroups that had changed but it rapidly got rather complicated so we keep to a simple @@ -366,13 +346,10 @@ def copy(self): def dataByFragment(self, fragmentType): """Return the data for a specific fragmentType within the block""" - if fragmentType not in self.fragmentTypes(): - logger.debug( - f"Cannot find fragmentType '{fragmentType}' in block: {self.fragmentTypes()}" - ) - return None + coords = [] symbols = [] + atomCount = 0 fbondRen = {} for i in range(len(self._dataMap)): @@ -381,15 +358,14 @@ def dataByFragment(self, fragmentType): coords.append(self.coord(i)) symbols.append(self.symbol(i)) atomCount += 1 - if atomCount > 0: - bonds = [ - (fbondRen[b1], fbondRen[b2]) - for ftype, (b1, b2) in self._bondsByFragmentType - if ftype == fragmentType - ] - return coords, symbols, bonds - else: - return None + + bonds = [ + (fbondRen[b1], fbondRen[b2]) + for ftype, (b1, b2) in self._bondsByFragmentType + if ftype == fragmentType + ] + + return coords, symbols, bonds def deleteBond(self, bond, root=None): """root is an optional fragment which we want to stay in this block""" @@ -584,7 +560,8 @@ def dihedrals(self, atom1Idx, atom2Idx, bondOnly=False): return dindices def flip(self, fvector): - """Rotate perpendicular to fvector so we facing the opposite way along the fvector""" + """Rotate perpendicular to fvector so we facing the opposite way along the fvector + """ # Find vector perpendicular to the bond axis # Dot product needs to be 0 @@ -631,7 +608,7 @@ def fragmentTypeDict(self): return self._fragmentTypeDict def fragmentTypes(self): - return list(self._fragmentTypeDict.keys()) + return self._fragmentTypeDict.keys() def freeEndGroups(self, endGroupTypes=None, fragment=None): """Return the list of free endGroups. @@ -667,7 +644,8 @@ def freeEndGroupTypes(self): return self._endGroupType2EndGroups.keys() def isEndGroup(self, idxAtom): - """Return True if this atom is a free endGroup""" + """Return True if this atom is a free endGroup + """ # No need to do conversion as atomEndGroups is external interface if self.atomEndGroups(idxAtom): return True @@ -683,33 +661,31 @@ def maxAtomRadius(self): assert self._maxAtomRadius > 0 return self._maxAtomRadius - def newBondPosition(self, staticEndGroup, growEndGroup): - """Return the position where a bond to an atom of type growEndGroup - would be placed if bonding to the staticEndGroup + def newBondPosition(self, endGroup, symbol): + """Return the position where a bond to an atom of type 'symbol' + would be placed if bonding to the target endgroup + I'm sure this algorithm is clunky in the extreme... """ + targetEndGroup = self.coord(endGroup.blockEndGroupIdx) + targetSymbol = self.symbol(endGroup.blockEndGroupIdx) + targetCapAtom = self.coord(endGroup.blockCapIdx) + # Get the bond length between these two atoms - growBlock = growEndGroup.block() - growAtomType = growBlock.type(growEndGroup.blockEndGroupIdx) - staticAtomType = self.type(staticEndGroup.blockEndGroupIdx) - bondLength = xyz_util.bondLength(staticAtomType, growAtomType) - if bondLength < 0: - growAtomSymbol = growBlock.symbol(growEndGroup.blockEndGroupIdx) - staticSymbol = self.symbol(staticEndGroup.blockEndGroupIdx) - bondLength = xyz_util.bondLength(staticSymbol, growAtomSymbol) - - # Find unit vector pointing from targetAngleAtom to staticEndGroup - staticCoord = self.coord(staticEndGroup.blockEndGroupIdx) - staticCapAtomCoord = self.coord(staticEndGroup.blockCapIdx) - # vector from staticEndGroup to staticCapAtom - # v1 = staticEndGroup - staticCapAtom - v1 = staticCapAtomCoord - staticCoord + bondLength = xyz_util.bondLength(targetSymbol, symbol) + + # Find unit vector pointing from targetAngleAtom to targetEndGroup + + # vector from targetEndgroup to targetCapAtom + # v1 = targetEndGroup - targetCapAtom + v1 = targetCapAtom - targetEndGroup # Now get unit vector uv = v1 / np.linalg.norm(v1) # Multiply unit vector by bond length to get the component to add on - newPosition = staticCoord + (uv * bondLength) + newPosition = targetEndGroup + (uv * bondLength) + return newPosition def numFreeEndGroups(self): @@ -731,10 +707,15 @@ def positionGrowBlock(self, endGroup, growEndGroup, dihedral=None): endGroupAtom = self.coord(endGroup.blockEndGroupIdx) capAtom = self.coord(endGroup.blockCapIdx) refVector = endGroupAtom - capAtom + growBlock = growEndGroup.block() # get the coord where the next block should bond - bondPos = self.newBondPosition(endGroup, growEndGroup) + # symbol of endGroup tells us the sort of bond we are making which determines + # the bond length + symbol = growBlock.symbol(growEndGroup.blockEndGroupIdx) + bondPos = self.newBondPosition(endGroup, symbol) + # print "got bondPos for {0}: {1}".format( symbol, bondPos ) # Align along the staticBlock bond growBlock.alignAtoms( @@ -776,8 +757,8 @@ def positionGrowBlock(self, endGroup, growEndGroup, dihedral=None): def randomRotate(self, origin=[0, 0, 0], atOrigin=False): """Randomly rotate a block. - Args: - atOrigin -- flag to indicate if the block is already positioned at the origin + Args: + atOrigin -- flag to indicate if the block is already positioned at the origin """ if not atOrigin: position = self.centroid() @@ -824,7 +805,7 @@ def solvent(self): def selectEndGroup(self, endGroupTypes=None, random=True): """Return a random free endGroup in the block""" - if endGroupTypes is None: + if endGroupTypes == None: if random: # We pick a random endGroup endGroup = _random.choice(self.freeEndGroups()) @@ -863,7 +844,7 @@ def selectEndGroup(self, endGroupTypes=None, random=True): return endGroup def translate(self, tvector): - """translate the molecule by the given vector""" + """ translate the molecule by the given vector""" # CHANGE SO WE CHECK IF IS A NUMPY ARRAY if isinstance(tvector, list): tvector = np.array(tvector) @@ -881,7 +862,8 @@ def translateCentroid(self, position): return def _update(self): - """Set the list of _endGroups & update data for new block""" + """Set the list of _endGroups & update data for new block + """ # Now build up the dataMap listing where each fragment starts in the block and linking the # overall block atom index to the fragment and fragment atom index self._dataMap = [] diff --git a/ambuild/ab_body.py b/ambuild/ab_body.py index b955e44..2c59030 100644 --- a/ambuild/ab_body.py +++ b/ambuild/ab_body.py @@ -56,8 +56,6 @@ def coordsRelativeToCom(self): @property def diameters(self): - #return self.radii * 2.0 - # We currently return a diameter of 0.1 as we don't want this influencing the simulation return [xyz_util.DUMMY_DIAMETER] * self.natoms @property @@ -88,10 +86,6 @@ def masses(self): def principalMoments(self): return xyz_core.principalMoments(self.coords, self.masses) - @property - def radii(self): - return np.compress(self.indexes, self.fragment._radii, axis=0) - @property def rigidConfigStr(self): """return the type of this body based on the endGroup configuration""" @@ -100,12 +94,10 @@ def rigidConfigStr(self): @property def static(self): if hasattr(self.fragment, "static") and self.fragment.static: - return True - return False - - @property - def staticAtoms(self): - return [self.static] * self.natoms + v = True + else: + v = False + return [v] * self.natoms @property def symbols(self): diff --git a/ambuild/ab_cell.py b/ambuild/ab_cell.py index 52048f1..e871d73 100644 --- a/ambuild/ab_cell.py +++ b/ambuild/ab_cell.py @@ -53,7 +53,7 @@ def __init__( """Construct an empty cell: Args: - boxDim - a list with three numbers specifying the size of the cell A,B and C dimensions (Angstroms) + boxDim - a list with three numbers specifying the size of the cell A,B and C dimensions (angstroms) Dimensions are from 0 - A, B or C filePath - a (.car) file with the cell dimensions and the coordinates of molecules that will be kept static throught the simulation. @@ -62,41 +62,40 @@ def __init__( bondMargin - two atoms are considered close enough to bond if they are within the bond length defined for the two atoms +/- the bondMargin bondAngleMargin - the tolerance (in degrees) from the ideal of 180 that defines an acceptable bond - paramsDir - path to the directory holding the forcefield parameter csv files (default ../params) debugLog - True/False - specifies if a log will be created - not recommended as it generates lots of data and slows the program. - + paramsDir - path to the directory holding the forcefield parameter csv files (default ../params) """ - self.origin = np.array( - [0, 0, 0], dtype=np.float64 - ) # For time being origin always 0,0,0 - self.dim = None + # For time being origin always 0,0,0 + self.origin = np.array([0, 0, 0], dtype=np.float64) + self.dim = None # The cell dimensions self.pbc = [True, True, True] self.walls = [False, False, False] self.wallAtomType = None self.wallRadius = ( None # (determind from the covalent radius of the wallAtomType) ) - # additional distance to add on to the characteristic bond length when checking whether 2 atoms are close enough to bond + # additional distance to add on to the characteristic bond length + # when checking whether 2 atoms are close enough to bond self.bondMargin = bondMargin - # additional distance to add on to the atom covalent radii when checking if two atoms are close enough to clash + # additional distance to add on to the atom covalent radii when checking if two atoms + # are close enough to clash self.atomMargin = atomMargin # convert bondAngle and bondMargin to angstroms + # Could check values are in degrees and not radians? self.bondAngleMargin = math.radians(bondAngleMargin) self.targetDensity = 10 self.targetEndGroups = 100 # number of free endgroups left - self.box1 = ( - {} - ) # Dict mapping box key to a list of tuples defining the atoms in that box - self.box3 = {} # Dict mapping key to list of boxes surrounding the keyed box + # Dict mapping box key to a list of tuples defining the atoms in that box + self.box1 = {} + # Dict mapping key to list of boxes surrounding the keyed box + self.box3 = {} + # max atom radius - used to calculate box size self.boxSize = None - self.maxAtomRadius = -1 # max atom radius - used to calculate box size + self.maxAtomRadius = -1 self.rCut = 5.0 - self.numBoxes = [ - None, - None, - None, - ] # number of boxes in A,B,C axes - used for calculating PBC + # number of boxes in A,B,C axes - used for calculating PBC + self.numBoxes = [None, None, None] self._fragmentLibrary = {} # fragmentType -> parentFragment self._endGroup2LibraryFragment = {} # endGroup type -> parentFragment self.bondTypes = [] @@ -324,7 +323,8 @@ def attachBlock(self, growEndGroup, staticEndGroup, dihedral=None): return False def bondAllowed(self, endGroup1, endGroup2): - """Check if the given bond is permitted from the types of the two fragments""" + """Check if the given bond is permitted from the types of the two fragments + """ # logger.debug( "checking bondAllowed {0} {1}".format( ftype1, ftype2 ) ) if endGroup1.fragmentType() == "cap" or endGroup2.fragmentType() == "cap": assert False, "NEED TO FIX CAPS!" @@ -342,7 +342,8 @@ def bondAllowed(self, endGroup1, endGroup2): return False def bondBlock(self, bond, selfBond=True): - """Bond the second block1 to the first and update the data structures""" + """ Bond the second block1 to the first and update the data structures + """ logger.debug("cell bondBlock: {0}".format(bond)) # logger.debug("before bond: {0} - {1}".format( bond.idxBlock1, bond.block1._bondObjects) ) selfBond = True # HACK @@ -360,7 +361,8 @@ def bondBlock(self, bond, selfBond=True): return True def bondClash(self, bond, clashDist): - """Check if any atoms are clashDist from bond.""" + """Check if any atoms are clashDist from bond. + """ if clashDist > self.boxSize: raise RuntimeError( "clashDist needs to be less than the boxSize: {0}".format(self.boxSize) @@ -451,6 +453,7 @@ def _cat1Paf2(self, bond, fragmentTypes): assert len(endGroups) == 2, "Assumption is CAT only has 2 endGroups!" catEG1, catEG2 = endGroups logger.info("_cat1Paf2 processing bond %s" % bond) + logger.debug("_cat1Paf2 %s %s" % (catEG1, catEG2)) # Now get the two paf endGroups that are bonded to the cat EndGroups cat = catEG1.block() @@ -482,11 +485,7 @@ def _cat1Paf2(self, bond, fragmentTypes): ) ) return - return self._joinPaf( - fragmentTypes, - bond1, - bond2, - ) + return self._joinPaf(fragmentTypes, bond1, bond2) def _cat2Paf2(self, cc_bond, fragmentTypes): """Function to unbond a Ni-catalyst bonded to two PAF groups @@ -585,11 +584,9 @@ def _cat2Paf2(self, cc_bond, fragmentTypes): self.addBlock(cat1) # Run optimisation to move CAT away - logger.info("** _cat2Paf2 Optimisation **") + logger.info("_cat2Paf2 Optimisation") # self.dump() - self.optimiseGeometry( - rigidBody=True, quiet=True, dt=0.000001, optCycles=1000000 - ) + self.optimiseGeometry(rigidBody=True, quiet=True, dt=0.00001, optCycles=1000000) # Now dealing with a CAT bonded to two PAF groups # Need to select the other cat-paf bond @@ -685,10 +682,9 @@ def cat1Paf2(self, fragmentTypes, dt=0.00001, optCycles=1000000): logger.info( "cat1Paf2 undertook %d processes. Now running optimisation", nbonds ) - if optCycles > 0: - self.optimiseGeometry( - rigidBody=True, dt=dt, optCycles=optCycles, dump=False, max_tries=1 - ) + self.optimiseGeometry( + rigidBody=True, dt=dt, optCycles=optCycles, dump=False, max_tries=1 + ) self.clearUnbonded() return True return False @@ -705,13 +701,12 @@ def cat2Paf2(self, fragmentTypes, dt=0.00001, optCycles=1000000): # logger.info("cat2Paf2 got new bonds %s" % [str(b) for b in self.newBonds ]) nbonds = len([self._cat2Paf2(b, fragmentTypes) for b in self.newBonds]) if nbonds > 0: - if optCycles > 0: - logger.info( - "cat2Paf2 undertook %d processes. Now running optimisation", nbonds - ) - self.optimiseGeometry( - rigidBody=True, dt=dt, optCycles=optCycles, dump=False, max_tries=1 - ) + logger.info( + "cat2Paf2 undertook %d processes. Now running optimisation", nbonds + ) + self.optimiseGeometry( + rigidBody=True, dt=dt, optCycles=optCycles, dump=False, max_tries=1 + ) self.clearUnbonded() return True return False @@ -726,24 +721,17 @@ def canBond( bondMargin, bondAngleMargin, ): - # The check should have been made before this is called on whether the two atoms are endGroups + # The check should have been made before this is called on whether the two atoms are endGroup + # Check length bond_length = xyz_util.bondLength( - addBlock.type(idxAddAtom), staticBlock.type(idxStaticAtom) + addBlock.symbol(idxAddAtom), staticBlock.symbol(idxStaticAtom) ) if bond_length < 0: - bond_length = xyz_util.bondLength( - addBlock.symbol(idxAddAtom), staticBlock.symbol(idxStaticAtom) - ) - if bond_length < 0: - raise RuntimeError( - "Missing bond distance for: {0}-{1} or {2}-{3}".format( - addBlock.symbol(idxAddAtom), - staticBlock.symbol(idxStaticAtom), - addBlock.type(idxAddAtom), - staticBlock.type(idxStaticAtom), - ) + raise RuntimeError( + "Missing bond distance for: {0}-{1}".format( + addBlock.symbol(idxAddAtom), staticBlock.symbol(idxStaticAtom) ) - + ) # See if the distance between them is acceptable # print "CHECKING BOND ATOMS ",bond_length,self.distance( addCoord, staticCoord ) if not ( @@ -768,6 +756,10 @@ def canBond( ) ) continue + # print "Possible bond for {0} {1} {2} dist {3}".format( idxAddAtom, + # idxStaticBlock, + # idxStaticAtom, + # self.distance( addCoord, staticCoord ) ) addCoord = addBlock.coord(idxAddAtom) staticCoord = staticBlock.coord(idxStaticAtom) addCapAtom = addBlock.coord(addEndGroup.capIdx()) @@ -775,6 +767,11 @@ def canBond( staticCapAtom = staticBlock.coord(staticEndGroup.capIdx()) angle2 = self.angle(staticCapAtom, staticCoord, addCoord) + # print "CHECKING ANGLE BETWEEN: {0} | {1} | {2}".format( angleAtom, addCoord, staticCoord ) + # print "{} < {} < {}".format( (self.bondAngle-bondAngleMargin) * util.RADIANS2DEGREES, + # angle * util.RADIANS2DEGREES, + # (self.bondAngle+bondAngleMargin) * util.RADIANS2DEGREES ) + # Check if atoms are in line (zero degrees) within margin if not ( ( @@ -862,7 +859,7 @@ def _checkMove(self, idxAddBlock): if wallClashes: logger.debug("_checkMove got clash with wall") return 1 - if len(close) == 0: + if len(close) is 0: logger.debug("_checkMove no close contacts") return 0 addBlock = self.blocks[idxAddBlock] @@ -1144,14 +1141,11 @@ def cellEndGroupPair(self, cellEndGroups=None): def cellData( self, rigidBody=True, + periodic=True, center=False, fragmentType=None, noRigidParticles=False, ): - d = ab_celldata.CellData(self.dim) # Object to hold the cell data - if fragmentType is not None: - return self._cellDataByFragment(d, fragmentType) - RIGIDPARTICLES = ( rigidBody and ab_util.HOOMDVERSION and ab_util.HOOMDVERSION[0] > 1 ) @@ -1159,7 +1153,24 @@ def cellData( RIGIDPARTICLES = False if RIGIDPARTICLES: self.rigidParticleMgr.reset() - + # Object to hold the cell data + d = ab_celldata.CellData() + d.cell = self.dim + if fragmentType is not None: + if RIGIDPARTICLES: + assert False, "Need to update for HOOMD2" + # Only returning data for one type of fragment + assert ( + fragmentType in self.fragmentTypes() + ), "FragmentType {0} not in cell!".format(fragmentType) + atomIdx = 0 + for b in self.blocks.values(): + coords, symbols, bonds = b.dataByFragment(fragmentType) + d.coords += coords + d.symbols += symbols + d.bonds += [(b1 + atomIdx, b2 + atomIdx) for (b1, b2) in bonds] + atomIdx += len(coords) + return d atomIdx = 0 bodyIdx = 0 for block in self.blocks.values(): @@ -1248,12 +1259,11 @@ def cellData( for body in frag.bodies(): d.natoms += body.natoms atomIdx += body.natoms - if RIGIDPARTICLES and not body.static: + if RIGIDPARTICLES: d.rigidParticles.append( self.rigidParticleMgr.createParticle(body) ) else: - # Either not using RIGIDPARTICLES or using RIGIDPARTICLES and this body is static coords = body.coords coords, images = xyz_core.wrapCoord3( coords, dim=self.dim, center=center @@ -1263,10 +1273,10 @@ def cellData( d.atomTypes += body.atomTypes d.bodies += [bodyIdx] * body.natoms d.charges += body.charges - d.diameters += list(body.diameters) + d.diameters += body.diameters d.masked += body.masked d.masses += list(body.masses) - d.static += body.staticAtoms + d.static += body.static d.symbols += body.symbols bodyIdx += 1 if RIGIDPARTICLES: @@ -1274,26 +1284,6 @@ def cellData( d.rigidParticleMgr = self.rigidParticleMgr return d - def _cellDataByFragment(self, d, fragmentType): - """Return data just for the given type of fragment""" - # Only returning data for one type of fragment - assert ( - fragmentType in self.fragmentTypes() - ), "FragmentType {0} not in cell!".format(fragmentType) - atomIdx = 0 - for b in self.blocks.values(): - data = b.dataByFragment(fragmentType) - if data: - coords, symbols, bonds = data - d.coords += coords - d.symbols += symbols - d.bonds += [(b1 + atomIdx, b2 + atomIdx) for (b1, b2) in bonds] - atomIdx += len(coords) - if atomIdx == 0: - return None - else: - return d - def delBlock(self, blockId): """ Remove the block with the given index from the cell @@ -1303,7 +1293,7 @@ def delBlock(self, blockId): keys = [] for iatom, key in enumerate(block.atomCell): # Skip dummy atoms - if key is None: + if key == None: continue keys.append(key) self.box1[key].remove((blockId, iatom)) @@ -1324,7 +1314,7 @@ def deleteBlocks(self, fragmentTypes=None, indices=None, maxFrags=1, numBlocks=0 fragmentTypes: the fragmentType, or list of fragmentTypes of the blocks to remove indices: list of indices of the blocks in the overall list of blocks in the cell to be deleted maxFrags: only blocks containing <= this number of fragments will be removed (default = 1) - numBlocks: optional - number of blocks to remove, otherwise all blocks of specified fragmentTypes are removed + numBlocks: optional - the number of blocks to remove, otherwise all blocks of the specified fragmentTypes will be removed Returns: A list of the blocks that were removed - suitable for re-adding with restoreBlocks @@ -1889,7 +1879,7 @@ def joinBlocks(self, toJoin, cellEndGroups=None, dihedral=None, maxTries=100): moveEndGroup, staticEndGroup = self.cellEndGroupPair( cellEndGroups=cellEndGroups ) - if moveEndGroup is None or staticEndGroup is None: + if moveEndGroup == None or staticEndGroup == None: logger.critical("joinBlocks cannot join any more blocks") return added # Copy the original block so we can replace it if the join fails @@ -1946,9 +1936,9 @@ def libraryAddFragment( ) ) # Create fragment - frag = ab_fragment.fragmentFactory( - fragmentType, + frag = ab_fragment.Fragment( filename, + fragmentType, solvent=solvent, markBonded=markBonded, catalyst=catalyst, @@ -2135,16 +2125,12 @@ def optimiseGeometry( """ logger.info("Running optimisation") if not self.mdEngineCls: - import warnings - - warnings.warn("No mdEngine defined - cannot run MD.") - return - # raise RuntimeError("No mdEngine defined - cannot run MD.") + raise RuntimeError("No mdEngine defined - cannot run MD.") mdEngine = self.mdEngineCls(self.paramsDir) if doDihedral and doImproper: raise RuntimeError("Cannot have impropers and dihedrals at the same time") self.setRcut(rigidBody, mdEngine, kw) - data = self.cellData(center=True, rigidBody=rigidBody) + data = self.cellData(periodic=True, center=True, rigidBody=rigidBody) d = {} # for printing results ok = mdEngine.optimiseGeometry( data, @@ -2276,8 +2262,8 @@ def positionBlock( self, block, margin=None, point=None, radius=None, zone=None, random=True ): """Randomly move the given block - If margin is given, use this as a buffer from the edges of the cell - when selecting the coord + If margin is given, use this as a buffer from the edges of the cell + when selecting the coord """ if point: assert len(point) == 3, "Point needs to be a list of three floats!" @@ -2414,7 +2400,7 @@ def runMD( raise RuntimeError("Cannot have impropers and dihedrals at the same time") self.setRcut(rigidBody, mdEngine, kw) d = {} - data = self.cellData(center=True, rigidBody=rigidBody) + data = self.cellData(periodic=True, center=True, rigidBody=rigidBody) ok = mdEngine.runMD( data, xmlFilename=xmlFilename, @@ -2485,7 +2471,7 @@ def seed( zone=None, random=True, ): - """Seed a cell with nblocks of type fragmentType. + """ Seed a cell with nblocks of type fragmentType. Args: nblocks - the number of blocks to add. @@ -2609,8 +2595,8 @@ def setBoxSize(self, boxDim): def setStaticBlock(self, filePath, replace=False): # Create fragment - fragmentType = os.path.splitext(os.path.basename(filePath))[0] - f = ab_fragment.fragmentFactory(fragmentType, filePath=filePath, static=True) + name = os.path.splitext(os.path.basename(filePath))[0] + f = ab_fragment.Fragment(filePath, fragmentType=name, static=True) p = f.cellParameters() if not p: raise RuntimeError( @@ -2725,7 +2711,8 @@ def setRcut(self, rigidBody, mdEngine, kw): return def setWall(self, XOY=False, XOZ=False, YOZ=False, wallAtomType="c"): - """Create walls along the specified sides.""" + """Create walls along the specified sides. + """ self.walls = [XOY, XOZ, YOZ] self.pbc = [not b for b in self.walls] @@ -2870,7 +2857,8 @@ def writePickle(self, fileStem, compress=True): return fileName def writeCar(self, ofile="ambuild.car", data=None, periodic=True, skipDummy=False): - """Car File""" + """Car File + """ if not data: data = self.cellData() car = "!BIOSYM archive 3\n" @@ -2949,7 +2937,7 @@ def writeXyz(self, ofile, data=None, periodic=False, atomTypes=False): If label is true we write out the atom label and block, otherwise the symbol """ if data is None: - d = self.cellData() + d = self.cellData(periodic=periodic, fragmentType=None) else: d = data if periodic: @@ -3115,7 +3103,16 @@ def zipBlocks( if not len(self._possibleBonds): logger.info("zipBlocks: No bonds remaining after clash checks") return 0 - logger.info("zipBlocks found %d potential bonds", len(self._possibleBonds)) + logger.info("zipBlocks found %d additional bonds", len(self._possibleBonds)) + # for b in self._possibleBonds: + # print "Attempting to bond: {0} {1} {2} -> {3} {4} {5}".format( b.block1.id(), + # b.endGroup1.blockEndGroupIdx, + # b.block1.atomCoord( b.endGroup1.blockEndGroupIdx), + # b.block2.id(), + # b.endGroup2.blockEndGroupIdx, + # b.block2.atomCoord( b.endGroup2.blockEndGroupIdx), + # ) + # todo = len(self._possibleBonds) bondsMade = self.processBonds(selfBond=selfBond) if bondsMade != todo: @@ -3142,7 +3139,7 @@ def __getstate__(self): return d def __setstate__(self, d): - """Called when we are unpickled""" + """Called when we are unpickled """ self.__dict__.update(d) if "logfile" in d: # Hack for older versions with no logfile attribute logfile = ab_util.newFilename(d["logfile"]) diff --git a/ambuild/ab_celldata.py b/ambuild/ab_celldata.py index b6e827f..c10eaef 100644 --- a/ambuild/ab_celldata.py +++ b/ambuild/ab_celldata.py @@ -1,6 +1,6 @@ class CellData(object): - def __init__(self, cell): - self.cell = cell + def __init__(self): + self.cell = [] self.natoms = 0 self.atomTypes = [] self.bodies = [] diff --git a/ambuild/ab_endgroup.py b/ambuild/ab_endgroup.py index 187c8cc..da26711 100644 --- a/ambuild/ab_endgroup.py +++ b/ambuild/ab_endgroup.py @@ -4,12 +4,8 @@ @author: abbietrewin """ import logging -import warnings - import numpy as np -from ambuild import xyz_util - ENDGROUPBONDED = "*" logger = logging.getLogger() @@ -32,13 +28,17 @@ def __init__(self): self.blockDihedralIdx = -1 self.fragmentUwIdx = -1 self.blockUwIdx = -1 - self.triAtoms = None - self.triDistances = None return def block(self): - if self.fragment is None or self.fragment.block is None: - raise RuntimeError("endGroup without fragment or block") + f = True + b = True + if self.fragment.block is None: + b = False + if self.fragment is None: + f = False + if not b and f: + raise RuntimeError("None Block {0} Fragment {1}\n{2}".format(b, f, self)) return self.fragment.block def capIdx(self): @@ -75,15 +75,18 @@ def setBonded(self, bond): return def unBond(self, bondEndGroup): - """See: - https://github.com/akshayb6/trilateration-in-3d/blob/master/trilateration.py - """ - - def unBondSimple(bondEndGroup): - """Fallback code - to be avoided as it will calculate incorrect positions - if the bond vector had moved due to the blocks moving through MD""" - # + self.bonded = False + # HACK WE REMOVE ALL SUFFIXES + for eg in self.fragment.endGroups(): + if eg._endGroupType.endswith(ENDGROUPBONDED): + logger.debug("unBond ENDGROUPBONDED") + eg._endGroupType = eg._endGroupType.rstrip(ENDGROUPBONDED) + self.fragment.delBond(self.type()) + # Unmask cap and set the coordinate to the coordinate of the last block atom + # NOTE - NEED TO SCALE BY CORRECT LENGTH + if hasattr(bondEndGroup, "coord"): # Reposition the cap atom based on the bond vector + # self.fragment._coords[self.fragmentCapIdx] = bondEndGroup.coord() # Get vector from this endGroup to the other endGroup egPos = self.fragment._coords[self.fragmentEndGroupIdx] v1 = bondEndGroup.coord() - egPos @@ -93,36 +96,6 @@ def unBondSimple(bondEndGroup): self.fragment._coords[self.fragmentCapIdx] = egPos + ( uv * self.capBondLength ) - return - - self.bonded = False - # HACK WE REMOVE ALL SUFFIXES - for eg in self.fragment.endGroups(): - if eg._endGroupType.endswith(ENDGROUPBONDED): - eg._endGroupType = eg._endGroupType.rstrip(ENDGROUPBONDED) - self.fragment.delBond(self.type()) - if hasattr(bondEndGroup, "coord"): - # Reposition the cap atom - if self.triAtoms is not None: - # Get positions of triAtoms - triAtoms = [self.fragment._coords[i] for i in self.triAtoms] - # Calculate the new position - try: - self.fragment._coords[self.fragmentCapIdx] = xyz_util.trilaterate3D( - self.triDistances, triAtoms - ) - except Exception as err: - logger.warning( - "Failed to trilaterate position for:\n" - + f"Distances: {self.triDistances}\n" - + f"Positions: {triAtoms}\n" - + f"Error was: {err}" - ) - logger.warning("** Reverting to simple unbond **") - unBondSimple(bondEndGroup) - else: - warnings.warn(f"Cannot properly unbond {bondEndGroup} as no triAtoms") - unBondSimple(bondEndGroup) # Unhide the cap atom self.fragment.masked[self.fragmentCapIdx] = False diff --git a/ambuild/ab_ffield.py b/ambuild/ab_ffield.py index 3a52383..5c08f3a 100755 --- a/ambuild/ab_ffield.py +++ b/ambuild/ab_ffield.py @@ -7,7 +7,6 @@ logger = logging.getLogger(__name__) PARAM_SEP = "," -ANGLE_SEP = "-" # Needs to be defined here to avoid circular import in util.py def read_bond_params(bond_file): @@ -97,9 +96,6 @@ def _processAngles(self, angle_file): k = float(row[1]) t0 = math.radians(float(row[2])) angles[angle] = {"k": k, "t0": t0} - # Angles are symmetrical so we also need to add the opposite notation - angle_reversed = ANGLE_SEP.join(angle.split(ANGLE_SEP)[::-1]) - angles[angle_reversed] = {"k": k, "t0": t0} except Exception as e: logger.critical( "Error reading angle parameters from file: {0}".format( @@ -268,7 +264,7 @@ def hasImproper(self, improper): return improper in self.impropers.keys() def pairParameter(self, p1, p2): - """Return whichever pair is defined""" + """ Return whichever pair is defined """ if (p1, p2) in self.pairs: return self.pairs[(p1, p2)] if (p2, p1) in self.pairs: @@ -278,7 +274,7 @@ def pairParameter(self, p1, p2): return def hasPair(self, p1, p2): - """dummy atoms treated specially""" + """ dummy atoms treated specially """ if p1.lower() == "x" or p2.lower() == "x": return True if (p1, p2) in self.pairs or (p2, p1) in self.pairs: diff --git a/ambuild/ab_fragment.py b/ambuild/ab_fragment.py index 880e94f..f8a82f0 100644 --- a/ambuild/ab_fragment.py +++ b/ambuild/ab_fragment.py @@ -21,71 +21,11 @@ logger = logging.getLogger() -def fragmentFactory( - fragmentType, - filePath=None, - catalyst=False, - markBonded=False, - solvent=False, - static=False, -): - """Dynamically create a fragment object. - - Dynamic classes allow us to have fragments where all fragments of the same type share - class variables such as symbols that are the same between all fragments (to save memory), - but have instance variables for holding attributes like coordinates. - - We need to define a __reduce__ method as otherwise the dynamic classes can't be pickled: - https://stackoverflow.com/questions/11658511/pickling-dynamically-generated-classes - """ - # Make sure fragments don't pollute the namespace - probably needs more thinking about - assert not fragmentType in [ - "fragmentFactory", - "Fragment", - ], f"Unacceptable fragmentType: {fragmentType}" - fobj = type(fragmentType, (Fragment,), dict()) - return fobj( - filePath=filePath, - fragmentType=fragmentType, - solvent=solvent, - markBonded=markBonded, - catalyst=catalyst, - static=static, - ) - - class Fragment(object): """ classdocs """ - #### These class variables are shared by all fragments of a particular type - ## Public variables - catalyst = None - fragmentType = None - markBonded = False - onbondFunction = None # A function to be called when we bond an endGroup - solvent = ( - None # True if this fragment is solvent and should be excluded from clashChecks - ) - static = False - - ## Private variables - _atomTypes = [] - _bodies = [] # a list of which body within this fragment each atom belongs to - _bonds = [] # List of internal fragment bonds - _bonded = [] # List of which atoms are bonded to which - _cellParameters = {} - _charges = [] - _coords = [] - _labels = [] - _masses = [] - _radii = [] - _maxBonds = {} - _radius = -1 - _symbols = [] # ordered array of symbols (in upper case) - _totalMass = -1 - def __init__( self, filePath=None, @@ -98,33 +38,66 @@ def __init__( """ Constructor """ + # This first set of variables are shared by all fragments + # When we copy a fragment the new fragment gets references to the variables + # that were created for the first fragment (see copy) + sharedAttrs = { + "_atomTypes": [], + "_bodies": [], # a list of which body within this fragment each atom belongs to + "_bonds": [], # List of internal fragment bonds + "_bonded": [], # List of which atoms are bonded to which + "catalyst": catalyst, # Whether this block is a catalyst + "_cellParameters": {}, + "fragmentType": fragmentType, + "_labels": [], + "_masses": [], + "markBonded": markBonded, + "onbondFunction": None, # A function to be called when we bond an endGroup + "_radii": [], + "_maxBonds": {}, + "_radius": -1, + "solvent": solvent, # True if this fragment is solvent and should be excluded from clashChecks + "static": static, + "_symbols": [], # ordered array of symbols (in upper case) + "_totalMass": -1, + "_individualAttrs": None, + "_sharedAttrs": None, + } + + # # The variables below here are specific to a fragment and change as the fragment moves - # and is involved in bonds; each fragment gets its own copy of these - self.block = None - self._coords = [] - self._centroid = None - self._centerOfMass = None - self._ext2int = collections.OrderedDict() - self._int2ext = collections.OrderedDict() - self._centerOfMass = None - self._maxAtomRadius = -1 - self._changed = ( - True # Flag for when we've been moved and need to recalculate things - ) - self.blockIdx = None # The index in the list of block data where the data for this fragment starts - self._endGroups = [] # A list of the endGroup objects - self._endGroupBonded = ( - [] - ) # A list of the number of each endGroup that are used in bonds - self.masked = [] # bool - whether the atoms is hidden (e.g. cap or uw atom) - self.unBonded = [] # bool - whether an atom has just been unbonded - - # Init variables - self.catalyst = catalyst - self.fragmentType = fragmentType - self.markBonded = markBonded - self.solvent = solvent - self.static = static + # and is involved in bonds - each fragment gets its own copy of these + individualAttrs = { + "block": None, + "config": None, + "configStr": None, + "_charges": [], + "_coords": [], + "_centroid": None, + "_centerOfMass": None, + "_ext2int": collections.OrderedDict(), + "_int2ext": collections.OrderedDict(), + "_centerOfMass": None, + "_maxAtomRadius": -1, + "_changed": True, # Flag for when we've been moved and need to recalculate things + "blockIdx": None, # The index in the list of block data where the data for this fragment starts + "_endGroups": [], # A list of the endGroup objects + "_endGroupBonded": [], # A list of the number of each endGroup that are used in bonds + "masked": [], # bool - whether the atoms is hidden (e.g. cap or uw atom) + "unBonded": [], # bool - whether an atom has just been unbonded + } + + # Set as attributes of self + for a, v in sharedAttrs.items(): + setattr(self, a, v) + + # Set as attributes of self + for a, v in individualAttrs.items(): + setattr(self, a, v) + + # Set these manually + self._individualAttrs = individualAttrs + self._sharedAttrs = sharedAttrs # Create from the file if filePath: @@ -194,6 +167,9 @@ def bonds(self): bonds = [] for b1, b2 in self._bonds: if not (self.masked[b1] or self.masked[b2]): + # Map to internal and then add blockIdx to get position in block + # b1 = b1 + self.blockIdx + # b2 = b2 + self.blockIdx bonds.append((self._int2ext[b1], self._int2ext[b2])) return bonds @@ -215,7 +191,8 @@ def _calcConfigStr(self): return self.fragmentType + "".join(["1" if c else "0" for c in self.config]) def _calcCenters(self): - """Calculate the center of mass and geometry for this fragment""" + """Calculate the center of mass and geometry for this fragment + """ self._centroid = np.sum(self._coords, axis=0) / np.size(self._coords, axis=0) self._totalMass = np.sum(self._masses) # Centre of mass is sum of the products of the individual coordinates multiplied by the mass, divded by the total mass @@ -225,7 +202,7 @@ def _calcCenters(self): def _calcRadius(self): """ Calculate a simple size metric of the block so that we can screen for whether - two blocks are within touching distance + two _blocks are within touching distance First try a simple approach with a loop just to get a feel for things - Find the largest distance between any atom and the center of geometry @@ -238,9 +215,8 @@ def _calcRadius(self): distances = [xyz_core.distance(self._centroid, coord) for coord in self._coords] imax = np.argmax(distances) dist = distances[imax] - self._radius = ( - dist + self.maxAtomRadius() - ) # Add on the radius of the largest atom + # Add on the radius of the largest atom + self._radius = dist + self.maxAtomRadius() return def _calcProperties(self): @@ -286,10 +262,21 @@ def copy(self): Those attributes in shared are just copied as references as they do not change between fragments of the same fragmentType Those in single are deep-copied as each fragment has its own""" - f = copy.deepcopy(self) - # Update fragment references in the endGroups - for e in f._endGroups: - e.fragment = f + + f = Fragment() + for a in f.__dict__: + if a in self._sharedAttrs.keys(): + setattr(f, a, getattr(self, a)) + elif a in self._individualAttrs.keys(): + setattr(f, a, copy.deepcopy(getattr(self, a))) + else: + # HACKS FOR DEALING WITH OLD FILES + msg = "Missing attribute in fragment copy: {0}".format(a) + logger.critical(msg) + raise RuntimeError(msg) + # Update fragment references in the endGroups + for e in f._endGroups: + e.fragment = f return f def endGroups(self): @@ -300,7 +287,7 @@ def endGroupTypes(self): return set(self._endGroupBonded.keys()) def fillData(self): - """Fill the data arrays from the label""" + """ Fill the data arrays from the label """ self.masked = np.array([False] * len(self._coords)) self.unBonded = [False] * len(self._coords) @@ -318,18 +305,12 @@ def fillData(self): self._maxAtomRadius = np.max(self._radii) return - # Can't use as in some cases the class name doesn't seem to get saved on pickling and - # the fragment has the class 'Fragment' - # @property - # def fragmentType(self): - # return str(self.__class__).split(".")[-1].rstrip("'>") - def freeEndGroups(self): return [eg for eg in self._endGroups if eg.free()] @staticmethod def fromCarFile(carFile): - """ "Abbie did this.""" + """"Abbie did this.""" labels = [] symbols = [] atomTypes = [] @@ -375,7 +356,8 @@ def fromCarFile(carFile): @staticmethod def fromXyzFile(xyzFile): - """ "Jens did this.""" + """"Jens did this. + """ labels = [] symbols = [] atomTypes = [] @@ -548,8 +530,9 @@ def processBodies(self, filepath): basename, suffix = os.path.splitext(filename) bodyFile = os.path.join(dirname, basename + ".ambody") if os.path.isfile(bodyFile): - with open(bodyFile) as f: - self._bodies = np.array([int(l.strip()) for l in f], dtype=int) + self._bodies = np.array( + [int(l.strip()) for l in open(bodyFile)], dtype=int + ) assert len(self._bodies) == len( self._coords ), "Must have as many bodies as coordinates: {0} - {1}!".format( @@ -571,7 +554,8 @@ def totalRadius(self): return self._radius def rotate(self, rotationMatrix, center): - """Rotate the molecule about the given axis by the angle in radians""" + """ Rotate the molecule about the given axis by the angle in radians + """ self._coords = self._coords - center # self._coords = np.array([ np.dot(rotationMatrix, c) for c in self._coords ]) # I don't actually undestand why this works at all... @@ -672,37 +656,6 @@ def setEndGroups(self, endGroupTypes, endGroups, capAtoms, dihedralAtoms, uwAtom assert ( len(eg.intersection(caps)) == 0 ), "Cap atom for one endGroup is another endGroup!" - - self.setupCapTrilateration() - return - - def setupCapTrilateration(self): - NUM_ANCHORS = 4 - caps = set([e.fragmentCapIdx for e in self._endGroups]) - # Get 4 atoms evenly spaced across the list of all of them - this - # should mean we get a good spacing for calculating distances - indexes = [i for i in range(len(self._coords)) if i not in caps] - if len(indexes) < NUM_ANCHORS: - # Cannot do trilateration without 4 non-cap atoms - return - triAtoms = np.round(np.linspace(0, len(indexes) - 1, NUM_ANCHORS)).astype(int) - for eg in self._endGroups: - eg.triAtoms = triAtoms - eg.triDistances = [ - xyz_core.distance(self._coords[eg.fragmentCapIdx], self._coords[a]) - for a in triAtoms - ] - # Check trilateration works for this set of atoms - check_pos = xyz_util.trilaterate3D( - eg.triDistances, [self._coords[p] for p in eg.triAtoms] - ) - - if not np.allclose( - check_pos, self._coords[eg.fragmentCapIdx], atol=1.0e-05 - ): - raise RuntimeError( - f"Incorrect trilateration check position: {check_pos} {self._coords[eg.fragmentCapIdx]}" - ) return def setMaxBond(self, endGroupType, maxBond): @@ -714,7 +667,7 @@ def symbol(self, idxAtom): return self._symbols[self._ext2int[idxAtom]] def translate(self, tvector): - """translate the molecule by the given vector""" + """ translate the molecule by the given vector""" self._coords += tvector self._changed = True return @@ -748,17 +701,6 @@ def updateCharges(self, charges): ) self._charges = np.array(charges) - def __reduce__(self): - """Required to allow dynamically created classes to be pickled. - See: https://docs.python.org/3/library/pickle.html#object.__reduce__ - """ - state = self.__dict__.copy() - return ( - fragmentFactory, - (self.fragmentType,), - state, - ) - def __str__(self): """List the data attributes of this object""" # me = {} diff --git a/ambuild/ab_util.py b/ambuild/ab_util.py index 4c2440f..178592f 100755 --- a/ambuild/ab_util.py +++ b/ambuild/ab_util.py @@ -92,6 +92,20 @@ def fixFragment(fragment): fragment._individualAttrs["configStr"] = fragment.configStr return + if not os.path.isfile(pickleFile): + raise RuntimeError("Cannot find file: {}".format(pickleFile)) + if pickleFile.endswith(GZIP_PKL_SUFFIX): + compressed = True + popen = gzip.open + elif pickleFile.endswith(PKL_SUFFIX): + compressed = False + popen = open + else: + raise RuntimeError( + "Unrecognised pkl file suffix for file {}. Use {} or {}".format( + pickleFile, PKL_SUFFIX, GZIP_PKL_SUFFIX + ) + ) # This is another terrible hack for dealing with old pkl files that used the old class names from ambuild import ab_block from ambuild import ab_bond @@ -106,8 +120,16 @@ def fixFragment(fragment): sys.modules["buildingBlock"] = ab_block sys.modules["cell"] = ab_cell sys.modules["fragment"] = ab_fragment - - myCell = unpickleObj(pickleFile) + mode = "r" + if compressed or PYTHONFLAVOUR == 3: + mode += "b" + # Renamed cell class so need to alias here for old files + with popen(pickleFile, mode) as f: + try: + myCell = pickle.load(f) + except UnicodeDecodeError: + # This probably indicates we are trying to open a filed pickled with Python2 with Python3 + myCell = pickle.load(f, encoding="latin1") del ab_block.Bond del sys.modules["buildingBlock"] del sys.modules["cell"] @@ -163,16 +185,13 @@ def dumpPkl(pickleFile, split=None, nonPeriodic=False, paramsDir=None): if split == "fragments": for t in mycell.fragmentTypes().keys(): data = mycell.cellData(fragmentType=t) - if data: - mycell.writeXyz( - "{0}_{1}_P.xyz".format(prefix, t), data=data, periodic=True - ) - mycell.writeCml( - "{0}_{1}_PV.cml".format(prefix, t), - data=data, - periodic=True, - pruneBonds=True, - ) + mycell.writeXyz("{0}_{1}_P.xyz".format(prefix, t), data=data, periodic=True) + mycell.writeCml( + "{0}_{1}_PV.cml".format(prefix, t), + data=data, + periodic=True, + pruneBonds=True, + ) elif split == "blocks": periodic = True for i, b in enumerate(mycell.blocks.values()): @@ -186,25 +205,21 @@ def dumpPkl(pickleFile, split=None, nonPeriodic=False, paramsDir=None): else: if nonPeriodic: data = mycell.cellData(rigidBody=False, periodic=False) - if data: - mycell.writeCml( - prefix + "_NP.cml", data, periodic=False, pruneBonds=False - ) - mycell.writeXyz(prefix + "_NP.xyz", data=data, periodic=False) - mycell.writeXyz( - prefix + "_NP_types.xyz", data=data, periodic=False, atomTypes=True - ) + mycell.writeCml(prefix + "_NP.cml", data, periodic=False, pruneBonds=False) + mycell.writeXyz(prefix + "_NP.xyz", data=data, periodic=False) + mycell.writeXyz( + prefix + "_NP_types.xyz", data=data, periodic=False, atomTypes=True + ) else: data = mycell.cellData(rigidBody=False) - if data: - mycell.writeXyz(prefix + "_P.xyz", data=data, periodic=True) - mycell.writeXyz( - prefix + "_P_types.xyz", data=data, periodic=True, atomTypes=True - ) - # self.writeCar(prefix+"_P.car",data=data,periodic=True) - mycell.writeCml( - prefix + "_PV.cml", data=data, periodic=True, pruneBonds=True - ) + mycell.writeXyz(prefix + "_P.xyz", data=data, periodic=True) + mycell.writeXyz( + prefix + "_P_types.xyz", data=data, periodic=True, atomTypes=True + ) + # self.writeCar(prefix+"_P.car",data=data,periodic=True) + mycell.writeCml( + prefix + "_PV.cml", data=data, periodic=True, pruneBonds=True + ) return @@ -368,34 +383,6 @@ def is_exe(fpath): return fpath and os.path.exists(fpath) and os.access(fpath, os.X_OK) -def unpickleObj(pickleFile): - if not os.path.isfile(pickleFile): - raise RuntimeError("Cannot find file: {}".format(pickleFile)) - if pickleFile.endswith(GZIP_PKL_SUFFIX): - compressed = True - popen = gzip.open - elif pickleFile.endswith(PKL_SUFFIX): - compressed = False - popen = open - else: - raise RuntimeError( - "Unrecognised pkl file suffix for file {}. Use {} or {}".format( - pickleFile, PKL_SUFFIX, GZIP_PKL_SUFFIX - ) - ) - mode = "r" - if compressed or PYTHONFLAVOUR == 3: - mode += "b" - # Renamed cell class so need to alias here for old files - with popen(pickleFile, mode) as f: - try: - myObj = pickle.load(f) - except UnicodeDecodeError: - # This probably indicates we are trying to open a filed pickled with Python2 with Python3 - myObj = pickle.load(f, encoding="latin1") - return myObj - - if __name__ == "__main__": import argparse diff --git a/ambuild/hoomd2.py b/ambuild/hoomd2.py index 7e72fda..e468ef7 100755 --- a/ambuild/hoomd2.py +++ b/ambuild/hoomd2.py @@ -39,7 +39,6 @@ def __init__(self, paramsDir): self.system = None self.rigidBody = False self.exclusions = set() # particle tags to be ignored in pair-pair interactions - self.idxStaticStart = 0 def checkParameters(self, skipDihedrals=False): assert self.ffield @@ -97,7 +96,6 @@ def createSnapshot(self, data, doCharges=True, doDihedral=True): """ # Reset exclusions here for time being self.exclusions = set() - self.idxStaticStart = 0 # Create snapshot # snap attributes: angles, bonds, box, constraints, dihedrals, impropers, pairs, particles if self.rigidBody: @@ -105,7 +103,7 @@ def createSnapshot(self, data, doCharges=True, doDihedral=True): nRigidParticles = len(data.rigidParticles) nparticles = data.natoms + nRigidParticles rigidCenters = set([r.type for r in data.rigidParticles]) - atomTypes = set(data.atomTypes) # include standard atomTypes in case any static atoms present + atomTypes = set() for r in data.rigidParticles: atomTypes.update(r.b_atomTypes) overlap = atomTypes.intersection(rigidCenters) @@ -208,19 +206,6 @@ def createSnapshot(self, data, doCharges=True, doDihedral=True): rp.b_atomTypes[j] ) idx += 1 - # Finally add in any static particles - self.idxStaticStart = idx - for i in range(len(data.coords)): - if doCharges: - snapshot.particles.charge[idx] = data.charges[i] - snapshot.particles.diameter[idx] = data.diameters[i] - snapshot.particles.image[idx] = data.images[i] - snapshot.particles.mass[idx] = data.masses[i] - snapshot.particles.position[idx] = data.coords[i] - snapshot.particles.typeid[idx] = self.particleTypes.index( - data.atomTypes[i] - ) - idx += 1 else: for i in range(nparticles): if doCharges: @@ -576,7 +561,7 @@ def setupGroups(self, data): # All atoms that are part of static fragments if any(data.static): self.groupStatic = hoomd.group.tag_list( - name="static", tags=[self.idxStaticStart + i for i, s in enumerate(data.static) if s] + name="static", tags=[i for i, s in enumerate(data.static) if s] ) if self.rigidBody: self.groupActive = hoomd.group.difference( @@ -683,28 +668,23 @@ def updateCell(self, cell): box = np.array([self.system.box.Lx, self.system.box.Ly, self.system.box.Lz]) snapshot = self.system.take_snapshot() if self.rigidBody: - idxFreeAtom = len(hoomd.group.rigid_center()) - idxStaticAtom = self.idxStaticStart + atomIdx = len(hoomd.group.rigid_center()) else: - idxFreeAtom = 0 + atomIdx = 0 for block in cell.blocks.values(): - idxBlockAtom = 0 - for frag in block.fragments: - for body in frag.bodies(): - for _ in range(len(body.coords)): - idx = None - if self.rigidBody and body.static: - idx = idxStaticAtom - idxStaticAtom += 1 - else: - idx = idxFreeAtom - idxFreeAtom += 1 - coord = snapshot.particles.position[idx] - coord = xyz_core.unWrapCoord3( - coord, snapshot.particles.image[idx], box, centered=True - ) - block.coord(idxBlockAtom, coord) - idxBlockAtom += 1 + for i in range(block.numAtoms()): + coord = snapshot.particles.position[atomIdx] + coord = xyz_core.unWrapCoord3( + coord, snapshot.particles.image[atomIdx], box, centered=True + ) + block.coord(i, coord) + atomIdx += 1 + if atomIdx != snapshot.particles.N: + raise RuntimeError( + "Read {0} positions but there were {1} particles!".format( + atomIdx, len(self.system.particles) + ) + ) # If we are running (e.g.) an NPT simulation, the cell size may have changed. In this case we need to update # our cell parameters. Repopulate cells will then update the halo cells and add the new blocks diff --git a/ambuild/xyz_core.py b/ambuild/xyz_core.py index 29b30ee..18cecb0 100755 --- a/ambuild/xyz_core.py +++ b/ambuild/xyz_core.py @@ -768,7 +768,8 @@ def orientationQuaternion(mobileCoords, refCoords): def principalMoments(coords, masses): - """http://farside.ph.utexas.edu/teaching/336k/Newtonhtml/node67.html""" + """http://farside.ph.utexas.edu/teaching/336k/Newtonhtml/node67.html + """ I = momentOfInertia(coords, masses) eigval, eigvec = np.linalg.eig(I) return np.sort(eigval) @@ -839,44 +840,18 @@ def rigid_rotate(A, B): """Return rotation matrix to rotate A to B. Assumes both sets of points are already centred. Uses SVD to calculate the rotation. Taken from: http://nghiaho.com/?page_id=671""" - - ## Original Code - # A = np.matrix(A) - # B = np.matrix(B) - # assert len(A) == len(B) - # # dot is matrix multiplication for array - # H = A.T * B - # U, S, Vt = np.linalg.svd(H) - # R = Vt.T * U.T - # # special reflection case - # if np.linalg.det(R) < 0: - # # Reflection detected - # Vt[2, :] *= -1 - # R = Vt.T * U.T - - assert A.shape == B.shape - - # Our matricies are in the opposite order so we need to transpose them - A = A.T - B = B.T - - num_rows, num_cols = A.shape - if num_rows != 3: - raise Exception(f"matrix A is not 3xN, it is {num_rows}x{num_cols}: {A}") - - num_rows, num_cols = B.shape - if num_rows != 3: - raise Exception(f"matrix B is not 3xN, it is {num_rows}x{num_cols}") - - H = A @ np.transpose(B) + A = np.matrix(A) + B = np.matrix(B) + assert len(A) == len(B) + # dot is matrix multiplication for array + H = A.T * B U, S, Vt = np.linalg.svd(H) - R = Vt.T @ U.T + R = Vt.T * U.T # special reflection case if np.linalg.det(R) < 0: # Reflection detected Vt[2, :] *= -1 - R = Vt.T @ U.T - + R = Vt.T * U.T return R @@ -902,12 +877,12 @@ def rotation_matrix(axis, angle): def vecDiff(v1, v2, dim=None, pbc=[True, True, True]): """Difference between vectors with numpy taking PBC into account - This works either with 2 points or a vector of any number of points - Adapted from: http://stackoverflow.com/questions/11108869/optimizing-python-distance-calculation-while-accounting-for-periodic-boundary-co - Changed so that it can cope with distances across more than one cell - Args: - dim - 3 element array with dimensions of unit cell - pbc - 3 element boolean array indicating if this dimension has periodic boundaries + This works either with 2 points or a vector of any number of points + Adapted from: http://stackoverflow.com/questions/11108869/optimizing-python-distance-calculation-while-accounting-for-periodic-boundary-co + Changed so that it can cope with distances across more than one cell + Args: + dim - 3 element array with dimensions of unit cell + pbc - 3 element boolean array indicating if this dimension has periodic boundaries """ # Currently (e.g. cell.dataDict return a list of coordinates rather than a numpy array, so we need to check if we have been given a list # or numpy array and convert accordingly @@ -945,7 +920,7 @@ def vecDiff(v1, v2, dim=None, pbc=[True, True, True]): def vectorAngle(v1, v2): - """Calculate the angle between two vectors + """ Calculate the angle between two vectors Return value in Radians A . B = |A|*|B|*cos(theta) so: theta = arccos( X.Y / |X||Y| ) @@ -968,7 +943,8 @@ def vectorAngle(v1, v2): def unWrapCoord3(coord, image, ldim, centered=False, inplace=False): - """Unwrap a coordinate back into a cell""" + """Unwrap a coordinate back into a cell + """ if not inplace: coord = np.copy(coord) if centered: diff --git a/ambuild/xyz_util.py b/ambuild/xyz_util.py index 0f2ac24..124017e 100755 --- a/ambuild/xyz_util.py +++ b/ambuild/xyz_util.py @@ -24,10 +24,10 @@ class BondLength(object): """Class to hold calculate bond lengths. - + This is required because we use data stored in the ATOM_TYPE_BOND_LENGTHS dict which is calculated from the parameter files, which is read at run time. Therefore - we initialise this object using a parameter file and set it as a module object for + we initialise this object using a parameter file and set it as a module object for use by the bondLength function """ @@ -43,8 +43,8 @@ def __init__(self, bond_param_file): self.ATOM_TYPE_BOND_LENGTHS[p.A][p.B] = float(p.r0) def bondLength(self, atomType1, atomType2): - """Get the characteristic lengths of single bonds as defined in: - Reference: CRC Handbook of Chemistry and Physics, 87th edition, (2006), Sec. 9 p. 46 + """ Get the characteristic lengths of single bonds as defined in: + Reference: CRC Handbook of Chemistry and Physics, 87th edition, (2006), Sec. 9 p. 46 """ # We first see if we can find the bond length in the ATOM_TYPE_BOND_LENGTHS table # If not we fall back to using the bonds calculated from element types @@ -76,7 +76,7 @@ def bondLength(self, atomType1, atomType2): # print "ELEMENT TYPE" return xyz_core.ELEMENT_TYPE_BOND_LENGTHS[symbol2][symbol1] warnings.warn("No data for bond length for %s-%s" % (atomType1, atomType2)) - return -1 + return 1.0 # This needs to be set to the bondLength function of the BondLength class @@ -94,7 +94,8 @@ def __STOP(x, y): def setModuleBondLength(paramFile): """Set the bondLength function of the util module""" global bondLength - bondLength = BondLength(paramFile).bondLength + BL = BondLength(paramFile) + bondLength = BL.bondLength return @@ -126,6 +127,9 @@ def calcBonds( ): """Calculate the bonds for the fragments. This is done at the start when the only coordinates are those in the fragment. + symbols can be chemical elements or atomTypes + If supplied cell is a list/numpy array with the dimensions of the simulation cell, in which case + PBC will be applied """ close = closeAtoms(coords, symbols, dim, maxAtomRadius, boxMargin) v1 = [] @@ -136,28 +140,24 @@ def calcBonds( distances = xyz_core.distance(np.array(v1), np.array(v2), dim=dim) bonds = [] + for i, (idxAtom1, idxAtom2) in enumerate(close): - symbol1 = symbols[idxAtom1] - symbol2 = symbols[idxAtom2] - bond_length = bondLength(symbol1, symbol2) + bond_length = bondLength(symbols[idxAtom1], symbols[idxAtom2]) if bond_length < 0: - warnings.warn( - "calcBonds: no data for bond length for %s-%s - setting to 1.0" - % (symbol1, symbol2) + continue + logger.debug( + "Dist:length {1}({0})-{3}({2}) {4} {5}".format( + symbols[idxAtom1], + idxAtom1, + symbols[idxAtom2], + idxAtom2, + distances[i], + bond_length, ) - bond_length = 1.0 - # logger.debug( - # "Dist:length {1}({0})-{3}({2}) {4} {5}".format( - # symbol1, - # idxAtom1, - # symbol2, - # idxAtom2, - # distances[i], - # bond_length, - # ) - # ) + ) if bond_length - bondMargin < distances[i] < bond_length + bondMargin: bonds.append((idxAtom1, idxAtom2)) + # We sort the bonds on return so that results are deterministic and can be tested return sorted(bonds) @@ -273,21 +273,16 @@ def haloCells(key, boxNum=None, pbc=[True, True, True]): # print "sKey ({},{},{})->({})".format(a,b,c,skey) # cells.add(skey) cells.add(skey) + return list(cells) def label2symbol(name): - """Determine the element type of an atom from its name, e.g. Co_2b -> Co - Returns a capitalised element name + """ Determine the element type of an atom from its name, e.g. Co_2b -> Co + Returns a capitalised element name """ origName = name name = name.strip().upper() - if not name[0].isalpha(): - raise RuntimeError( - "label2symbol first character of name is not a character: {0}".format( - origName - ) - ) # Determine the element from the first 2 chars of the name if len(name) > 2: name = name[0:2] @@ -301,6 +296,12 @@ def label2symbol(name): return name.capitalize() # If it was a valid 2 character symbol we should have picked it up so now only 1 symbol name = name[0] + if not name.isalpha(): + raise RuntimeError( + "label2symbol first character of name is not a character: {0}".format( + origName + ) + ) # Hack - for x return x if name.lower() == "x": return "x" @@ -316,126 +317,6 @@ def label2symbol(name): return -def trilaterate3D(distances, points): - """Taken from: https://github.com/akshayb6/trilateration-in-3d/blob/master/trilateration.py""" - - r1 = distances[0] - r2 = distances[1] - r3 = distances[2] - r4 = distances[3] - p1 = points[0] - p2 = points[1] - p3 = points[2] - p4 = points[3] - e_x = (p2 - p1) / np.linalg.norm(p2 - p1) - i = np.dot(e_x, (p3 - p1)) - e_y = (p3 - p1 - (i * e_x)) / (np.linalg.norm(p3 - p1 - (i * e_x))) - e_z = np.cross(e_x, e_y) - d = np.linalg.norm(p2 - p1) - j = np.dot(e_y, (p3 - p1)) - x = ((r1 ** 2) - (r2 ** 2) + (d ** 2)) / (2 * d) - y = (((r1 ** 2) - (r3 ** 2) + (i ** 2) + (j ** 2)) / (2 * j)) - ((i / j) * (x)) - with np.errstate(invalid="raise"): - z1 = np.sqrt(r1 ** 2 - x ** 2 - y ** 2) - z2 = np.sqrt(r1 ** 2 - x ** 2 - y ** 2) * (-1) - ans1 = p1 + (x * e_x) + (y * e_y) + (z1 * e_z) - ans2 = p1 + (x * e_x) + (y * e_y) + (z2 * e_z) - dist1 = np.linalg.norm(p4 - ans1) - dist2 = np.linalg.norm(p4 - ans2) - if np.abs(r4 - dist1) < np.abs(r4 - dist2): - return ans1 - else: - return ans2 - - -def trilaterate3D_2(distances, points, tolerance=0.0001): - """Taken from: https://math.stackexchange.com/questions/2272223/getting-wrong-coordinate-with-3d-trilateration-in-different-cases-noisy-environ - - See farmuaa6 - """ - r1 = distances[0] - r2 = distances[1] - r3 = distances[2] - r4 = distances[3] - p1 = points[0] - p2 = points[1] - p3 = points[2] - p4 = points[3] - - p1x = p1[0] - p1y = p1[1] - p1z = p1[2] - p2x = p2[0] - p2y = p2[1] - p2z = p2[2] - p3x = p3[0] - p3y = p3[1] - p3z = p3[2] - - x1 = p2x - p1x - y1 = p2y - p1y - z1 = p2z - p1z - n1 = np.sqrt(x1 ** 2 + y1 ** 2 + z1 ** 2) - if n1 <= 0.0: - raise RuntimeError(f"First and second anchor are at the same point: {p1} {p2}") - - x1 = x1 / n1 - y1 = y1 / n1 - z1 = z1 / n1 - - i = (p3x - p1x) * x1 + (p3y - p1y) * y1 + (p3z - p1z) * z1 - - x2 = p3x - p1x - i * x1 - y2 = p3y - p1y - i * y1 - z2 = p3z - p1z - i * z1 - - with np.errstate(invalid="raise"): - n2 = np.sqrt(x2 ** 2 + y2 ** 2 + z2 ** 2) - if n2 <= 0.0: - raise RuntimeError(f"The three anchors are collinear: {p1} {p2} {p3}") - x2 = x2 / n2 - y2 = y2 / n2 - z2 = z2 / n2 - - x3 = y1 * z2 - z1 * y2 - y3 = z1 * x2 - x1 * z2 - z3 = x1 * y2 - y1 * x2 - with np.errstate(invalid="raise"): - n3 = np.sqrt(x3 ** 2 + y3 ** 2 + z3 ** 2) - if n3 <= 0.0: - raise RuntimeError("Something is wrong with anchors") - - x3 = x3 / n3 - y3 = y3 / n3 - z3 = z3 / n3 - - d = np.sqrt((p2x - p1x) ** 2 + (p2y - p1y) ** 2 + (p2z - p1z) ** 2) - j = x2 * (p3x - p1x) + y2 * (p3y - p1y) + z2 * (p3z - p1z) - - xr = (r1 ** 2 - r2 ** 2 + d ** 2) / (2 * d) - yr = (r1 ** 2 - r3 ** 2 + i * i + j * j) / (2 * j) - i * xr / j - with np.errstate(invalid="raise"): - zr = np.sqrt(r1 ** 2 - xr ** 2 - yr ** 2) - - pos1 = [ - p1x + xr * x1 + yr * x2 + zr * x3, - p1y + xr * y1 + yr * y2 + zr * y3, - p1z + xr * z1 + yr * z2 + zr * z3, - ] - pos2 = [ - p1x + xr * x1 + yr * x2 - zr * x3, - p1y + xr * y1 + yr * y2 - zr * y3, - p1z + xr * z1 + yr * z2 - zr * z3, - ] - - dist1 = np.linalg.norm(p4 - pos1) - dist2 = np.linalg.norm(p4 - pos2) - if np.abs(r4 - dist1) < np.abs(r4 - dist2): - return pos1 - else: - return pos2 - - def writeCml( cmlFilename, coords, diff --git a/misc/Restart/Python_header.py b/misc/Restart/Python_header.py new file mode 100755 index 0000000..d05e7f4 --- /dev/null +++ b/misc/Restart/Python_header.py @@ -0,0 +1,27 @@ +#!/usr/bin/env python3 + +# Our imports +from ambuild import ab_cell +from ambuild import ab_util + +import sys, getopt + +inputfile = '' +try: + opts, args = getopt.getopt(sys.argv[1:],"hi:",["ifile="]) +except getopt.GetoptError: + print (sys.argv[0], ' -i ') + sys.exit(2) +for opt, arg in opts: + if opt == '-h': + print (sys.argv[0], ' -i ') + sys.exit(0) + elif opt in ("-i", "--ifile"): + inputfile = arg +sys.argv = [sys.argv[0]] +#print ('Restart file is', inputfile) + +paramsDir='./params' + + #Create Cell and seed it with the blocks +mycell = ab_util.cellFromPickle(inputfile, paramsDir=paramsDir ) diff --git a/misc/Restart/Testrestart b/misc/Restart/Testrestart new file mode 100755 index 0000000..fca4889 --- /dev/null +++ b/misc/Restart/Testrestart @@ -0,0 +1,19 @@ +#!/bin/bash + +# Test whether the code needs a restart + +# This is a demo that simply checks for the phrase 'Error computing cell list' +# in the last line of the output file + +if [ $# -ne 1 ]; then + echo Usage: $0 "[filename]" + exit 1 +fi + +tail -n 1 $1|grep "Error computing cell list" > /dev/null + +if [ $? -eq 0 ]; then + exit 1 +else + exit 0 +fi diff --git a/misc/Restart/run_ambuild_docker_restart.sh b/misc/Restart/run_ambuild_docker_restart.sh new file mode 100755 index 0000000..1545bd2 --- /dev/null +++ b/misc/Restart/run_ambuild_docker_restart.sh @@ -0,0 +1,49 @@ + GNU nano 4.8 run_ambuild_docker_restart.sh +#!/bin/bash + +# Usage: +# ambuild_docker.py script + +# Get root dir and script argumnts +run_dir="$PWD" +ambuild_dir="$( cd "$( dirname "${BASH_SOURCE[0]}" )/.." >/dev/null 2>&1 && pwd )" +script=$1 +extra_args="" +if [ ${#} -ge 2 ]; then + script="${@:(-1):1}" + extra_args="${@:1:$(($#-1))}" +fi + +MAXLOOP=100 + +LOOPCOUNT=0 +RESTART=1 +while [ $RESTART -eq 1 ] && [ $LOOPCOUNT -lt $MAXLOOP ]; do + echo $MAXLOOP + echo Restart file is: + ls -1t step_*.pkl.gz|head -1 + RFILE=`ls -1t step_*.pkl.gz|head -1` + #STDIR=${> restart.o.$LOOPCOUNT 2> restart.e.$LOOPCOUNT} + STDIR=" > restart.o.${LOOPCOUNT} 2> restart.e.${LOOPCOUNT}" + FULLSTRING="${script} -i ${RFILE}" + # Run + docker run \ + --rm \ + --runtime=nvidia \ + --volume $run_dir:$run_dir \ + --volume ${ambuild_dir}/ambuild:/usr/lib/python3/dist-packages/ambuild \ + --workdir $run_dir \ + $extra_args \ + glotzerlab/software:2020.11.18-cuda10 \ + python3 $FULLSTRING > restart.o.${LOOPCOUNT} 2> restart.e.${LOOPCOUNT} + mv runmd.log runmd.log.$LOOPCOUNT + ./testrestart restart.e.$LOOPCOUNT + RESTART=$? + LOOPCOUNT=$[$LOOPCOUNT+1] +done + +if [ $LOOPCOUNT -eq $MAXLOOP ]; then + echo Error: Maximum restart loop count value of $MAXLOOP reached +fi + + diff --git a/misc/run_ambuild_docker.sh b/misc/run_ambuild_docker.sh index 0ff2291..22a6b43 100755 --- a/misc/run_ambuild_docker.sh +++ b/misc/run_ambuild_docker.sh @@ -18,9 +18,8 @@ docker run \ --rm \ --runtime=nvidia \ --volume $run_dir:$run_dir \ ---volume ${ambuild_dir}/ambuild:/usr/lib/python3/dist-packages/ambuild \ ---env PYTHONPATH=/usr/lib/python3/dist-packages \ +--volume ${ambuild_dir}/ambuild:/home/abbie/ambuild/ambuild \ --workdir $run_dir \ $extra_args \ -glotzerlab/software \ +glotzerlab/software:2020.11.18-cuda10 \ python3 $script diff --git a/misc/run_ambuild_util_docker.sh b/misc/run_ambuild_util_docker.sh index 517e312..8d77362 100755 --- a/misc/run_ambuild_util_docker.sh +++ b/misc/run_ambuild_util_docker.sh @@ -6,18 +6,16 @@ # Get root dir and script argumnts ambuild_dir="$( cd "$( dirname "${BASH_SOURCE[0]}" )/.." >/dev/null 2>&1 && pwd )" -# Where python lives in the container -docker_python_install_path=/usr/lib/python3/dist-packages - # Run docker run \ -it \ --rm \ --runtime=nvidia \ --volume "$PWD":/home/glotzerlab \ ---volume ${ambuild_dir}/ambuild:${docker_python_install_path}/ambuild \ ---env PYTHONPATH=${docker_python_install_path} \ +--volume ${ambuild_dir}/ambuild:/usr/lib/python3/dist-packages/ambuild \ --workdir /home/glotzerlab \ -glotzerlab/software \ -python3 ${docker_python_install_path}/ambuild/ab_util.py $* +glotzerlab/software:2020.11.18-cuda10 \ +python3 /usr/lib/python3/dist-packages/ambuild/ab_util.py $* + + diff --git a/tests/blocks/hydrogen.car b/tests/blocks/hydrogen.car old mode 100644 new mode 100755 diff --git a/tests/blocks/hydrogen.csv b/tests/blocks/hydrogen.csv old mode 100644 new mode 100755 diff --git a/tests/clean.sh b/tests/clean.sh index 3eb5257..0a4e3ca 100755 --- a/tests/clean.sh +++ b/tests/clean.sh @@ -1,5 +1,4 @@ #!/bin/bash rm -f *.csv *.log CONFIG CONTROL *.xml *.tsv *.dcd *.gsd *.pyc *.cml *.xyz -rm -rf poreblazer_0 -rm -f *.pkl.gz +#rm -f *.pkl diff --git a/tests/test.py b/tests/test.py new file mode 100755 index 0000000..956bebf --- /dev/null +++ b/tests/test.py @@ -0,0 +1,26 @@ +#!/usr/bin/env python3 +""" +Created on 14 May 2016 + +@author: jmht +""" +import os +import sys +import unittest + + +ROOT_DIR = os.path.abspath(os.path.join(os.path.dirname(__file__), '..')) +sys.path.insert(0, ROOT_DIR) + +TEST_DIR = "." +VERBOSITY = 2 +suite = unittest.TestLoader().discover(TEST_DIR) +if int(suite.countTestCases()) <= 0: + msg = ( + "Could not find any tests to run in directory: {0}".format(TEST_DIR) + + os.linesep + ) + sys.stderr.write(msg) + sys.exit(1) + +unittest.TextTestRunner(verbosity=VERBOSITY, buffer=False).run(suite) diff --git a/tests/testBlock.py b/tests/testBlock.py old mode 100755 new mode 100644 index 1bf1063..92f91fd --- a/tests/testBlock.py +++ b/tests/testBlock.py @@ -1,4 +1,3 @@ -#!/usr/bin/env python3 import logging import math import os @@ -89,15 +88,7 @@ def testCX4(self): self.assertEqual( [0, 0, 0, 0], [e.blockEndGroupIdx for e in cx4_1.freeEndGroups()] ) - self.assertEqual( - [ - 1, - 2, - 3, - 4, - ], - [e.blockCapIdx for e in cx4_1.freeEndGroups()], - ) + self.assertEqual([1, 2, 3, 4,], [e.blockCapIdx for e in cx4_1.freeEndGroups()]) cx4_2 = Block(filePath=self.cx4Car, fragmentType="A") @@ -580,7 +571,7 @@ def testFreeEndGroups(self): self.assertEqual(idxs, ref_idxs) def testMaxBond(self): - f = ab_fragment.fragmentFactory("A", self.ch4_1Car) + f = ab_fragment.Fragment(filePath=self.ch4_1Car, fragmentType="A") f.setMaxBond("A:b", 1) block1 = Block(initFragment=f) block2 = Block(filePath=self.ch4Car, fragmentType="B") @@ -600,7 +591,7 @@ def testMaxBond(self): def testBondingFunction(self): carfile = os.path.join(BLOCKS_DIR, "benzene6.car") - f = ab_fragment.fragmentFactory("A", carfile) + f = ab_fragment.Fragment(filePath=carfile, fragmentType="A") def x(endGroup): fragment = endGroup.fragment @@ -648,7 +639,7 @@ def x(endGroup): def testMultiEndGroups(self): """Test we can move correctly""" # Try with no settings - f = ab_fragment.fragmentFactory("A", self.ch4_1Car) + f = ab_fragment.Fragment(filePath=self.ch4_1Car, fragmentType="A") m1 = Block(initFragment=f) m2 = m1.copy() eg1 = m1.freeEndGroups()[0] @@ -659,7 +650,7 @@ def testMultiEndGroups(self): self.assertEqual(6, len(m1.freeEndGroups())) # Try with specifying bond - f = ab_fragment.fragmentFactory("A", self.ch4_1Car) + f = ab_fragment.Fragment(filePath=self.ch4_1Car, fragmentType="A") m1 = Block(initFragment=f) f.setMaxBond("A:a", 1) m2 = m1.copy() @@ -698,7 +689,9 @@ def testPositionGrowBlock(self): endGroup2 = growBlock.freeEndGroups()[1] # Get position to check - newPos = blockS.newBondPosition(endGroup1, endGroup2) + newPos = blockS.newBondPosition( + endGroup1, growBlock.symbol(endGroup2.blockEndGroupIdx) + ) # Position the block blockS.positionGrowBlock(endGroup1, endGroup2) @@ -726,7 +719,9 @@ def testPositionGrowBlock2(self): endGroup2 = growBlock.freeEndGroups()[0] # Get position to check - newPos = staticBlock.newBondPosition(endGroup1, endGroup2) + newPos = staticBlock.newBondPosition( + endGroup1, growBlock.symbol(endGroup2.blockEndGroupIdx) + ) # staticBlock._symbols.append( 'N' ) # staticBlock._coords.append( newPos ) @@ -747,10 +742,10 @@ def testPositionGrowBlock2(self): return - @unittest.skip("Broken as changed bond lengths so previous check no longer works") def testPositionDihedral(self): staticBlock = Block(filePath=self.benzeneCar, fragmentType="A") + growBlock = staticBlock.copy() growBlock.translateCentroid([3, 4, 5]) @@ -761,7 +756,9 @@ def testPositionDihedral(self): endGroup2 = growBlock.freeEndGroups()[0] # Get position to check - staticBlock.newBondPosition(endGroup1, endGroup2) + staticBlock.newBondPosition( + endGroup1, growBlock.symbol(endGroup2.blockEndGroupIdx) + ) # Position the block staticBlock.positionGrowBlock(endGroup1, endGroup2, dihedral=math.pi / 2) @@ -772,6 +769,8 @@ def testPositionDihedral(self): np.allclose(hcheck, endGroupCoord, rtol=1e-9, atol=1e-7), msg="testCenterOfMass incorrect COM.", ) + + # self.catBlocks( [staticBlock, growBlock ], "both2.xyz") return def testRadius(self): @@ -838,7 +837,6 @@ def testRotate(self): return - @unittest.skip("Broken test - too senstive?") def testSplitFragment(self): """ Test the rotation @@ -879,7 +877,6 @@ def testSplitFragment(self): return - @unittest.skip("Broken test - too sensitive") def testWriteCml(self): """foo""" ch4_1 = Block(filePath=self.benzeneCar, fragmentType="A") diff --git a/tests/testBody.py b/tests/testBody.py index 546b7e0..7e41695 100644 --- a/tests/testBody.py +++ b/tests/testBody.py @@ -24,22 +24,15 @@ class Test(unittest.TestCase): def testBodyConfigStr(self): ch4ca = os.path.join(BLOCKS_DIR, "ch4Ca2.car") ftype = "A" - f1 = ab_fragment.fragmentFactory(ftype, ch4ca) + f1 = ab_fragment.Fragment(filePath=ch4ca, fragmentType=ftype) for i, b in enumerate(f1.bodies()): bstr = "{}{}{}".format(i, ftype, "0000") self.assertEqual(bstr, b.rigidConfigStr) - def testCoords(self): - ch4ca = os.path.join(BLOCKS_DIR, "ch4.car") - ftype = "A" - frag = ab_fragment.fragmentFactory(ftype, ch4ca) - body = list(frag.bodies())[0] - assert np.array_equal(body.coords, frag._coords) - def testPrincipalMoments(self): """Test principal moments don't change with rotation""" ch4 = os.path.join(BLOCKS_DIR, "ch4.car") - frag = ab_fragment.fragmentFactory('A', ch4) + frag = ab_fragment.Fragment(filePath=ch4, fragmentType="A") body = list(frag.bodies())[0] pi1 = body.principalMoments diff --git a/tests/testCell.py b/tests/testCell.py index 8ae67c4..0f27f9c 100644 --- a/tests/testCell.py +++ b/tests/testCell.py @@ -157,8 +157,110 @@ def XtestCapBlocks(self): mycell.blocks[mycell.blocks.keys()[0]] return + @unittest.skipUnless(ab_util.HOOMDVERSION is not None, "Need HOOMD-BLUE to run") + def testCat1Paf2(self): + boxDim = [40, 40, 40] + mycell = Cell(boxDim, paramsDir=PARAMS_DIR) + mycell.libraryAddFragment(filename=self.ch4Car, fragmentType="PAF") + mycell.libraryAddFragment( + filename=self.nh4Car, fragmentType="cat", catalyst=True + ) + mycell.addBondType("PAF:a-cat:a") + mycell.addBondType("PAF:a-PAF:a") + + # Add PAF and grow so that we have multi-PAF blocks + mycell.seed(2, fragmentType="PAF", point=[10.0, 10.0, 10.0], radius=5.0) + mycell.growBlocks( + toGrow=5, cellEndGroups="PAF:a", libraryEndGroups="PAF:a", maxTries=500 + ) + + # Get the ids of the blocks + paf1id, paf2id = mycell.blocks.keys() + paf1 = mycell.blocks[paf1id] + paf2 = mycell.blocks[paf2id] + # Remove from the cell + mycell.delBlock(paf1id) + mycell.delBlock(paf2id) + + # Add catalyst + mycell.seed(1, fragmentType="cat", center=True) + # Now join the pafs to the cat + cat = list(mycell.blocks.values())[0] + paf1eg = paf1.freeEndGroups()[0] + paf2eg = paf2.freeEndGroups()[0] + cat1eg = cat.freeEndGroups()[0] + cat.positionGrowBlock(cat1eg, paf1eg) + bond = Bond(cat1eg, paf1eg) + bond.engage() + + cat2eg = cat.freeEndGroups()[0] + cat.positionGrowBlock(cat2eg, paf2eg) + bond = Bond(cat2eg, paf2eg) + bond.engage() + mycell.newBonds = [bond] # Hack to set newBonds + + self.assertEqual(len(mycell.blocks), 1) + mycell.cat1Paf2(["PAF"], dt=0.00001, optCycles=10000) + self.assertEqual(len(mycell.blocks), 2) + return + + @unittest.skipUnless(ab_util.HOOMDVERSION is not None, "Need HOOMD-BLUE to run") + def testCat2Paf2(self): + """Given two catalysts bonded to each other, each with PAF blocks bonded, break the bond + between the catalysts, move the PAFS from one catalysts to the other, and then join the PAFS + on that catalyst with all the PAFS""" + boxDim = [40, 40, 40] + mycell = Cell(boxDim, paramsDir=PARAMS_DIR) + mycell.libraryAddFragment(filename=self.ch4Car, fragmentType="PAF") + mycell.libraryAddFragment( + filename=self.nh4Car, fragmentType="cat", markBonded=True, catalyst=True + ) + mycell.addBondType("PAF:a-PAF:a") + mycell.addBondType("PAF:a-cat:a") + mycell.addBondType("cat:a*-cat:a*") + + # Add a cat block and bond it to a PAF block + mycell.seed(1, fragmentType="cat", center=True) + mycell.growBlocks( + toGrow=1, cellEndGroups="cat:a", libraryEndGroups="PAF:a", maxTries=500 + ) + # Add three PAF blocks to the PAF + mycell.growBlocks( + toGrow=3, cellEndGroups="PAF:a", libraryEndGroups="PAF:a", maxTries=500 + ) + + # copy the block and get the two endGroups that we will use to position the two catalysts so they can bond + # newblock = self.getLibraryBlock(fragmentType=fragmentType) # Create new block + b1 = list(mycell.blocks.values())[0] + + # Make copy of first block + b2 = b1.copy() + + # Get the two cat* endGroups + endGroup1 = b1.freeEndGroups(endGroupTypes="cat:a*")[0] + endGroup2 = b2.freeEndGroups(endGroupTypes="cat:a*")[0] + + # Position so they can bond + b1.positionGrowBlock(endGroup1, endGroup2) + + # Put b2 in the cell and bond them together + mycell.addBlock(b2) + mycell.checkMove(b2.id) + mycell.processBonds() + # End setup + + self.assertEqual(len(mycell.blocks), 1) + + # Now see if we can split off the two cat blocks and join the two PAF blocks + mycell.cat2Paf2(["PAF"], dt=0.00001, optCycles=10000) + # mycell.dump() + + self.assertEqual(len(mycell.blocks), 3) + return + def testCellIO(self): - """Check we can write out and then read in a cell""" + """Check we can write out and then read in a cell + """ # Remember a coordinate for checking mycell = self.createTestCell() test_coord = mycell.blocks[list(mycell.blocks.keys())[0]].coord(4) @@ -475,9 +577,7 @@ def testDistance(self): dc1 = mycell.distance(nv1, nv2) dn = np.linalg.norm(nv2 - nv1) - self.assertAlmostEqual( - dc1, dn, 11, "Distance within cell:{} | {}".format(dc1, dn) - ) + self.assertEqual(dc1, dn, "Distance within cell:{} | {}".format(dc1, dn)) x = v2[0] + 2 * CELLA y = v2[1] + 2 * CELLB @@ -519,45 +619,23 @@ def testDihedral(self): self.assertEqual(ref, mycell.dihedral(p1, p2, p3, p4)) return + @unittest.skipUnless(ab_util.HOOMDVERSION is not None, "Need HOOMD-BLUE to run") def testDump(self): - """Test we can dump a cell and read it back in correctly""" + """Test we can dump a cell""" boxDim = [30, 30, 30] mycell = Cell(boxDim, paramsDir=PARAMS_DIR) mycell.libraryAddFragment(filename=self.ch4Car, fragmentType="A") - mycell.libraryAddFragment(filename=self.ch4Car, fragmentType="B") - mycell.addBondType("A:a-B:a") + mycell.addBondType("A:a-A:a") - toSeed = 1 + toSeed = 2 + added = mycell.seed(toSeed, center=True, random=False) + self.assertEqual(added, toSeed, "seed") nblocks = 2 - ftype1 = "A" - ftype2 = "B" - ftypes = set([ftype1, ftype2]) - mycell.seed(toSeed, fragmentType=ftype1) - mycell.seed(toSeed, fragmentType=ftype2) - - # Run initial tests - self.assertEqual(len(mycell.blocks), nblocks) - cell_ftypes = set(mycell.fragmentTypes().keys()) - self.assertEqual(ftypes, cell_ftypes) - for b in mycell.blocks.values(): - self.assertEqual(len(b.fragments), 1) - f = b.fragments[0].fragmentType - self.assertIn(f, cell_ftypes) - - # Dump and then read back in - dumpfile = mycell.dump() - mycell = ab_util.cellFromPickle(dumpfile) - - # Run same tests again - self.assertEqual(len(mycell.blocks), nblocks) - cell_ftypes = set(mycell.fragmentTypes().keys()) - self.assertEqual(ftypes, cell_ftypes) - for b in mycell.blocks.values(): - self.assertEqual(len(b.fragments), 1) - f = b.fragments[0].fragmentType - self.assertIn(f, cell_ftypes) - - os.unlink(dumpfile) + added = mycell.growBlocks(nblocks, endGroupType=None, maxTries=1, random=False) + self.assertEqual(added, nblocks, "growBlocks did not return ok") + mycell.optimiseGeometry(rigidBody=True, optCycles=10) + mycell.dump() + os.unlink("step_1" + ab_util.GZIP_PKL_SUFFIX) return def testEndGroupTypes(self): @@ -705,7 +783,7 @@ def testGrowBlocks(self): self.assertEqual(natoms2, (natoms * nblocks) - (nblocks - 1) * 2) block = list(mycell.blocks.values())[0] self.assertTrue( - np.allclose(block.centroid(), [12.98361984, 17.32539607, 11.75529531]) + np.allclose(block.centroid(), [12.91963557, 17.39975016, 11.65120767]) ) self.assertFalse(self.clashes(mycell)) return @@ -795,7 +873,7 @@ def testGrowBlocks2(self): ) self.assertFalse(self.clashes(mycell)) self.assertTrue( - np.allclose(block1.centroid(), [12.76268068, 14.93950509, 14.96662624]) + np.allclose(block1.centroid(), [12.69119085, 14.93774993, 14.96541503]) ) return @@ -1059,7 +1137,8 @@ def testRunMDAll(self): @unittest.skipUnless(ab_util.HOOMDVERSION is not None, "Need HOOMD-BLUE to run") def testOptimiseGeometryDihedral(self): - """ """ + """ + """ mycell = self.createTestCell() mycell.optimiseGeometry(doDihedral=True, quiet=True, optCycles=1000) self.assertFalse(self.clashes(mycell)) @@ -1079,30 +1158,15 @@ def testOptimiseGeometryStatic(self): toGrow=2, cellEndGroups=None, libraryEndGroups=["A:a"], maxTries=10 ) ok = mycell.optimiseGeometry( - rigidBody=False, doDihedral=True, optCycles=1000, dump=False, quiet=True - ) - self.assertFalse(self.clashes(mycell)) - return - - @unittest.skipUnless(ab_util.HOOMDVERSION is not None, "Need HOOMD-BLUE to run") - def testOptimiseGeometryStaticRigid(self): - """Test reading in a static structure defined in the cell""" - mycell = Cell(filePath=self.graphiteCar, paramsDir=PARAMS_DIR) - mycell.libraryAddFragment(filename=self.ch4Car, fragmentType="A") - mycell.addBondType("A:a-A:a") - mycell.seed(3, fragmentType="A", center=True) - mycell.growBlocks( - toGrow=2, cellEndGroups=None, libraryEndGroups=["A:a"], maxTries=10 - ) - ok = mycell.optimiseGeometry( - rigidBody=True, doDihedral=True, optCycles=1000, dump=False, quiet=True + rigidBody=False, doDihedral=True, optCycles=1000, dump=False, quiet=False ) self.assertFalse(self.clashes(mycell)) return @unittest.skipUnless(ab_util.HOOMDVERSION is not None, "Need HOOMD-BLUE to run") def testRunMD(self): - """ """ + """ + """ mycell = self.createTestCell() mycell.runMD( doDihedral=True, @@ -1121,7 +1185,8 @@ def testRunMD(self): @unittest.skipUnless(ab_util.HOOMDVERSION is not None, "Need HOOMD-BLUE to run") def testRunMdNpt(self): - """ """ + """ + """ mycell = self.createTestCell() mycell.runMD( doDihedral=True, @@ -1149,7 +1214,8 @@ def testRunMdNpt(self): "Need HOOMD-BLUE 1 to run", ) def testRunMDAndOptimise(self): - """ """ + """ + """ mycell = self.createTestCell() mycell.runMDAndOptimise(doDihedral=True, quiet=True) self.assertFalse(self.clashes(mycell)) @@ -1409,11 +1475,12 @@ def testSubunit(self): sumall = sum(count) current = [float(c) / float(sumall) for c in count] wratio = [float(r) / float(sum(ratio)) for r in ratio] - self.assertTrue(np.allclose(np.array(wratio), np.array(current), rtol=0.1)) + self.assertTrue(np.allclose(np.array(wratio), np.array(current), rtol=0.05)) return def testSurroundBoxes(self): - """ """ + """ + """ boxDim = [5, 5, 5] mycell = Cell(boxDim, paramsDir=PARAMS_DIR) # box size=1 - need to set manually as not reading in a block @@ -1767,7 +1834,6 @@ def testZipClashPBC2(self): self.assertEqual(made, 1) return - @unittest.skip("Broken test") def testWriteCml(self): """ write out cml diff --git a/tests/testFragment.py b/tests/testFragment.py old mode 100755 new mode 100644 index afc1cff..7b65237 --- a/tests/testFragment.py +++ b/tests/testFragment.py @@ -1,10 +1,13 @@ -#!/usr/bin/env python3 +""" +Created on 4 Mar 2018 + +@author: jmht +""" import os import unittest import context from context import ab_fragment -from context import ab_util from context import xyz_util BLOCKS_DIR = context.BLOCKS_DIR @@ -14,36 +17,10 @@ class Test(unittest.TestCase): - def testFragmentType(self): - ch4 = os.path.join(BLOCKS_DIR, "ch4.car") - ftype = "A" - f1 = ab_fragment.fragmentFactory(ftype, ch4) - self.assertEqual(ftype, f1.fragmentType) - - def testDisallowedFragmentType(self): - ch4 = os.path.join(BLOCKS_DIR, "ch4.car") - ftype = "Fragment" - with self.assertRaises(AssertionError): - f1 = ab_fragment.fragmentFactory(ftype, ch4) - - def testPickleRoundtrip(self): - ch4 = os.path.join(BLOCKS_DIR, "ch4.car") - ftype = "A" - f1 = ab_fragment.fragmentFactory(ftype, ch4) - f2 = f1.copy() - labels = ["A", "A", "A", "A", "A"] - f2._labels = labels - fname = "f1.pkl.gz" - ab_util.pickleObj(f2, fname) - f3 = ab_util.unpickleObj(fname) - self.assertEqual(f3.fragmentType, ftype) - self.assertEqual(f3._labels, labels) - os.unlink(fname) - - def testConfigStr(self): + def testFragmentConfigStr(self): ch4 = os.path.join(BLOCKS_DIR, "ch4.car") ftype = "A" - f1 = ab_fragment.fragmentFactory(ftype, ch4) + f1 = ab_fragment.Fragment(filePath=ch4, fragmentType=ftype) cstr = ftype + "0000" self.assertEqual(cstr, f1.configStr) @@ -56,7 +33,7 @@ def testConfigStr(self): def testMultipleCapAtoms(self): h3car = os.path.join(BLOCKS_DIR, "H3.car") with self.assertRaises(RuntimeError): - ab_fragment.fragmentFactory("A", h3car) + ab_fragment.Fragment(filePath=h3car, fragmentType="A") if __name__ == "__main__": diff --git a/tests/testRigidParticle.py b/tests/testRigidParticle.py index 9908fd3..7eae621 100644 --- a/tests/testRigidParticle.py +++ b/tests/testRigidParticle.py @@ -29,7 +29,7 @@ def testConfig(self): rigidParticleMgr = ab_rigidparticle.RigidParticleManager() ch4ca = os.path.join(BLOCKS_DIR, "ch4Ca2.car") ftype = "A" - f1 = ab_fragment.fragmentFactory(ftype, ch4ca) + f1 = ab_fragment.Fragment(filePath=ch4ca, fragmentType=ftype) rigidParticles = [] for body in f1.bodies(): rigidParticles.append(rigidParticleMgr.createParticle(body)) diff --git a/tests/testXyzCore.py b/tests/testXyzCore.py index c2b16fe..d15cad9 100644 --- a/tests/testXyzCore.py +++ b/tests/testXyzCore.py @@ -110,7 +110,7 @@ def testRigidRotate2(self): return def testRotMat2Quat(self): - M = np.array([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]) + M = np.matrix([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]) Q = xyz_core.quaternion_from_matrix(M) ref = np.array([1.0, 0.0, 0.0, 0.0]) self.assertTrue(np.allclose(Q, ref, rtol=0.0001)) @@ -122,7 +122,7 @@ def testWrapCoord3_1(self): center = True c1u, image = xyz_core.wrapCoord3(c1, dim, center=center) cref = np.array([-4.0, -8.0, -12.0]) - iref = np.array([10, 10, 10], dtype=int) + iref = np.array([10, 10, 10], dtype=np.int) self.assertTrue(np.allclose(c1u, cref)) self.assertTrue(np.array_equal(image, iref)) @@ -132,7 +132,7 @@ def testWrapCoord3_2(self): center = True c1u, image = xyz_core.wrapCoord3(c1, dim, center=center) cref = np.array([4.0, 8.0, 12.0]) - iref = np.array([-11, -11, -11], dtype=int) + iref = np.array([-11, -11, -11], dtype=np.int) self.assertTrue(np.allclose(c1u, cref)) self.assertTrue(np.array_equal(image, iref)) @@ -142,23 +142,23 @@ def testWrapCoord3_list(self): center = True c1u, image = xyz_core.wrapCoord3(c1, dim, center=center) cref = np.array([[-4.0, -8.0, -12.0], [4.0, 8.0, 12.0]]) - iref = np.array([[10, 10, 10], [-11, -11, -11]], dtype=int) + iref = np.array([[10, 10, 10], [-11, -11, -11]], dtype=np.int) self.assertTrue(np.allclose(c1u, cref)) self.assertTrue(np.array_equal(image, iref)) def testUnWrapCoord3_1(self): cin = np.array([-4.0, -8.0, -12.0]) - idxin = np.array([10, 10, 10], dtype=int) + idxin = np.array([10, 10, 10], dtype=np.int) dim = np.array([10, 20, 30]) coord = xyz_core.unWrapCoord3(cin, idxin, dim, centered=True) self.assertTrue(np.allclose(coord, np.array([101.0, 202.0, 303.0]))) def testVectorAngle(self): """Test we can measure angles""" - v1 = np.array([1, 1, 0]) - v2 = np.array([0, 1, 1]) + v1 = np.array([0, 0, 0]) + v2 = np.array([1, 0, 0]) theta = xyz_core.vectorAngle(v1, v2) - self.assertAlmostEqual(60, math.degrees(theta)) + self.assertEqual(180, math.degrees(theta)) return def testVecDiff(self): diff --git a/tests/testXyzUtil.py b/tests/testXyzUtil.py index d307461..7b8a0e4 100644 --- a/tests/testXyzUtil.py +++ b/tests/testXyzUtil.py @@ -5,7 +5,6 @@ # external imports import numpy as np -from context import xyz_core from context import xyz_util from context import PARAMS_DIR @@ -76,50 +75,6 @@ def testCloseAtoms(self): # order of atoms doesn't 'matter self.assertEqual(set(close), set(ref_close)) - def testTrilateration(self): - """ - Test Trilateration to determine atom position - - 3[0,1,0] - 2[-1,0,0] 1[0,0,0] 4[1,0,0] - t[0,-1,0] - """ - - p1 = np.array([0.0, 0.0, 0.0]) - p2 = np.array([-1.0, 0.0, 0.0]) - p3 = np.array([0.0, 1.0, 0.0]) - p4 = np.array([1.0, 0.0, 0.0]) - points = [p1, p2, p3, p4] - - t = np.array([0.0, -1.0, 0.0]) - p1t = xyz_core.distance(t, p1) - p2t = xyz_core.distance(t, p2) - p3t = xyz_core.distance(t, p3) - p4t = xyz_core.distance(t, p4) - distances = [p1t, p2t, p3t, p4t] - - position = xyz_util.trilaterate3D(distances, points) - self.assertTrue(np.allclose(position, t)) - - # - # Real example with Benzene atoms - # - p1 = np.array([12.5632, 6.03464, -1.74782]) - p2 = np.array([13.86152, 4.39379, 0.51081]) - p3 = np.array([14.74447, 6.33331, -2.51304]) - p4 = np.array([11.74412, 4.14417, 1.20077]) - points = [p1, p2, p3, p4] - - t = np.array([12.12731, 6.58007, -2.50777]) - p1t = xyz_core.distance(t, p1) - p2t = xyz_core.distance(t, p2) - p3t = xyz_core.distance(t, p3) - p4t = xyz_core.distance(t, p4) - distances = [p1t, p2t, p3t, p4t] - - position = xyz_util.trilaterate3D(distances, points) - self.assertTrue(np.allclose(position, t)) - if __name__ == "__main__": """