Merge pull request #12 from njzjz/dev2

fix: add test and common fixes

.github/main.workflow

@@ -0,0 +1,25 @@
+workflow "Release to pypi" {
+  on = "release"
+  resolves = ["upload"]
+}
+
+action "check" {
+  uses = "ross/python-actions/setup-py/3.7@627646f618c3c572358bc7bc4fc413beb65fa50f"
+  args = "check"
+}
+
+action "sdist" {
+  uses = "ross/python-actions/setup-py/3.7@627646f618c3c572358bc7bc4fc413beb65fa50f"
+  args = "sdist"
+  needs = "check"
+}
+
+action "upload" {
+  uses = "ross/python-actions/twine@627646f618c3c572358bc7bc4fc413beb65fa50f"
+  args = "upload ./dist/aimdfragmentation-*.tar.gz"
+  secrets = [
+    "TWINE_USERNAME",
+    "TWINE_PASSWORD",
+  ]
+  needs = "sdist"
+}

.gitignore

@@ -101,3 +101,5 @@ ENV/
 .mypy_cache/
 
 .sw*
+.pytest_cache/
+testfiles/

.travis.yml

@@ -0,0 +1,25 @@
+language: c
+os:
+  - linux
+sudo: false
+cache: pip
+
+env:
+  global:
+    - PYTHON_VERSION=3.7
+    - PIP_DEPENDENCIES="tox-conda"
+    - secure: "ryUltJvieQGY43Cd7nmSzKBWEZbywksbPOxJZ7Ig7Y428P8ML9u8zdKd9u2VT6JIdzpcfO0aKLuJfN4x+gL9G6DwJ2zIWDeou7fEcKtTcDc5SmVxhnErTFcgQ+7NTiLS2hufuQVZSpHGQ+GdSaDMa4y0/GRt/yxpKhBMa7dhd5L71hSi0Hs9BH6gJAuNfwaYhMOGy3GVBSvjKKFuH99u5KxaOUJ/R8Dd8evuQh5UgvWT1tpy0boPXCNLOEJbsXZtxqsTwezVoHoFHnuxiIHWABgAVye0DnIwjOCNQyGCey6HGYbKAx+lU7ldwmftEC8Tyody+VCbjz5DScEyxIEOLStzVFD6IFcJ6DgizXAFqZh2CXrBDo2e8IdkWgwZKrfDIh2YQntKc90isPimWamaYxHRh8rLorumWjeSt8D3CjdVoAOinTKLK+W3zK6YFymw1hgfj3h1aGxHBFlcuw4L7CzIc9u4duq5+VUZRQ4qSqPFIddd5N5vE0YBUJKCDm+LIk/YBY9udbM0vJDL4X9sKc02igbY4Mj6nm5TxuYMqHsos0CiDfhUJ8q3JEP6HKXunKIRw/AzF7BayYvVdHeZdSk8g1OhBIy4xi21d9aAs/6UpgTBw8lNVTC+SABtzH0iG574NYI9WmSt5r00qamQSsDnjJjZ5V1eM0LIf3Nk77A="
+
+before_install:
+  - git clone --depth 1 https://github.com/circulosmeos/gdown.pl;
+  - ./gdown.pl/gdown.pl $GAUSSIANURL 'g16.tgz'>/dev/null 2>&1;
+  - tar -xzvf g16.tgz;
+  - export g16root=$(pwd)
+  - source g16/bsd/g16.profile
+
+install:
+  - git clone --depth 1 git://github.com/astropy/ci-helpers.git
+  - source ci-helpers/travis/setup_conda.sh
+
+script:
+  - tox

aimdfragmentation/__init__.py

@@ -0,0 +1,20 @@
+# Copyright 2018 East China Normal University
+"""Init."""
+import logging
+
+import coloredlogs
+from pkg_resources import DistributionNotFound, get_distribution
+
+from .frag import AIMDFragmentation
+
+__all__ = ['AIMDFragmentation']
+
+try:
+    __version__ = get_distribution(__name__).version
+except DistributionNotFound:
+    # package is not installed
+    __version__ = ''
+
+coloredlogs.install(
+    fmt=f'%(asctime)s - AIMDFragmentation {__version__} - %(levelname)s: %(message)s',
+    level=logging.INFO, milliseconds=True)

AIMDFragmentation/frag.py → aimdfragmentation/frag.py

@@ -1,5 +1,5 @@
 #!/usr/bin/env python3
-"""Fragmentation for AIMD"""
+"""Fragmentation for AIMD."""
 
 __author__ = "Jinzhe Zeng"
 __email__ = "[email protected]"
@@ -8,26 +8,35 @@
 
 
 import itertools
+import logging
 import os
 import subprocess as sp
-import sys
 import time
 from collections import Counter, defaultdict
 from functools import partial
 from multiprocessing import Pool, cpu_count
 
 import numpy as np
 import openbabel
-from ase import Atoms
 from ase.geometry import get_distances
 from ase.io import read as readxyz
+
 from gaussianrunner import GaussianAnalyst, GaussianRunner
 
 
 class AIMDFragmentation(object):
-    def __init__(self, nproc_sum=None, nproc=4, cutoff=3.5, xyzfilename="comb.xyz", pdbfilename="comb.pdb", qmmethod="mn15", qmbasis="6-31g(d)", addkw="", qmmem="400MW", atombondnumber=None, logfile=None, outputfile="force.dat", outputenergyfile="energy.dat", unit=1, energyunit=1, pbc=False, cell=[0, 0, 0], gaussian_dir="gaussian_files", command="g16", gaussiancommand=None, jobfile="gaussianjobs", onebodykeyword="scf=(xqc,MaxConventionalCycles=256)", twobodykeyword="scf=(maxcyc=256)", kbodyfile="kbforce.dat", fg=True, kmax=3):
+    def __init__(
+            self, nproc_sum=None, nproc=4, cutoff=3.5, xyzfilename="comb.xyz",
+            pdbfilename=None, qmmethod="mn15", qmbasis="6-31g(d)",
+            addkw="", qmmem="400MW", atombondnumber=None, logfile=None,
+            outputfile="force.dat", outputenergyfile="energy.dat", unit=1,
+            energyunit=1, pbc=False, cell=None, gaussian_dir="gaussian_files",
+            command="g16", gaussiancommand=None, jobfile="gaussianjobs",
+            onebodykeyword="scf=(xqc,MaxConventionalCycles=256)",
+            twobodykeyword="scf=(maxcyc=256)", kbodyfile="kbforce.dat",
+            fg=True, kmax=3):
         self.nproc_sum = nproc_sum if nproc_sum else cpu_count()
-        self.nproc = nproc
+        self.nproc = min(nproc, self.nproc_sum)
         self.cutoff = cutoff
         self.xyzfilename = xyzfilename
         self.qmmethod = qmmethod
@@ -39,7 +48,7 @@ def __init__(self, nproc_sum=None, nproc=4, cutoff=3.5, xyzfilename="comb.xyz",
         self.unit = unit
         self.energyunit = energyunit
         self.pbc = pbc
-        self.cell = cell
+        self.cell = cell if cell else [0, 0, 0]
         self._atomid = {}
         self.jobs = []
         self.gaussian_dir = gaussian_dir
@@ -62,18 +71,16 @@ def run(self):
 
     def _rungaussian(self):
         if not self.gaussiancommand:
-            GaussianRunner(command=self.command, cpu_num=self.nproc_sum, nproc=self.nproc).runGaussianInParallel(
-                'gjf', [os.path.join(self.gaussian_dir, f"{job}.gjf") for job in self.jobs])
+            GaussianRunner(command=self.command, cpu_num=self.nproc_sum,
+                           nproc=self.nproc).runGaussianInParallel('gjf',
+                                                                   [os.path.join(self.gaussian_dir, f"{job}.gjf")
+                                                                    for job in self.jobs])
         else:
             with open(self.jobfile, 'w') as f:
                 print(*[os.path.join(self.gaussian_dir, f"{job}.gjf")
                         for job in self.jobs], file=f)
             sp.call(self.gaussiancommand.split())
 
-    def _logging(self, *message):
-        localtime = time.asctime(time.localtime(time.time()))
-        print(localtime, 'AIMDFragmentation', *message)
-
     def _getjobname(self, *molid):
         molid = sorted(molid)
         return f'{len(molid)}b{"-".join((str(x) for x in molid))}'
@@ -110,8 +117,11 @@ def _printgjf(self, jobname, selected_atomsid):
         buff = []
         # only supported for C, H, and O
         selected_atoms = [self._atoms[atomsid] for atomsid in selected_atomsid]
-        multiplicities = list((3 if atoms.get_chemical_symbols() == ["O", "O"] else (
-            Counter(atoms.get_chemical_symbols())['H'] % 2+1)) for atoms in selected_atoms)
+        multiplicities = list(
+            (3
+             if atoms.get_chemical_symbols() == ["O", "O"]
+             else(Counter(atoms.get_chemical_symbols())['H'] % 2 + 1))
+            for atoms in selected_atoms)
         atoms_whole, multiplicity_whole, kbodykeyword = self._atoms[sum(selected_atomsid, [])], sum(
             multiplicities)-len(selected_atomsid)+1, (self.onebodykeyword if len(selected_atomsid) == 1 else self.twobodykeyword)
         nproc = f'%nproc={self.nproc}'
@@ -135,7 +145,10 @@ def _printgjf(self, jobname, selected_atomsid):
             kw1 = f'# {self.qmmethod}/{self.qmbasis} {kbodykeyword} {self.addkw} guess=fragment={len(selected_atomsid)}'
             kw2 = f'# force {self.qmmethod}/{self.qmbasis} {kbodykeyword} {self.addkw} geom=chk guess=read'
             multiplicities_str = ' '.join(
-                (f'0 {multiplicity}' for multiplicity in itertools.chain((multiplicity_whole,), multiplicities)))
+                (f'0 {multiplicity}'
+                 for multiplicity in itertools.chain(
+                     (multiplicity_whole,),
+                     multiplicities)))
             buff.extend((chk, nproc, mem, kw1, title, multiplicities_str))
             for index, atoms in enumerate(selected_atoms, 1):
                 buff.extend(('{}(Fragment={}) {} {} {}'.format(
@@ -154,14 +167,19 @@ def _processjob(self, molids):
 
     def _printkb(self, k):
         for molids in itertools.combinations(range(1, len(self._mols)+1), k):
-            if all((self._isclose(molida, molidb) for molida, molidb in itertools.combinations(molids, 2))):
+            if all((self._isclose(molida, molidb) for molida,
+                    molidb in itertools.combinations(molids, 2))):
                 self._processjob(molids)
 
     def _isclose(self, molid1, molid2):
         name = f'{molid1}-{molid2}'
-        if not name in self._distances:
-            self._distances[name] = np.min(get_distances(self._atoms[self._mols[molid1-1]].positions, self._atoms[self._mols[molid2-1]].positions,
-                                                         cell=self._atoms.get_cell(), pbc=self._atoms.get_pbc())[1]) <= self.cutoff
+        if name not in self._distances:
+            self._distances[name] = np.min(
+                get_distances(
+                    self._atoms[self._mols[molid1 - 1]].positions, self.
+                    _atoms[self._mols[molid2 - 1]].positions,
+                    cell=self._atoms.get_cell(),
+                    pbc=self._atoms.get_pbc())[1]) <= self.cutoff
         return self._distances[name]
 
     def _readbond(self):
@@ -179,7 +197,7 @@ def _readforce(self, mols):
             os.path.join(self.gaussian_dir, f'{jobname}.log'))
         qm_force = results['force']
         qm_energy = results['energy']
-        if not qm_force is None:
+        if qm_force is not None:
             forces = np.zeros((self._natom, 3))
             forces[self._atomid[jobname]] = qm_force
             forces *= self.unit
@@ -194,10 +212,10 @@ def fold(self):
                 loadingfold = np.loadtxt(self.kbodyfile)
                 if loadingfold.shape == (self._natom, 3*self.kmax):
                     self._fold = loadingfold
-                    self._logging("Load old forces.")
+                    logging.info("Load old forces.")
         if self._fold is None:
             self._fold = np.zeros((self._natom, 3*self.kmax))
-            self._logging("No old forces found. Use 0 instead.")
+            logging.warn("No old forces found. Use 0 instead.")
         return self._fold
 
     def _takeforce(self):
@@ -207,21 +225,35 @@ def _takeforce(self):
         molsforces = defaultdict(partial(np.zeros, (self._natom, 3)))
         molsenergies = defaultdict(float)
         with Pool(self.nproc_sum) as pool:
-            kbodyresults = [pool.imap(
-                self._readforce, [molids for molids in itertools.combinations(range(1, len(self._mols)+1), i+1) if self._getjobname(*molids) in self.jobs]) for i in range(self.kmax)]
+            kbodyresults = [
+                pool.imap(
+                    self._readforce,
+                    [molids
+                     for molids in itertools.combinations(
+                         range(1, len(self._mols) + 1),
+                         i + 1) if self._getjobname(*molids) in self.jobs])
+                for i in range(self.kmax)]
             for i, results in enumerate(kbodyresults):
                 for force, energy, mols in results:
-                    if not force is None:
-                        molsforces[mols] = force - np.sum((np.sum(
-                            (molsforces[klessmols] for klessmols in itertools.combinations(mols, j)), axis=0) for j in range(i+1)), axis=0)
-                        molsenergies[mols] = energy - np.sum((np.sum(
-                            (molsenergies[klessmols] for klessmols in itertools.combinations(mols, j))) for j in range(i+1)))
+                    if force is not None:
+                        molsforces[mols] = force - np.sum(tuple(
+                            np.sum(tuple(
+                                molsforces[klessmols]
+                                 for klessmols in itertools.combinations(
+                                    mols, j)),
+                                axis=0) for j in range(i + 1)),
+                            axis=0)
+                        molsenergies[mols] = energy - np.sum(np.fromiter(
+                            (np.sum(np.fromiter(
+                                (molsenergies[klessmols]
+                                 for klessmols in itertools.combinations(
+                                    mols, j)), float)) for j in range(i + 1)), float))
                         kbodyforces[i] += molsforces[mols]
                         kbodyenergies[i] += molsenergies[mols]
                     else:
                         jobname = self._getjobname(*mols)
-                        self._logging(
-                            f'WARNING: No forces of {jobname} found. Use the old forces instead.')
+                        logging.warn(
+                            f'No forces of {jobname} found. Use the old forces instead.')
                         self.errorfiles.append(os.path.join(
                             self.gaussian_dir, f'{jobname}.log'))
                         if i == 0:
@@ -233,8 +265,8 @@ def _takeforce(self):
                 if kbodyerroratoms[i]:
                     kbodyforces[i][kbodyerroratoms[i]
                                    ] = self.fold[kbodyerroratoms[i]][:, 3:6]
-                    self._logging("Atom", *kbodyerroratoms[i],
-                                  f"use(s) the old {i}-body forces.")
+                    logging.info(
+                        f"Atom {' '.join(kbodyerroratoms[i])} use(s) the old {i}-body forces.")
         finalforces = np.sum(kbodyforces, axis=0)
         finalenergies = np.sum(kbodyenergies)
         # Make the resultant force equal to 0
@@ -244,10 +276,10 @@ def _takeforce(self):
         np.savetxt(self.outputfile, finalforces, fmt='%16.9f')
         with open(self.outputenergyfile, 'w') as f:
             f.write(f'{finalenergies:16.9f}')
-        np.savetxt(self.kbodyfile, np.hstack(
-            (kbodyforces[i] for i in range(self.kmax))), fmt='%16.9f')
+        np.savetxt(self.kbodyfile, np.hstack(tuple(
+            kbodyforces[i] for i in range(self.kmax))), fmt='%16.9f')
         forcesum = np.sum(finalforces, axis=0)
         forcesumdis = np.linalg.norm(forcesum)
-        self._logging(f"Energy: {finalenergies:16.9f}")
-        self._logging("Resultant force:", *(f"{x:16.9f}" for x in forcesum))
-        self._logging(f"Magnitude: {forcesumdis:16.9f}")
+        logging.info(f"Energy: {finalenergies:16.9f}")
+        logging.info("Resultant force: {:16.9f} {:16.9f} {:16.9f}".format(*forcesum))
+        logging.info(f"Magnitude: {forcesumdis:16.9f}")

aimdfragmentation/test/__init__.py

@@ -0,0 +1 @@
+"""Test."""

aimdfragmentation/test/test.xyz

@@ -0,0 +1,14 @@
+12
+Only for test.
+C 7.512907 12.696144 1.222267
+H 7.293064 13.280151 2.138578
+H 6.574419 12.213078 0.885761
+H 7.898877 13.370358 0.431661
+H 8.280431 11.922762 1.429509
+C 7.574584 15.462496 1.102391
+H 7.385296 15.849984 2.124015
+H 7.986105 16.278961 0.474207
+H 8.302545 14.626973 1.150956
+H 6.625172 15.092911 0.663442
+O 7.183741 9.090479 0.638752
+O 6.104619 9.359606 0.936382

aimdfragmentation/test/test_fragment.py

@@ -0,0 +1,25 @@
+import os
+import unittest
+import logging
+
+import pkg_resources
+
+from aimdfragmentation import AIMDFragmentation
+
+
+class Test_all(unittest.TestCase):
+    def test_fragment(self):
+        xyzfilename = 'test_fragment.xyz'
+        with open(xyzfilename, 'w') as f:
+            f.write(pkg_resources.resource_string(
+                __name__, 'test.xyz').decode())
+        af = AIMDFragmentation(cutoff=6.0, xyzfilename=xyzfilename)
+        af.run()
+        for filename in (af.outputenergyfile, af.outputfile):
+            with open(filename) as f:
+                logging.info(filename)
+                print(f.read())
+
+
+if __name__ == '__main__':
+    unittest.main()