Repository: cabb99/open3spn2 (language: python)

Issue:
Adds analysis script



Discussion / hints:

Fix patch:
diff --git a/.gitattributes b/.gitattributes
new file mode 100644
index 0000000..f89001e
--- /dev/null
+++ b/.gitattributes
@@ -0,0 +1,1 @@
+open3SPN2/_version.py export-subst
diff --git a/.github/workflows/main.yaml b/.github/workflows/main.yaml
new file mode 100644
index 0000000..d3574b3
--- /dev/null
+++ b/.github/workflows/main.yaml
@@ -0,0 +1,56 @@
+name: CI
+
+on:
+  push:
+    branches:
+      - dev
+      - main
+  pull_request:
+    branches:
+      - dev
+      - main
+
+jobs:
+  test:
+
+    runs-on: ubuntu-latest
+
+    strategy:
+      matrix:
+        python-version: ["3.11"]  # Note: Make sure all dependencies are compatible with this Python version
+
+    steps:
+    - uses: actions/checkout@v4
+
+    - name: Set up Conda
+      uses: conda-incubator/setup-miniconda@v2
+      with:
+        miniconda-version: "latest"
+        channels: conda-forge,defaults
+
+    - name: Create and activate Conda environment
+      shell: bash -l {0}
+      run: |
+        conda env create -f environment.yml
+        
+
+    - name: Install package
+      shell: bash -l {0}
+      run: |
+        conda activate open3spn2
+        pip install -e .
+
+    - name: Run tests with pytest
+      shell: bash -l {0}
+      run: |
+        conda activate open3spn2
+        pytest --skip-platform=opencl Tests/
+
+    - name: Slack Notification (on failure)
+      if: failure()
+      uses: act10ns/slack@v2
+      with:
+        status: ${{ job.status }}
+        steps: ${{ toJson(steps) }}
+      env:
+        SLACK_WEBHOOK_URL: ${{ secrets.SLACK_WEBHOOK_URL }}
\ No newline at end of file
diff --git a/.gitignore b/.gitignore
index 0b3c0d1..e841d22 100644
--- a/.gitignore
+++ b/.gitignore
@@ -119,3 +119,10 @@ dmypy.json
 
 #Scratch files
 *Scratch.ipynb
+
+#IDE
+.vscode
+.idea
+
+#Benchmark tests
+Tests/Benchmark
diff --git a/.travis.yml b/.travis.yml
deleted file mode 100644
index 1e900d8..0000000
--- a/.travis.yml
+++ /dev/null
@@ -1,55 +0,0 @@
-language: python
-python:
-  - "3.4"
-  - "3.5"
-  - "3.6"
-  - "3.7"
-  - "3.8"
-  - "3.8-dev"
-  - "nightly"  # nightly build
-env:
-  OPENCL=true
-
-matrix:
-  include:
-    - sudo: required
-      dist: xenial
-      env: ==CPU_OPENCL==
-           OPENCL=true
-           CUDA=false
-      addons: {apt: {packages: []}}
-
-
-before_install:
-  - wget http://s3.amazonaws.com/omnia-ci/AMD-APP-SDKInstaller-v3.0.130.135-GA-linux64.tar.bz2;
-  - tar -xjf AMD-APP-SDK*.tar.bz2;
-  - AMDAPPSDK=${HOME}/AMDAPPSDK;
-  - export OPENCL_VENDOR_PATH=${AMDAPPSDK}/etc/OpenCL/vendors;
-  - mkdir -p ${OPENCL_VENDOR_PATH};
-  - sh AMD-APP-SDK*.sh --tar -xf -C ${AMDAPPSDK};
-  - echo libamdocl64.so > ${OPENCL_VENDOR_PATH}/amdocl64.icd;
-  - export LD_LIBRARY_PATH=${AMDAPPSDK}/lib/x86_64:${LD_LIBRARY_PATH};
-  - chmod +x ${AMDAPPSDK}/bin/x86_64/clinfo;
-  - ${AMDAPPSDK}/bin/x86_64/clinfo;
-  - sudo apt-get install -y libgl1-mesa-dev;
-
-
-install:
-- wget https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh -O miniconda.sh;
-- bash miniconda.sh -b -p $HOME/miniconda
-- export PATH="$HOME/miniconda/bin:$PATH"
-- hash -r
-- conda config --set always_yes yes --set changeps1 no
-- conda install -c anaconda setuptools
-- conda update -q conda
-- conda info -a
-- conda env create -f environment.yaml
-- source activate py36
-script:
-- conda list
-- conda env list
-- python -c "help('modules')"
-- nosetests open3SPN2/ff3SPN2.py
-notifications:
-  slack:
-    secure: AI0yyURG2eMEKuQResaLfVwj0G+Rb49MIw/PNacgOJqUiYJ7bEPWbSuQ61ItInKHhV3KWCtFSsgRxp96RYghSOcHwH8TGlVDQ4ZzjheFP+DGQRN3ORQsr/jq2QhjTwcADi4otUqRKzeVMq+WHaGsXq2mahvlcXkuhjTk0h1vH+9B+/v+IF1qZfdKOjhrGaeu2vwZ4m+eIuigHLZICC6kXBxZ7LVBfrJt2+lcdwgnnpVOdjNwYPLXnDBrkw/qqrX/gHkklUgEH89ZuhSIh2fwGPHi1mIgxTxiiwOpNSTTA6mEXEv98Z4VBf1P7U22qMCRwvcS9ESZHsE5wyZqxFgNGjp9xZJliSGtc8X2aBNrOOwSCEY6LusWSgJtYLYzbgJhQL526qVyS8Rc3hgh4P3JfJT5TtvSA4SmSkErVtXOOrej4ipYbOJf1DvKE/VjAVKurX/OXi0obmkpwrKhbYR3nufIPchQMDjElh3At2/jOB256TMpHcKMpwgae4gwagkkZ9HjanVOwW/B3FVM3ENEliyVEJtdUmko//66NvrsFRU4hzAHANzJulxQ1xO2JyDANBlQCG684V3TPkYYcImAPspzsbzA6fbiup1eeqUCUBY8/jK/2ZimXwDm8CLpOYuPDqnB9L9klKX2Xn1LR2MfhGYp2lWMu5EuIn0O61rcUTc=
diff --git a/README.md b/README.md
index 3fbf0f0..8567dc6 100644
--- a/README.md
+++ b/README.md
@@ -1,4 +1,4 @@
-[![Build Status](https://travis-ci.org/cabb99/open3spn2.svg?branch=master)](https://travis-ci.org/cabb99/open3spn2?branch=master)
+[![Actions Status](https://github.com/cabb99/open3spn2/workflows/build/badge.svg)](https://github.com/cabb99/open3spn2/actions)
 [![Documentation Status](https://readthedocs.org/projects/open3spn2/badge/?version=latest)](https://open3spn2.readthedocs.io/en/latest/?badge=latest)
 [![Anaconda-Server Badge](https://anaconda.org/wolynes-lab/open3spn2/badges/installer/conda.svg)](https://conda.anaconda.org/wolynes-lab)
 
@@ -11,20 +11,14 @@ A Implementation of the 3SPN.2 and 3SPN.2C coarse-grained molecular model of DNA
 
 Installation of the open3SPN2 repository is available through anaconda. Some of the dependencies  are contained in the omnia channel (openmm, pdbfixer) or in the conda-forge channel (mdtraj).
 
-```conda config --append channels omnia```
-
 ```conda config --append channels conda-forge```
 
 ```conda install -c wolynes-lab open3spn2```
 
-It is also necessary to install [X3DNA](http://x3dna.org/) >= 2.4 and set the environment variable `$X3DNA` to the location of the installation folder.For protein-DNA simulations you will also need to install [openAWSEM](https://github.com/npschafer/openawsem) and add the openAWSEM path to the `$PYTHONPATH` environment variable. In linux you can set the path variables on `~/.bashrc`.
-
-```export X3DNA=/path/to/x3dna-v2.4```
+For protein-DNA simulations you will also need to install [openAWSEM](https://github.com/npschafer/openawsem) and add the openAWSEM path to the `$PYTHONPATH` environment variable. In linux you can set the path variables on `~/.bashrc`.
 
 ```export PYTHONPATH=/path/to/openAWSEM:$PYTHONPATH```
 
-Note: open3SPN2 requires the installation of openMM. As of Feb 10 2020, openMM requires 3.5 <= python <= 3.7 and 7.5 <= CUDA <= 10.1.
-
 ## Documentation
 
 Further documentation and tutorials are hosted in [readthedocs](https://open3spn2.readthedocs.io/en/latest/).
diff --git a/docs/source/tutorial.rst b/docs/source/tutorial.rst
index d02dec9..03a1528 100644
--- a/docs/source/tutorial.rst
+++ b/docs/source/tutorial.rst
@@ -63,8 +63,8 @@ To set up the simulation you will need the openmm package as detailed in :doc:`i
 
 .. testcode::
 
-    import simtk.openmm
-    import simtk.openmm.app
+    import openmm
+    import openmm.app
     import simtk.unit
     import sys
     import numpy as np
@@ -76,7 +76,7 @@ To set up the simulation you will need the openmm package as detailed in :doc:`i
     #Set initial positions
     simulation.context.setPositions(s.coord.getPositions())
     
-    energy_unit=simtk.openmm.unit.kilojoule_per_mole
+    energy_unit=openmm.unit.kilojoule_per_mole
     #Total energy
     state = simulation.context.getState(getEnergy=True)
     energy = state.getPotentialEnergy().value_in_unit(energy_unit)
@@ -109,8 +109,8 @@ Please make sure that the energies obtained coincide with the energies shown her
 .. testcode::
     
     #Add simulation reporters
-    dcd_reporter=simtk.openmm.app.DCDReporter(f'output.dcd', 1000)
-    energy_reporter=simtk.openmm.app.StateDataReporter(sys.stdout, 1000, step=True,time=True,
+    dcd_reporter=openmm.app.DCDReporter(f'output.dcd', 1000)
+    energy_reporter=openmm.app.StateDataReporter(sys.stdout, 1000, step=True,time=True,
                                                        potentialEnergy=True, temperature=True)
     simulation.reporters.append(dcd_reporter)
     simulation.reporters.append(energy_reporter)
@@ -208,10 +208,10 @@ You can find this example on the `examples/Protein_DNA <https://github.com/cabb9
 
     #Create the merged system
 
-    pdb=simtk.openmm.app.PDBFile('clean.pdb')
+    pdb=openmm.app.PDBFile('clean.pdb')
     top=pdb.topology
     coord=pdb.positions
-    forcefield=simtk.openmm.app.ForceField(ffAWSEM.xml,open3SPN2.xml)
+    forcefield=openmm.app.ForceField(ffAWSEM.xml,open3SPN2.xml)
     s=forcefield.createSystem(top)
 
 .. code:: ipython3
@@ -232,7 +232,7 @@ You can find this example on the `examples/Protein_DNA <https://github.com/cabb9
     protein=ffAWSEM.Protein.fromCoarsePDB('clean.pdb',sequence=protein_sequence_one)
     dna.periodic=False
     protein.periodic=False
-    
+
     #Copy the AWSEM parameter files
     ffAWSEM.copy_parameter_files()
 
@@ -242,7 +242,7 @@ You can find this example on the `examples/Protein_DNA <https://github.com/cabb9
     keepCMMotionRemover=True
     j=0
     for i, f in enumerate(s.getForces()):
-        if keepCMMotionRemover and i == 0 and f.__class__ == simtk.openmm.CMMotionRemover:
+        if keepCMMotionRemover and i == 0 and f.__class__ == openmm.CMMotionRemover:
             # print('Kept ', f.__class__)
             j += 1
             continue
@@ -253,10 +253,10 @@ You can find this example on the `examples/Protein_DNA <https://github.com/cabb9
         assert len(s.getForces()) == 0, 'Not all the forces were removed'
     else:
         assert len(s.getForces()) <= 1, 'Not all the forces were removed'
-    
+
 .. code:: ipython3
 
-    #Initialize the force dictionary    
+    #Initialize the force dictionary
     forces={}
     for i in range(s.getNumForces()):
         force = s.getForce(i)
@@ -276,14 +276,15 @@ You can find this example on the `examples/Protein_DNA <https://github.com/cabb9
     openAWSEMforces = dict(Connectivity=ffAWSEM.functionTerms.basicTerms.con_term,
                            Chain=ffAWSEM.functionTerms.basicTerms.chain_term,
                            Chi=ffAWSEM.functionTerms.basicTerms.chi_term,
-                           Excl=ffAWSEM.functionTerms.basicTerms.excl_term,
+                           Excl=ffAWSEM.functionTerms.basicTerms.excl_term_v2,
                            rama=ffAWSEM.functionTerms.basicTerms.rama_term,
                            rama_pro=ffAWSEM.functionTerms.basicTerms.rama_proline_term,
+                           #rama_ss=ffAWSEM.functionTerms.basicTerms.rama_ssweight_term,
                            contact=ffAWSEM.functionTerms.contactTerms.contact_term,
-                           frag  = partial(ffAWSEM.functionTerms.templateTerms.fragment_memory_term, 
-                                           frag_file_list_file="./single_frags.mem", 
-                                           npy_frag_table="./single_frags.npy", 
-                                           UseSavedFragTable=False, 
+                           frag  = partial(ffAWSEM.functionTerms.templateTerms.fragment_memory_term,
+                                           frag_file_list_file="./single_frags.mem",
+                                           npy_frag_table="./single_frags.npy",
+                                           UseSavedFragTable=False,
                                            k_fm=0.04184/3),
                            beta1 = ffAWSEM.functionTerms.hydrogenBondTerms.beta_term_1,
                            beta2 = ffAWSEM.functionTerms.hydrogenBondTerms.beta_term_2,
@@ -348,8 +349,8 @@ You can find this example on the `examples/Protein_DNA <https://github.com/cabb9
     frag
     Loading Fragment files(Gro files)
     Saving fragment table as npy file to speed up future calculation.
-    All gro files information have been stored in the ./single_frags.npy.             
-    You might want to set the 'UseSavedFragTable'=True to speed up the loading next time.             
+    All gro files information have been stored in the ./single_frags.npy.
+    You might want to set the 'UseSavedFragTable'=True to speed up the loading next time.
     But be sure to remove the .npy file if you modify the .mem file. otherwise it will keep using the old frag memeory.
     beta1
     beta_1 term ON
@@ -365,17 +366,19 @@ You can find this example on the `examples/Protein_DNA <https://github.com/cabb9
     No ssweight given, assume all zero
 
 
+Then set-up the simulation
+
 .. code:: ipython3
 
     # Set up the simulation
-    temperature=300 * simtk.openmm.unit.kelvin
+    temperature=300 * openmm.unit.kelvin
     platform_name='OpenCL' #'Reference','CPU','CUDA', 'OpenCL'
 
-    integrator = simtk.openmm.LangevinIntegrator(temperature, 1 / simtk.openmm.unit.picosecond, 2 * simtk.openmm.unit.femtoseconds)
-    platform = simtk.openmm.Platform.getPlatformByName(platform_name)
-    simulation = simtk.openmm.app.Simulation(top,s, integrator, platform)
+    integrator = openmm.LangevinIntegrator(temperature, 1 / openmm.unit.picosecond, 2 * openmm.unit.femtoseconds)
+    platform = openmm.Platform.getPlatformByName(platform_name)
+    simulation = openmm.app.Simulation(top,s, integrator, platform)
     simulation.context.setPositions(coord)
-    energy_unit=simtk.openmm.unit.kilojoule_per_mole
+    energy_unit=openmm.unit.kilojoule_per_mole
     state = simulation.context.getState(getEnergy=True)
     energy = state.getPotentialEnergy().value_in_unit(energy_unit)
     print(energy)
@@ -389,7 +392,7 @@ You can find this example on the `examples/Protein_DNA <https://github.com/cabb9
 
     #Obtain total energy
 
-    energy_unit=simtk.openmm.unit.kilojoule_per_mole
+    energy_unit=openmm.unit.kilojoule_per_mole
     state = simulation.context.getState(getEnergy=True)
     energy = state.getPotentialEnergy().value_in_unit(energy_unit)
     print('TotalEnergy',round(energy,6),energy_unit.get_symbol())
@@ -464,7 +467,7 @@ You can find this example on the `examples/Protein_DNA <https://github.com/cabb9
     8000,16.00000000000201,-3296.966064453125,305.38198987465074
     9000,18.000000000000902,-3182.54443359375,316.07187422798313
     10000,19.999999999999794,-3229.941650390625,309.6002450725328
-        
+
 .. code:: ipython3
 
     #Get the detailed energy after the simulation
diff --git a/environment.yaml b/environment.yml
similarity index 82%
rename from environment.yaml
rename to environment.yml
index 2d7a512..afd77bd 100644
--- a/environment.yaml
+++ b/environment.yml
@@ -1,6 +1,5 @@
-name: py36
+name: open3spn2
 channels:
-- omnia
 - conda-forge
 - defaults
 dependencies:
@@ -9,8 +8,8 @@ dependencies:
 - pandas
 - numpy
 - scipy
-- fftw3f
 - mdtraj
 - openmm
 - pdbfixer
 - nose
+- pytest
diff --git a/examples/from_sequence/Simulation_from_sequence.ipynb b/examples/from_sequence/Simulation_from_sequence.ipynb
index a7905f5..1558edb 100644
--- a/examples/from_sequence/Simulation_from_sequence.ipynb
+++ b/examples/from_sequence/Simulation_from_sequence.ipynb
@@ -19,7 +19,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": 7,
    "metadata": {},
    "outputs": [
     {
@@ -66,9 +66,17 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": 8,
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "Warning: importing 'simtk.openmm' is deprecated.  Import 'openmm' instead.\n"
+     ]
+    }
+   ],
    "source": [
     "import sys\n",
     "import simtk.openmm\n",
@@ -90,7 +98,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": 9,
    "metadata": {},
    "outputs": [
     {
@@ -101,7 +109,7 @@
       "Bond 2.1e-05 kJ/mol\n",
       "Angle 0.002319 kJ/mol\n",
       "Stacking -793.345703 kJ/mol\n",
-      "Dihedral -839.999756 kJ/mol\n",
+      "Dihedral -839.999817 kJ/mol\n",
       "BasePair -522.301514 kJ/mol\n",
       "CrossStacking -138.427765 kJ/mol\n",
       "Exclusion 0.219572 kJ/mol\n",
@@ -131,7 +139,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": 10,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -144,6 +152,27 @@
     "simulation.reporters.append(energy_reporter)"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": 15,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "system_xml=simtk.openmm.XmlSerializer.serialize(s)\n",
+    "with open('temp.xml','w+') as f:\n",
+    "    f.write(system_xml)\n",
+    "    \n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "simtk.openmm.System."
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": 5,
@@ -177,9 +206,9 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3",
+   "display_name": "Python [conda env:open3spn2]",
    "language": "python",
-   "name": "python3"
+   "name": "conda-env-open3spn2-py"
   },
   "language_info": {
    "codemirror_mode": {
@@ -191,7 +220,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.10"
+   "version": "3.6.15"
   },
   "toc": {
    "base_numbering": 1,
diff --git a/examples/from_sequence/output.dcd b/examples/from_sequence/output.dcd
index 8e66828..e69de29 100644
Binary files a/examples/from_sequence/output.dcd and b/examples/from_sequence/output.dcd differ
diff --git a/examples/from_sequence/temp_template.pdb b/examples/from_sequence/temp_template.pdb
index 6b78b6a..8a6c933 100644
--- a/examples/from_sequence/temp_template.pdb
+++ b/examples/from_sequence/temp_template.pdb
@@ -1,4 +1,4 @@
-REMARK    3DNA v2.3.1-2017jun24, created and maintained by Xiang-Jun Lu (PhD)
+REMARK    3DNA v2.4.6-2021nov01, created and maintained by xiangjun@x3dna.org
 ATOM      1  P    DA A   1      -0.193   9.323  -1.389  1.00  1.00           P  
 ATOM      2  OP1  DA A   1      -0.207  10.666  -2.008  1.00  1.00           O  
 ATOM      3  OP2  DA A   1       0.673   9.164  -0.199  1.00  1.00           O  
diff --git a/meta.yaml b/meta.yaml
index afa703f..9eefac4 100644
--- a/meta.yaml
+++ b/meta.yaml
@@ -1,10 +1,10 @@
 package:
   name: "open3spn2"
-  version: "0.2.2"
+  version: "0.3.2"
 
 source:
   git_url: https://github.com/cabb99/open3spn2.git
-  git_rev: 0.2.2
+  git_rev: 0.3.2
   git_depth: -1 # (Defaults to -1/not shallow)
 
 build:
diff --git a/open3SPN2/3SPN2.conf b/open3SPN2/3SPN2.conf
index 55bfa76..0682a65 100644
--- a/open3SPN2/3SPN2.conf
+++ b/open3SPN2/3SPN2.conf
@@ -1,258 +1,261 @@
 [Particles]
-name    =   DNA name    epsilon	radius	mass	charge
-row1	=	A	P	0.239006	4.50	94.96960	-0.6
-row2	=	A	S	0.239006	6.20	83.11040	0
-row3	=	A	A	0.239006	4.46	134.1220	0
-row4	=	A	T	0.239006	5.50	125.1078	0
-row5	=	A	G	0.239006	4.20	150.1214	0
-row6	=	A	C	0.239006	5.70	110.0964	0
-row7	=	B	P	0.239006	4.50	94.96960	-0.6
-row8	=	B	S	0.239006	6.20	83.11040	0
-row9	=	B	A	0.239006	5.40	134.1220	0
-row10	=	B	T	0.239006	7.10	125.1078	0
-row11	=	B	G	0.239006	4.90	150.1214	0
-row12	=	B	C	0.239006	6.40	110.0964	0
-row13	=	B_curved	P	0.239006	4.50	94.96960	-0.6
-row14	=	B_curved	S	0.239006	6.20	83.11040	0
-row15	=	B_curved	A	0.239006	5.40	134.1220	0
-row16	=	B_curved	T	0.239006	7.10	125.1078	0
-row17	=	B_curved	G	0.239006	4.90	150.1214	0
-row18	=	B_curved	C	0.239006	6.40	110.0964	0
+name    =    DNA         name epsilon radius mass        charge
+#unit   =    nan         nan  kJ/mol  nm     Da          e
+row1    =    A           P    1      .450    94.96960    -0.6
+row2    =    A           S    1      .620    83.11040    0
+row3    =    A           A    1      .446    134.1220    0
+row4    =    A           T    1      .550    125.1078    0
+row5    =    A           G    1      .420    150.1214    0
+row6    =    A           C    1      .570    110.0964    0
+row7    =    B           P    1      .450    94.96960    -0.6
+row8    =    B           S    1      .620    83.11040    0
+row9    =    B           A    1      .540    134.1220    0
+row10   =    B           T    1      .710    125.1078    0
+row11   =    B           G    1      .490    150.1214    0
+row12   =    B           C    1      .640    110.0964    0
+row13   =    B_curved    P    1      .450    94.96960    -0.6
+row14   =    B_curved    S    1      .620    83.11040    0
+row15   =    B_curved    A    1      .540    134.1220    0
+row16   =    B_curved    T    1      .710    125.1078    0
+row17   =    B_curved    G    1      .490    150.1214    0
+row18   =    B_curved    C    1      .640    110.0964    0
 
 [Bonds]
-name = DNA i j s1 r0 Kb2 Kb3 Kb4
-#unit= nan nan bp nm kJ/nm2/mol kJ/nm3/mol kJ/nm4/mol
+name    =    DNA      i  j  s1 r0     Kb2        Kb3        Kb4
+#unit   =    nan      bp bp bp nm     kJ/nm2/mol kJ/nm3/mol kJ/nm4/mol
 #A DNA
-row1 = A P S 0 0.4157 60 0 600000 #Phosphate - Sugar 3-4,6-7
-row2 = A S P 1 0.3780 60 0 600000 #Sugar - Phosphate 1-3,4-6
-row3 = A S A 0 0.4697 60 0 600000 #Sugar - Adenine
-row4 = A S T 0 0.4220 60 0 600000 #Sugar - Thymine
-row5 = A S G 0 0.4852 60 0 600000 #Sugar - Guanine
-row6 = A S C 0 0.4066 60 0 600000 #Sugar - Cytosine
+row1    =    A        P  S  0  0.4157 60         0          600000     #Phosphate - Sugar 3-4,6-7
+row2    =    A        S  P  1  0.3780 60         0          600000     #Sugar - Phosphate 1-3,4-6
+row3    =    A        S  A  0  0.4697 60         0          600000     #Sugar - Adenine
+row4    =    A        S  T  0  0.4220 60         0          600000     #Sugar - Thymine
+row5    =    A        S  G  0  0.4852 60         0          600000     #Sugar - Guanine
+row6    =    A        S  C  0  0.4066 60         0          600000     #Sugar - Cytosine
 #B DNA
-row7 = B P S 0 0.3899 60 0 600000
-row8 = B S P 1 0.3559 60 0 600000
-row9 = B S A 0 0.4670 60 0 600000
-row10= B S T 0 0.4189 60 0 600000
-row11= B S G 0 0.4829 60 0 600000
-row12= B S C 0 0.4112 60 0 600000
+row7    =    B        P  S  0  0.3899 60         0          600000    
+row8    =    B        S  P  1  0.3559 60         0          600000    
+row9    =    B        S  A  0  0.4670 60         0          600000    
+row10   =    B        S  T  0  0.4189 60         0          600000    
+row11   =    B        S  G  0  0.4829 60         0          600000    
+row12   =    B        S  C  0  0.4112 60         0          600000    
 #B DNA
-row13 = B_curved P S 0 0.3899 60 0 600000
-row14 = B_curved S P 1 0.3559 60 0 600000
-row15 = B_curved S A 0 0.4670 60 0 600000
-row16= B_curved S T 0 0.4189 60 0 600000
-row17= B_curved S G 0 0.4829 60 0 600000
-row18= B_curved S C 0 0.4112 60 0 600000
+row13   =    B_curved P  S  0  0.3899 60         0          600000    
+row14   =    B_curved S  P  1  0.3559 60         0          600000    
+row15   =    B_curved S  A  0  0.4670 60         0          600000    
+row16   =    B_curved S  T  0  0.4189 60         0          600000    
+row17   =    B_curved S  G  0  0.4829 60         0          600000    
+row18   =    B_curved S  C  0  0.4112 60         0          600000    
 
 [Base Stackings]
-name	=	DNA	i	j	k	s1	s2	epsilon	sigma	t0	alpha	rng
-#unit	=	nan	nan	nan	nan	bp	bp	kJ/mol	nm	degrees	1/nm	nan
-#A	DNA
-row1	=	A	S	A	A	0	1	14.39	0.4022	108.32	30	6
-row2	=	A	S	A	T	0	1	14.34	0.3344	96.74	30	6
-row3	=	A	S	A	G	0	1	13.25	0.4261	111.32	30	6
-row4	=	A	S	A	C	0	1	14.51	0.3737	97.36	30	6
-row5	=	A	S	T	A	0	1	10.37	0.4794	103.33	30	6
-row6	=	A	S	T	T	0	1	13.36	0.4031	94.85	30	6
-row7	=	A	S	T	G	0	1	10.34	0.5064	105.36	30	6
-row8	=	A	S	T	C	0	1	12.89	0.4445	94.51	30	6
-row9	=	A	S	G	A	0	1	14.81	0.3855	108.25	30	6
-row10	=	A	S	G	T	0	1	15.57	0.3217	95.59	30	6
-row11	=	A	S	G	G	0	1	14.93	0.4077	111.66	30	6
-row12	=	A	S	G	C	0	1	15.39	0.3592	96.71	30	6
-row13	=	A	S	C	A	0	1	11.42	0.4499	111.39	30	6
-row14	=	A	S	C	T	0	1	12.79	0.3708	102.73	30	6
-row15	=	A	S	C	G	0	1	10.52	0.4772	113.47	30	6
-row16	=	A	S	C	C	0	1	13.24	0.4116	102.14	30	6
-#B	DNA
-row17	=	B	S	A	A	0	1	14.39	0.3716	101.15	30	6
-row18	=	B	S	A	T	0	1	14.34	0.3675	85.94	30	6
-row19	=	B	S	A	G	0	1	13.25	0.3827	105.26	30	6
-row20	=	B	S	A	C	0	1	14.51	0.3744	89.00	30	6
-row21	=	B	S	T	A	0	1	10.37	0.4238	101.59	30	6
-row22	=	B	S	T	T	0	1	13.36	0.3984	89.5	30	6
-row23	=	B	S	T	G	0	1	10.34	0.4416	104.31	30	6
-row24	=	B	S	T	C	0	1	12.89	0.4141	91.28	30	6
-row25	=	B	S	G	A	0	1	14.81	0.3576	100.89	30	6
-row26	=	B	S	G	T	0	1	15.57	0.3598	84.83	30	6
-row27	=	B	S	G	G	0	1	14.93	0.3664	105.48	30	6
-row28	=	B	S	G	C	0	1	15.39	0.3635	88.28	30	6
-row29	=	B	S	C	A	0	1	11.42	0.4038	106.49	30	6
-row30	=	B	S	C	T	0	1	12.79	0.3798	93.31	30	6
-row31	=	B	S	C	G	0	1	10.52	0.4208	109.54	30	6
-row32	=	B	S	C	C	0	1	13.24	0.3935	95.46	30	6
-#B_curved	DNA
-row33	=	B_curved	S	A	A	0	1	13.82	0.358	100.13	30	6
-row34	=	B_curved	S	A	T	0	1	15.05	0.356	90.48	30	6
-row35	=	B_curved	S	A	G	0	1	13.32	0.385	104.39	30	6
-row36	=	B_curved	S	A	C	0	1	15.82	0.345	93.23	30	6
-row37	=	B_curved	S	T	A	0	1	9.15	0.415	102.59	30	6
-row38	=	B_curved	S	T	T	0	1	12.44	0.393	93.32	30	6
-row39	=	B_curved	S	T	G	0	1	9.58	0.432	103.7	30	6
-row40	=	B_curved	S	T	C	0	1	13.11	0.387	94.55	30	6
-row41	=	B_curved	S	G	A	0	1	13.76	0.351	95.45	30	6
-row42	=	B_curved	S	G	T	0	1	14.59	0.347	87.63	30	6
-row43	=	B_curved	S	G	G	0	1	14.77	0.367	106.36	30	6
-row44	=	B_curved	S	G	C	0	1	15.17	0.342	83.12	30	6
-row45	=	B_curved	S	C	A	0	1	9.25	0.415	102.69	30	6
-row46	=	B_curved	S	C	T	0	1	12.42	0.399	96.05	30	6
-row47	=	B_curved	S	C	G	0	1	8.83	0.434	100.46	30	6
-row48	=	B_curved	S	C	C	0	1	14.01	0.384	100.68	30	6
+name    =    DNA      i    j    k    s1   s2   epsilon  sigma    t0        alpha rng
+#unit   =    nan      nan  nan  nan  bp   bp   kJ/mol   nm       degrees   1/nm  nan
+#A    DNA
+row1    =    A        S    A    A    0    1    14.39    0.4022   108.32    30    6
+row2    =    A        S    A    T    0    1    14.34    0.3344    96.74    30    6
+row3    =    A        S    A    G    0    1    13.25    0.4261   111.32    30    6
+row4    =    A        S    A    C    0    1    14.51    0.3737    97.36    30    6
+row5    =    A        S    T    A    0    1    10.37    0.4794   103.33    30    6
+row6    =    A        S    T    T    0    1    13.36    0.4031    94.85    30    6
+row7    =    A        S    T    G    0    1    10.34    0.5064   105.36    30    6
+row8    =    A        S    T    C    0    1    12.89    0.4445    94.51    30    6
+row9    =    A        S    G    A    0    1    14.81    0.3855   108.25    30    6
+row10   =    A        S    G    T    0    1    15.57    0.3217    95.59    30    6
+row11   =    A        S    G    G    0    1    14.93    0.4077   111.66    30    6
+row12   =    A        S    G    C    0    1    15.39    0.3592    96.71    30    6
+row13   =    A        S    C    A    0    1    11.42    0.4499   111.39    30    6
+row14   =    A        S    C    T    0    1    12.79    0.3708   102.73    30    6
+row15   =    A        S    C    G    0    1    10.52    0.4772   113.47    30    6
+row16   =    A        S    C    C    0    1    13.24    0.4116   102.14    30    6
+#B    DNA
+row17   =    B        S    A    A    0    1    14.39    0.3716   101.15    30    6
+row18   =    B        S    A    T    0    1    14.34    0.3675    85.94    30    6
+row19   =    B        S    A    G    0    1    13.25    0.3827   105.26    30    6
+row20   =    B        S    A    C    0    1    14.51    0.3744    89.00    30    6
+row21   =    B        S    T    A    0    1    10.37    0.4238   101.59    30    6
+row22   =    B        S    T    T    0    1    13.36    0.3984    89.5     30    6
+row23   =    B        S    T    G    0    1    10.34    0.4416   104.31    30    6
+row24   =    B        S    T    C    0    1    12.89    0.4141    91.28    30    6
+row25   =    B        S    G    A    0    1    14.81    0.3576   100.89    30    6
+row26   =    B        S    G    T    0    1    15.57    0.3598    84.83    30    6
+row27   =    B        S    G    G    0    1    14.93    0.3664   105.48    30    6
+row28   =    B        S    G    C    0    1    15.39    0.3635    88.28    30    6
+row29   =    B        S    C    A    0    1    11.42    0.4038   106.49    30    6
+row30   =    B        S    C    T    0    1    12.79    0.3798    93.31    30    6
+row31   =    B        S    C    G    0    1    10.52    0.4208   109.54    30    6
+row32   =    B        S    C    C    0    1    13.24    0.3935    95.46    30    6
+#B_curved    DNA
+row33   =    B_curved S    A    A    0    1    13.82    0.358    100.13    30    6
+row34   =    B_curved S    A    T    0    1    15.05    0.356     90.48    30    6
+row35   =    B_curved S    A    G    0    1    13.32    0.385    104.39    30    6
+row36   =    B_curved S    A    C    0    1    15.82    0.345     93.23    30    6
+row37   =    B_curved S    T    A    0    1     9.15    0.415    102.59    30    6
+row38   =    B_curved S    T    T    0    1    12.44    0.393     93.32    30    6
+row39   =    B_curved S    T    G    0    1     9.58    0.432    103.7     30    6
+row40   =    B_curved S    T    C    0    1    13.11    0.387     94.55    30    6
+row41   =    B_curved S    G    A    0    1    13.76    0.351     95.45    30    6
+row42   =    B_curved S    G    T    0    1    14.59    0.347     87.63    30    6
+row43   =    B_curved S    G    G    0    1    14.77    0.367    106.36    30    6
+row44   =    B_curved S    G    C    0    1    15.17    0.342     83.12    30    6
+row45   =    B_curved S    C    A    0    1     9.25    0.415    102.69    30    6
+row46   =    B_curved S    C    T    0    1    12.42    0.399     96.05    30    6
+row47   =    B_curved S    C    G    0    1     8.83    0.434    100.46    30    6
+row48   =    B_curved S    C    C    0    1    14.01    0.384    100.68    30    6
 
 [Harmonic Angles]
-name	=	DNA	i	j	k	s1	s2	epsilon	t0	Base1	Base2	sB
-#unit   =   nan nan nan nan bp  bp  kJ/mol  degrees nan nan     bp
-row0	=	A	S	P	S	1	1	200	92.77	*	*	0
-row1	=	A	P	S	P	0	1	200	91.24	*	*	0
-row2	=	A	A	S	P	0	1	200	104.86	*	*	0
-row3	=	A	T	S	P	0	1	200	110.58	*	*	0
-row4	=	A	G	S	P	0	1	200	103.86	*	*	0
-row5	=	A	C	S	P	0	1	200	106.94	*	*	0
-row6	=	A	P	S	A	0	0	200	103.71	*	*	0
-row7	=	A	P	S	T	0	0	200	93.27	*	*	0
-row8	=	A	P	S	G	0	0	200	107.49	*	*	0
-row9	=	A	P	S	C	0	0	200	97.58	*	*	0
-row10	=	B	S	P	S	1	1	200	94.49	*	*	0
-row11	=	B	P	S	P	0	1	200	120.15	*	*	0
-row12	=	B	A	S	P	0	1	200	112.07	*	*	0
-row13	=	B	T	S	P	0	1	200	116.68	*	*	0
-row14	=	B	G	S	P	0	1	200	110.12	*	*	0
-row15	=	B	C	S	P	0	1	200	114.34	*	*	0
-row16	=	B	P	S	A	0	0	200	103.53	*	*	0
-row17	=	B	P	S	T	0	0	200	92.06	*	*	0
-row18	=	B	P	S	G	0	0	200	107.4	*	*	0
-row19	=	B	P	S	C	0	0	200	96.96	*	*	0
-row20	=	B_curved	S	P	S	1	1	355	nan	A	A	1
-row21	=	B_curved	S	P	S	1	1	464	nan	A	C	1
-row22	=	B_curved	S	P	S	1	1	368	nan	A	G	1
-row23	=	B_curved	S	P	S	1	1	147	nan	A	T	1
-row24	=	B_curved	S	P	S	1	1	273	nan	C	A	1
-row25	=	B_curved	S	P	S	1	1	165	nan	C	C	1
-row26	=	B_curved	S	P	S	1	1	478	nan	C	G	1
-row27	=	B_curved	S	P	S	1	1	368	nan	C	T	1
-row28	=	B_curved	S	P	S	1	1	442	nan	G	A	1
-row29	=	B_curved	S	P	S	1	1	228	nan	G	C	1
-row30	=	B_curved	S	P	S	1	1	165	nan	G	G	1
-row31	=	B_curved	S	P	S	1	1	464	nan	G	T	1
-row32	=	B_curved	S	P	S	1	1	230	nan	T	A	1
-row33	=	B_curved	S	P	S	1	1	442	nan	T	C	1
-row34	=	B_curved	S	P	S	1	1	273	nan	T	G	1
-row35	=	B_curved	S	P	S	1	1	355	nan	T	T	1
-row36	=	B_curved	P	S	P	0	1	300	nan	*	*	0
-row37	=	B_curved	A	S	P	0	1	460	nan	A	A	1
-row38	=	B_curved	A	S	P	0	1	442	nan	A	C	1
-row39	=	B_curved	A	S	P	0	1	358	nan	A	G	1
-row40	=	B_curved	A	S	P	0	1	370	nan	A	T	1
-row41	=	B_curved	T	S	P	0	1	120	nan	T	A	1
-row42	=	B_curved	T	S	P	0	1	383	nan	T	C	1
-row43	=	B_curved	T	S	P	0	1	206	nan	T	G	1
-row44	=	B_curved	T	S	P	0	1	460	nan	T	T	1
-row45	=	B_curved	G	S	P	0	1	383	nan	G	A	1
-row46	=	B_curved	G	S	P	0	1	336	nan	G	C	1
-row47	=	B_curved	G	S	P	0	1	278	nan	G	G	1
-row48	=	B_curved	G	S	P	0	1	442	nan	G	T	1
-row49	=	B_curved	C	S	P	0	1	206	nan	C	A	1
-row50	=	B_curved	C	S	P	0	1	278	nan	C	C	1
-row51	=	B_curved	C	S	P	0	1	278	nan	C	G	1
-row52	=	B_curved	C	S	P	0	1	358	nan	C	T	1
-row53	=	B_curved	P	S	A	0	0	460	nan	A	A	-1
-row54	=	B_curved	P	S	A	0	0	206	nan	A	C	-1
-row55	=	B_curved	P	S	A	0	0	383	nan	A	G	-1
-row56	=	B_curved	P	S	A	0	0	120	nan	A	T	-1
-row57	=	B_curved	P	S	T	0	0	370	nan	T	A	-1
-row58	=	B_curved	P	S	T	0	0	358	nan	T	C	-1
-row59	=	B_curved	P	S	T	0	0	442	nan	T	G	-1
-row60	=	B_curved	P	S	T	0	0	460	nan	T	T	-1
-row61	=	B_curved	P	S	G	0	0	358	nan	G	A	-1
-row62	=	B_curved	P	S	G	0	0	278	nan	G	C	-1
-row63	=	B_curved	P	S	G	0	0	278	nan	G	G	-1
-row64	=	B_curved	P	S	G	0	0	206	nan	G	T	-1
-row65	=	B_curved	P	S	C	0	0	442	nan	C	A	-1
-row66	=	B_curved	P	S	C	0	0	278	nan	C	C	-1
-row67	=	B_curved	P	S	C	0	0	336	nan	C	G	-1
-row68	=	B_curved	P	S	C	0	0	383	nan	C	T	-1
+name    =    DNA      i    j    k    s1   s2   epsilon t0      Base1 Base2 sB
+#unit   =   nan       nan  nan  nan  bp   bp   kJ/mol  degrees nan   nan   bp
+row0    =    A        S    P    S    1    1    200      92.77  *     *      0
+row1    =    A        P    S    P    0    1    200      91.24  *     *      0
+row2    =    A        A    S    P    0    1    200     104.86  *     *      0
+row3    =    A        T    S    P    0    1    200     110.58  *     *      0
+row4    =    A        G    S    P    0    1    200     103.86  *     *      0
+row5    =    A        C    S    P    0    1    200     106.94  *     *      0
+row6    =    A        P    S    A    0    0    200     103.71  *     *      0
+row7    =    A        P    S    T    0    0    200      93.27  *     *      0
+row8    =    A        P    S    G    0    0    200    107.49   *     *      0
+row9    =    A        P    S    C    0    0    200     97.58   *     *      0
+row10   =    B        S    P    S    1    1    200     94.49   *     *      0
+row11   =    B        P    S    P    0    1    200    120.15   *     *      0
+row12   =    B        A    S    P    0    1    200    112.07   *     *      0
+row13   =    B        T    S    P    0    1    200    116.68   *     *      0
+row14   =    B        G    S    P    0    1    200    110.12   *     *      0
+row15   =    B        C    S    P    0    1    200    114.34   *     *      0
+row16   =    B        P    S    A    0    0    200    103.53   *     *      0
+row17   =    B        P    S    T    0    0    200     92.06   *     *      0
+row18   =    B        P    S    G    0    0    200    107.4    *     *      0
+row19   =    B        P    S    C    0    0    200     96.96   *     *      0
+row20   =    B_curved S    P    S    1    1    355    nan      A     A      1
+row21   =    B_curved S    P    S    1    1    464    nan      A     C      1
+row22   =    B_curved S    P    S    1    1    368    nan      A     G      1
+row23   =    B_curved S    P    S    1    1    147    nan      A     T      1
+row24   =    B_curved S    P    S    1    1    273    nan      C     A      1
+row25   =    B_curved S    P    S    1    1    165    nan      C     C      1
+row26   =    B_curved S    P    S    1    1    478    nan      C     G      1
+row27   =    B_curved S    P    S    1    1    368    nan      C     T      1
+row28   =    B_curved S    P    S    1    1    442    nan      G     A      1 
+row29   =    B_curved S    P    S    1    1    228    nan      G     C      1
+row30   =    B_curved S    P    S    1    1    165    nan      G     G      1
+row31   =    B_curved S    P    S    1    1    464    nan      G     T      1
+row32   =    B_curved S    P    S    1    1    230    nan      T     A      1
+row33   =    B_curved S    P    S    1    1    442    nan      T     C      1
+row34   =    B_curved S    P    S    1    1    273    nan      T     G      1 
+row35   =    B_curved S    P    S    1    1    355    nan      T     T      1
+row36   =    B_curved P    S    P    0    1    300    nan      *     *      0
+row37   =    B_curved A    S    P    0    1    460    nan      A     A      1
+row38   =    B_curved A    S    P    0    1    442    nan      A     C      1
+row39   =    B_curved A    S    P    0    1    358    nan      A     G      1
+row40   =    B_curved A    S    P    0    1    370    nan      A     T      1
+row41   =    B_curved T    S    P    0    1    120    nan      T     A      1
+row42   =    B_curved T    S    P    0    1    383    nan      T     C      1
+row43   =    B_curved T    S    P    0    1    206    nan      T     G      1
+row44   =    B_curved T    S    P    0    1    460    nan      T     T      1
+row45   =    B_curved G    S    P    0    1    383    nan      G     A      1
+row46   =    B_curved G    S    P    0    1    336    nan      G     C      1
+row47   =    B_curved G    S    P    0    1    278    nan      G     G      1
+row48   =    B_curved G    S    P    0    1    442    nan      G     T      1
+row49   =    B_curved C    S    P    0    1    206    nan      C     A      1
+row50   =    B_curved C    S    P    0    1    278    nan      C     C      1
+row51   =    B_curved C    S    P    0    1    278    nan      C     G      1
+row52   =    B_curved C    S    P    0    1    358    nan      C     T      1
+row53   =    B_curved P    S    A    0    0    460    nan      A     A     -1
+row54   =    B_curved P    S    A    0    0    206    nan      A     C     -1
+row55   =    B_curved P    S    A    0    0    383    nan      A     G     -1
+row56   =    B_curved P    S    A    0    0    120    nan      A     T     -1
+row57   =    B_curved P    S    T    0    0    370    nan      T     A     -1
+row58   =    B_curved P    S    T    0    0    358    nan      T     C     -1
+row59   =    B_curved P    S    T    0    0    442    nan      T     G     -1
+row60   =    B_curved P    S    T    0    0    460    nan      T     T     -1
+row61   =    B_curved P    S    G    0    0    358    nan      G     A     -1
+row62   =    B_curved P    S    G    0    0    278    nan      G     C     -1
+row63   =    B_curved P    S    G    0    0    278    nan      G     G     -1
+row64   =    B_curved P    S    G    0    0    206    nan      G     T     -1
+row65   =    B_curved P    S    C    0    0    442    nan      C     A     -1
+row66   =    B_curved P    S    C    0    0    278    nan      C     C     -1
+row67   =    B_curved P    S    C    0    0    336    nan      C     G     -1
+row68   =    B_curved P    S    C    0    0    383    nan      C     T     -1 
 
 [Dihedrals]
-name	=	DNA	i	j	k	l	s1	s2	s3	dihedral_type	K_dihedral	K_gaussian	sigma	t0
-#unit	=	nan	nan	nan	nan	nan	bp	bp	bp	nan	kJ/mol	kJ/mol	nan	degrees
-row1	=	A	S	P	S	P	1	1	2	gaussian	0	6	0.3	-9.58
-row2	=	A	P	S	P	S	0	1	1	gaussian	0	6	0.3	-328.4
-row3	=	B	S	P	S	P	1	1	2	gaussian	0	6	0.3	-359.17
-row4	=	B	P	S	P	S	0	1	1	gaussian	0	6	0.3	-334.79
-row5	=	B_curved	S	P	S	P	1	1	2	gaussian_periodic	2	7	0.3
-row6	=	B_curved	P	S	P	S	0	1	1	gaussian_periodic	2	7	0.3
-row7	=	B_curved	A	S	P	S	0	1	1	periodic	2	0	0.3
-row8	=	B_curved	T	S	P	S	0	1	1	periodic	2	0	0.3
-row9	=	B_curved	G	S	P	S	0	1	1	periodic	2	0	0.3
-row10	=	B_curved	C	S	P	S	0	1	1	periodic	2	0	0.3
-row11	=	B_curved	S	P	S	A	1	1	1	periodic	2	0	0.3
-row12	=	B_curved	S	P	S	T	1	1	1	periodic	2	0	0.3
-row13	=	B_curved	S	P	S	G	1	1	1	periodic	2	0	0.3
-row14	=	B_curved	S	P	S	C	1	1	1	periodic	2	0	0.3
+name    =    DNA      i    j    k    l    s1   s2   s3   dihedral_type     K_dihedral K_gaussian sigma t0
+#unit   =    nan      nan  nan  nan  nan  bp   bp   bp   nan               kJ/mol     kJ/mol     nan   degrees
+row1    =    A        S    P    S    P    1    1    2    gaussian          0          6          0.3     -9.58
+row2    =    A        P    S    P    S    0    1    1    gaussian          0          6          0.3   -328.4
+row3    =    B        S    P    S    P    1    1    2    gaussian          0          6          0.3   -359.17
+row4    =    B        P    S    P    S    0    1    1    gaussian          0          6          0.3   -334.79
+row5    =    B_curved S    P    S    P    1    1    2    gaussian_periodic 2          7          0.3
+row6    =    B_curved P    S    P    S    0    1    1    gaussian_periodic 2          7          0.3
+row7    =    B_curved A    S    P    S    0    1    1    periodic          2          0          0.3
+row8    =    B_curved T    S    P    S    0    1    1    periodic          2          0          0.3
+row9    =    B_curved G    S    P    S    0    1    1    periodic          2          0          0.3
+row10   =    B_curved C    S    P    S    0    1    1    periodic          2          0          0.3
+row11   =    B_curved S    P    S    A    1    1    1    periodic          2          0          0.3
+row12   =    B_curved S    P    S    T    1    1    1    periodic          2          0          0.3
+row13   =    B_curved S    P    S    G    1    1    1    periodic          2          0          0.3
+row14   =    B_curved S    P    S    C    1    1    1    periodic          2          0          0.3
 
 
 [Base Pairs]
-name	=	DNA	Base1	Base2	sigma	alpha	rang	epsilon	torsion	t1	t2
-row1	=	A	A	T	5.861	2	12	3.99874	50.17	160.91	140.49
-row2	=	A	G	C	5.528	2	12	5.06241	38.33	165.25	147.11
-row3	=	B	A	T	5.941	2	12	3.99874	-38.35	156.54	135.78
-row4	=	B	G	C	5.53	2	12	5.06241	-42.98	159.81	141.16
-row5	=	B_curved	A	T	5.82	2	12	3.44292	-38.18	153.17	133.51
-row6	=	B_curved	G	C	5.52	2	12	4.35873	-35.75	159.5	138.08
+name    =    DNA      Base1 Base2 sigma  alpha rang epsilon torsion t1     t2
+#unit   =    nan      nan   nan   nm    1/nm  nan  kJ/mol  deg     deg    deg
+row1    =    A        A     T     .5861  20    12   16.731   50.17  160.91 140.49
+row2    =    A        G     C     .5528  20    12   21.181   38.33  165.25 147.11
+row3    =    B        A     T     .5941  20    12   16.731  -38.35  156.54 135.78
+row4    =    B        G     C     .553   20    12   21.181  -42.98  159.81 141.16
+row5    =    B_curved A     T     .582   20    12   14.4052 -38.18  153.17 133.51
+row6    =    B_curved G     C     .552   20    12   18.2368 -35.75  159.5  138.08
 
 [Cross Stackings]
-name = DNA Base_d1 Base_a1 Base_a3 t03 T0CS_1 T0CS_2 Sigma_1 Sigma_2 alpha_cs1 alpha_cs2 rng_cs1 rng_cs2 rng_bp eps_cs1 eps_cs2
-row1  = A	A	T   A	126.57	147.44	130.5	7.344	4.624	4	4	8	8	12	0.522452	0.522452
-row2  = A	A	T   T	126.57	148.97	138.73	8.081	5.095	4	4	8	8	12	0.662942	0.662942
-row3  = A	A	T   G	126.57	146.21	126.68	7.187	4.464	4	4	8	8	12	0.677075	0.712197
-row4  = A	A	T   C	126.57	150.17	134.18	7.990	5.162	4	4	8	8	12	0.466340	0.606830
-row5  = A	T	A   A	126.57	138.42	130.41	8.081	5.095	4	4	8	8	12	0.662942	0.662942
-row6  = A	T	A   T	126.57	141.67	134.68	8.755	5.693	4	4	8	8	12	0.522452	0.522452
-row7  = A	T	A   G	126.57	136.64	127.69	7.952	4.896	4	4	8	8	12	0.606830	0.466340
-row8  = A	T	A   C	126.57	141.64	131.38	8.697	5.724	4	4	8	8	12	0.712197	0.677075
-row9  = A	G	C   A	134.71	147.67	130.57	7.187	4.464	4	4	8	8	12	0.677075	0.712197
-row10 = A	G	C   T	134.71	148.28	140.17	7.952	4.896	4	4	8	8	12	0.606830	0.466340
-row11 = A	G	C   G	134.71	146.84	126.44	7.019	4.315	4	4	8	8	12	0.901943	1.147800
-row12 = A	G	C   C	134.71	150.02	135.31	7.844	4.968	4	4	8	8	12	0.269738	0.269738
-row13 =	A	C	G   A	134.71	145.83	132.69	7.990	5.162	4	4	8	8	12	0.466340	0.606830
-row14 =	A	C	G   T	134.71	148.39	138.21	8.697	5.724	4	4	8	8	12	0.712197	0.677075
-row15 =	A	C	G   G	134.71	144.24	129.73	7.844	4.968	4	4	8	8	12	0.269738	0.269738
-row16 =	A	C	G   C	134.71	148.74	134.45	8.630	5.759	4	4	8	8	12	1.147800	0.901943
-row17 =	B	A	T	A	116.09	154.38	116.88	6.208	5.435	4	4	8	8	12	0.522452	0.522452
-row18 =	B	A	T	T	116.09	159.1	121.74	6.876	6.295	4	4	8	8	12	0.662942	0.662942
-row19 =	B	A	T	G	116.09	152.46	114.23	6.072	5.183	4	4	8	8	12	0.677075	0.712197
-row20 =	B	A	T	C	116.09	158.38	119.06	6.811	6.082	4	4	8	8	12	0.466340	0.60683
-row21 =	B	T	A	A	116.09	147.1	109.42	6.876	6.295	4	4	8	8	12	0.662942	0.662942
-row22 =	B	T	A	T	116.09	153.79	112.95	7.48	7.195	4	4	8	8	12	0.522452	0.522452
-row23 =	B	T	A	G	116.09	144.44	107.32	6.771	6.028	4	4	8	8	12	0.606830	0.46634
-row24 =	B	T	A	C	116.09	151.48	110.56	7.453	6.981	4	4	8	8	12	0.712197	0.677075
-row25 =	B	G	C	A	124.93	154.69	119.34	6.072	5.183	4	4	8	8	12	0.677075	0.712197
-row26 =	B	G	C	T	124.93	157.83	124.72	6.771	6.028	4	4	8	8	12	0.606830	0.46634
-row27 =	B	G	C	G	124.93	153.43	116.51	5.921	4.934	4	4	8	8	12	0.901943	1.1478
-row28 =	B	G	C	C	124.93	158.04	121.98	6.688	5.811	4	4	8	8	12	0.269738	0.269738
-row29 =	B	C	G	A	124.93	152.99	114.6	6.811	6.082	4	4	8	8	12	0.466340	0.60683
-row30 =	B	C	G	T	124.93	159.08	118.26	7.453	6.981	4	4	8	8	12	0.712197	0.677075
-row31 =	B	C	G	G	124.93	150.53	112.45	6.688	5.811	4	4	8	8	12	0.269738	0.269738
-row32 =	B	C	G	C	124.93	157.17	115.88	7.409	6.757	4	4	8	8	12	1.147800	0.901943
-row33	=	B_curved	A	T	A	110.92	154.04	116.34	6.42	5.58	4	4	8	8	12	0.449831	0.449831
-row34	=	B_curved	A	T	T	110.92	158.77	119.61	6.77	6.14	4	4	8	8	12	0.570793	0.570793
-row35	=	B_curved	A	T	G	110.92	153.88	115.19	6.27	5.63	4	4	8	8	12	0.582961	0.613202
-row36	=	B_curved	A	T	C	110.92	157.69	120.92	6.84	6.18	4	4	8	8	12	0.401518	0.52248
-row37	=	B_curved	T	A	A	110.92	148.62	107.40	6.77	6.14	4	4	8	8	12	0.570793	0.570793
-row38	=	B_curved	T	A	T	110.92	155.05	110.76	7.21	6.80	4	4	8	8	12	0.449831	0.449831
-row39	=	B_curved	T	A	G	110.92	147.54	106.33	6.53	6.07	4	4	8	8	12	0.522480	0.401518
-row40	=	B_curved	T	A	C	110.92	153.61	111.57	7.08	6.64	4	4	8	8	12	0.613202	0.582961
-row41	=	B_curved	G	C	A	120.45	153.91	121.61	6.27	5.63	4	4	8	8	12	0.582961	0.613202
-row42	=	B_curved	G	C	T	120.45	155.72	124.92	6.53	6.07	4	4	8	8	12	0.522480	0.401518
-row43	=	B_curved	G	C	G	120.45	151.84	120.52	5.74	5.87	4	4	8	8	12	0.776573	0.988256
-row44	=	B_curved	G	C	C	120.45	157.80	124.88	6.86	5.66	4	4	8	8	12	0.232244	0.232244
-row45	=	B_curved	C	G	A	120.45	152.04	112.45	6.84	6.18	4	4	8	8	12	0.401518	0.52248
-row46	=	B_curved	C	G	T	120.45	157.72	115.43	7.08	6.64	4	4	8	8	12	0.613202	0.582961
-row47	=	B_curved	C	G	G	120.45	151.65	110.51	6.86	5.66	4	4	8	8	12	0.232244	0.232244
-row48	=	B_curved	C	G	C	120.45	154.49	115.80	6.79	6.80	4	4	8	8	12	0.988256	0.776573
+name    =    DNA      Base_d1 Base_a1 Base_a3 t03       T0CS_1    T0CS_2   Sigma_1 Sigma_2 alpha_cs1 alpha_cs2 rng_cs1 rng_cs2 rng_bp eps_cs1   eps_cs2
+#unit   =    nan      nan     nan     nan     deg       deg       deg      kJ      kJ      1/nm      1/nm      nan     nan     nan    kJ/mol    kJ/mol
+row1    =    A        A       T       A       126.57    147.44    130.5    .7344   .4624   40        40        8       8       12     2.185939  2.185939
+row2    =    A        A       T       T       126.57    148.97    138.73   .8081   .5095   40        40        8       8       12     2.773749  2.773749
+row3    =    A        A       T       G       126.57    146.21    126.68   .7187   .4464   40        40        8       8       12     2.832882  2.979832
+row4    =    A        A       T       C       126.57    150.17    134.18   .7990   .5162   40        40        8       8       12     1.951167  2.538977
+row5    =    A        T       A       A       126.57    138.42    130.41   .8081   .5095   40        40        8       8       12     2.773749  2.773749
+row6    =    A        T       A       T       126.57    141.67    134.68   .8755   .5693   40        40        8       8       12     2.185939  2.185939
+row7    =    A        T       A       G       126.57    136.64    127.69   .7952   .4896   40        40        8       8       12     2.538977  1.951167
+row8    =    A        T       A       C       126.57    141.64    131.38   .8697   .5724   40        40        8       8       12     2.979832  2.832882
+row9    =    A        G       C       A       134.71    147.67    130.57   .7187   .4464   40        40        8       8       12     2.832882  2.979832
+row10   =    A        G       C       T       134.71    148.28    140.17   .7952   .4896   40        40        8       8       12     2.538977  1.951167
+row11   =    A        G       C       G       134.71    146.84    126.44   .7019   .4315   40        40        8       8       12     3.773730  4.802395
+row12   =    A        G       C       C       134.71    150.02    135.31   .7844   .4968   40        40        8       8       12     1.128584  1.128584
+row13   =    A        C       G       A       134.71    145.83    132.69   .7990   .5162   40        40        8       8       12     1.951167  2.538977
+row14   =    A        C       G       T       134.71    148.39    138.21   .8697   .5724   40        40        8       8       12     2.979832  2.832882
+row15   =    A        C       G       G       134.71    144.24    129.73   .7844   .4968   40        40        8       8       12     1.128584  1.128584
+row16   =    A        C       G       C       134.71    148.74    134.45   .8630   .5759   40        40        8       8       12     4.802395  3.773730
+row17   =    B        A       T       A       116.09    154.38    116.88   .6208   .5435   40        40        8       8       12     2.185939  2.185939
+row18   =    B        A       T       T       116.09    159.10    121.74   .6876   .6295   40        40        8       8       12     2.773749  2.773749
+row19   =    B        A       T       G       116.09    152.46    114.23   .6072   .5183   40        40        8       8       12     2.832882  2.979832
+row20   =    B        A       T       C       116.09    158.38    119.06   .6811   .6082   40        40        8       8       12     1.951167  2.538977
+row21   =    B        T       A       A       116.09    147.10    109.42   .6876   .6295   40        40        8       8       12     2.773749  2.773749
+row22   =    B        T       A       T       116.09    153.79    112.95   .748    .7195   40        40        8       8       12     2.185939  2.185939
+row23   =    B        T       A       G       116.09    144.44    107.32   .6771   .6028   40        40        8       8       12     2.538977  1.951167
+row24   =    B        T       A       C       116.09    151.48    110.56   .7453   .6981   40        40        8       8       12     2.979832  2.832882
+row25   =    B        G       C       A       124.93    154.69    119.34   .6072   .5183   40        40        8       8       12     2.832882  2.979832
+row26   =    B        G       C       T       124.93    157.83    124.72   .6771   .6028   40        40        8       8       12     2.538977  1.951167
+row27   =    B        G       C       G       124.93    153.43    116.51   .5921   .4934   40        40        8       8       12     3.773730  4.802395
+row28   =    B        G       C       C       124.93    158.04    121.98   .6688   .5811   40        40        8       8       12     1.128584  1.128584
+row29   =    B        C       G       A       124.93    152.99    114.6    .6811   .6082   40        40        8       8       12     1.951167  2.538977
+row30   =    B        C       G       T       124.93    159.08    118.26   .7453   .6981   40        40        8       8       12     2.979832  2.832882
+row31   =    B        C       G       G       124.93    150.53    112.45   .6688   .5811   40        40        8       8       12     1.128584  1.128584
+row32   =    B        C       G       C       124.93    157.17    115.88   .7409   .6757   40        40        8       8       12     4.802395  3.773730
+row33   =    B_curved A       T       A       110.92    154.04    116.34   .642    .558    40        40        8       8       12     1.882093  1.882093
+row34   =    B_curved A       T       T       110.92    158.77    119.61   .677    .614    40        40        8       8       12     2.388198  2.388198
+row35   =    B_curved A       T       G       110.92    153.88    115.19   .627    .563    40        40        8       8       12     2.439109  2.565637
+row36   =    B_curved A       T       C       110.92    157.69    120.92   .684    .618    40        40        8       8       12     1.679951  2.186056
+row37   =    B_curved T       A       A       110.92    148.62    107.40   .677    .614    40        40        8       8       12     2.388198  2.388198
+row38   =    B_curved T       A       T       110.92    155.05    110.76   .721    .680    40        40        8       8       12     1.882093  1.882093
+row39   =    B_curved T       A       G       110.92    147.54    106.33   .653    .607    40        40        8       8       12     2.186056  1.679951
+row40   =    B_curved T       A       C       110.92    153.61    111.57   .708    .664    40        40        8       8       12     2.565637  2.439109
+row41   =    B_curved G       C       A       120.45    153.91    121.61   .627    .563    40        40        8       8       12     2.439109  2.565637
+row42   =    B_curved G       C       T       120.45    155.72    124.92   .653    .607    40        40        8       8       12     2.186056  1.679951
+row43   =    B_curved G       C       G       120.45    151.84    120.52   .574    .587    40        40        8       8       12     3.249181  4.134863
+row44   =    B_curved G       C       C       120.45    157.80    124.88   .686    .566    40        40        8       8       12     0.971709  0.971709
+row45   =    B_curved C       G       A       120.45    152.04    112.45   .684    .618    40        40        8       8       12     1.679951  2.186056
+row46   =    B_curved C       G       T       120.45    157.72    115.43   .708    .664    40        40        8       8       12     2.565637  2.439112
+row47   =    B_curved C       G       G       120.45    151.65    110.51   .686    .566    40        40        8       8       12     0.971709  0.971709
+row48   =    B_curved C       G       C       120.45    154.49    115.80   .679    .680    40        40        8       8       12     4.134863  3.249181
 
 #This section contains the base pair and base step parameters used in 3SPN.2C.
-#The parameters come from the following references:
+#The parameters originate from the following references:
 
 #Morozov, Alexandre V., et al. "Using DNA mechanics to predict in vitro 
 #nucleosome positions and formation energies." Nucleic acids research
@@ -262,89 +265,91 @@ row48	=	B_curved	C	G	C	120.45	154.49	115.80	6.79	6.80	4	4	8	8	12	0.988256	0.7765
 #structures." Computational studies of RNA and DNA. Springer Netherlands, 2006.
 #235-257.
 
+#Updated to correspong to the values from reference:
+#Morozov, Alexandre V., Karissa Fortney, Daria A. Gaykalova, Vasily M. Studitsky, Jonathan Widom,
+#and Eric D. Siggia. 2009. "Using DNA Mechanics to Predict in Vitro Nucleosome Positions and Formation Energies."
+#Nucleic Acids Research 37 (14): 4707-22. https://doi.org/10.1093/nar/gkp475.
+
 [Base Pair Geometry]
-name  = stea steb  buckle  propeller  opening  shear  stretch  stagger
-row1  = A    T     1.8    -15         1.5      0.07  -0.19     0.07
-row2  = T    A    -1.8    -15         1.5     -0.07  -0.19     0.07
-row3  = C    G    -4.9    -8.7       -0.6      0.16  -0.17     0.15
-row4  = G    C     4.9    -8.7       -0.6     -0.16  -0.17     0.15
+name    = stea steb  buckle  propeller  opening  shear  stretch  stagger
+#unit   = nan   nan  deg     deg        deg      A      A        A
+row1    = A    T     1.8    -15         1.5      0.07  -0.19     0.07
+row2    = T    A    -1.8    -15         1.5     -0.07  -0.19     0.07
+row3    = C    G    -4.9    -8.7       -0.6      0.16  -0.17     0.15
+row4    = G    C     4.9    -8.7       -0.6     -0.16  -0.17     0.15
 
 [Base Step Geometry]
-name  = stea  steb  twist   roll   tilt   shift   slide   rise
-row1  = A     A     35.31   0.76  -1.84  -0.05   -0.21    3.27
-row2  = A     T     31.21  -1.39   0.00   0.00   -0.56    3.39
-row3  = A     C     31.52   0.91  -0.64   0.21   -0.54    3.39
-row4  = A     G     33.05   3.15  -1.48   0.12   -0.27    3.38
-row5  = T     A     36.20   5.25   0.00   0.00    0.03    3.34
-#row6  = T     T     35.31   0.91   1.84   0.05   -0.21    3.27 #Previous value
-row6  = T     T     35.31   0.76   1.84   0.05   -0.21    3.27
-row7  = T     C     34.80   3.87   1.52   0.27   -0.03    3.35
-row8  = T     G     35.02   5.95   0.05   0.16    0.18    3.38
-row9  = C     A     35.02   5.95  -0.05  -0.16    0.18    3.38
-row10 = C     T     33.05   3.15   1.48  -0.12   -0.27    3.38
-row11 = C     C     33.17   3.86   0.40   0.02   -0.47    3.28
-#row12 = C     G     34.38   4.29   0.00   0.00    0.57    3.49 #Previous value
-row12 = C     G     35.30   4.29   0.00   0.00    0.57    3.49
-row13 = G     A     34.80   3.87  -1.52  -0.27   -0.03    3.35
-row14 = G     T     31.52   0.91   0.64  -0.21   -0.54    3.39
-row15 = G     C     34.38   0.67   0.00   0.00   -0.07    3.38
-row16 = G     G     33.17   3.86  -0.40  -0.02   -0.47    3.28
-
-#Updated using values from reference:
-#Morozov, Alexandre V., Karissa Fortney, Daria A. Gaykalova, Vasily M. Studitsky, Jonathan Widom,
-#and Eric D. Siggia. 2009. “Using DNA Mechanics to Predict in Vitro Nucleosome Positions and Formation Energies.”
-#Nucleic Acids Research 37 (14): 4707–22. https://doi.org/10.1093/nar/gkp475.
+name    = stea  steb  twist   roll   tilt   shift   slide   rise
+#unit   = nan   nan   deg     deg    deg    A       A       A
+row1    = A     A     35.31   0.76  -1.84  -0.05   -0.21    3.27
+row2    = A     T     31.21  -1.39   0.00   0.00   -0.56    3.39
+row3    = A     C     31.52   0.91  -0.64   0.21   -0.54    3.39
+row4    = A     G     33.05   3.15  -1.48   0.12   -0.27    3.38
+row5    = T     A     36.20   5.25   0.00   0.00    0.03    3.34
+#row6   = T     T     35.31   0.91   1.84   0.05   -0.21    3.27 #Previous value
+row6    = T     T     35.31   0.76   1.84   0.05   -0.21    3.27
+row7    = T     C     34.80   3.87   1.52   0.27   -0.03    3.35
+row8    = T     G     35.02   5.95   0.05   0.16    0.18    3.38
+row9    = C     A     35.02   5.95  -0.05  -0.16    0.18    3.38
+row10   = C     T     33.05   3.15   1.48  -0.12   -0.27    3.38
+row11   = C     C     33.17   3.86   0.40   0.02   -0.47    3.28
+#row12  = C     G     34.38   4.29   0.00   0.00    0.57    3.49 #Previous value
+row12   = C     G     35.30   4.29   0.00   0.00    0.57    3.49
+row13   = G     A     34.80   3.87  -1.52  -0.27   -0.03    3.35
+row14   = G     T     31.52   0.91   0.64  -0.21   -0.54    3.39
+row15   = G     C     34.38   0.67   0.00   0.00   -0.07    3.38
+row16   = G     G     33.17   3.86  -0.40  -0.02   -0.47    3.28
 
 
 # DNA-AWSEM interaction parameters
 [Protein-DNA particles]
-name	=	molecule	DNA	name	epsilon	radius	cutoff	mass	charge
-row1	=	DNA	A	P	0.01991714	6.2	15.5	94.9696	-0.6
-row2	=	DNA	A	S	0.01991714	6.2	15.5	83.1104	0
-row3	=	DNA	A	A	0.01991714	6.2	15.5	134.122	0
-row4	=	DNA	A	T	0.01991714	6.2	15.5	125.1078	0
-row5	=	DNA	A	G	0.01991714	6.2	15.5	150.1214	0
-row6	=	DNA	A	C	0.01991714	6.2	15.5	110.0964	0
-row7	=	DNA	B	P	0.01991714	6.2	15.5	94.9696	-0.6
-row8	=	DNA	B	S	0.01991714	6.2	15.5	83.1104	0
-row9	=	DNA	B	A	0.01991714	6.2	15.5	134.122	0
-row10	=	DNA	B	T	0.01991714	6.2	15.5	125.1078	0
-row11	=	DNA	B	G	0.01991714	6.2	15.5	150.1214	0
-row12	=	DNA	B	C	0.01991714	6.2	15.5	110.0964	0
-row13	=	DNA	B_curved	P	0.01991714	6.2	15.5	94.9696	-0.6
-row14	=	DNA	B_curved	S	0.01991714	6.2	15.5	83.1104	0
-row15	=	DNA	B_curved	A	0.01991714	6.2	15.5	134.122	0
-row16	=	DNA	B_curved	T	0.01991714	6.2	15.5	125.1078	0
-row17	=	DNA	B_curved	G	0.01991714	6.2	15.5	150.1214	0
-row18	=	DNA	B_curved	C	0.01991714	6.2	15.5	110.0964	0
-row19	=	Protein	Protein	CA	0.01991714	6.2	15.5	12	0	#	15
-row20	=	Protein	Protein	N	0	6.2	15.5	0	0	#	16
-row21	=	Protein	Protein	O	0.01991714	6.2	15.5	16	0	#	17
-row22	=	Protein	Protein	CB	0.01991714	6.2	15.5	12	1	#	18	19
-row23	=	Protein	Protein	H	0	6.2	15.5	0	0
-row24	=	Protein	Protein	C	0	6.2	15.5	0	0
+name    = molecule DNA      name epsilon radius cutoff mass     charge
+#unit   = nan      nan      nan  kJ/mol  nm     nm     Da        e
+row1    = DNA      A        P    0.8333  0.62   1.55    94.9696 -0.6
+row2    = DNA      A        S    0.8333  0.62   1.55    83.1104  0
+row3    = DNA      A        A    0.8333  0.62   1.55   134.122   0
+row4    = DNA      A        T    0.8333  0.62   1.55   125.1078  0
+row5    = DNA      A        G    0.8333  0.62   1.55   150.1214  0
+row6    = DNA      A        C    0.8333  0.62   1.55   110.0964  0
+row7    = DNA      B        P    0.8333  0.62   1.55    94.9696 -0.6
+row8    = DNA      B        S    0.8333  0.62   1.55    83.1104  0
+row9    = DNA      B        A    0.8333  0.62   1.55   134.122   0
+row10   = DNA      B        T    0.8333  0.62   1.55   125.1078  0
+row11   = DNA      B        G    0.8333  0.62   1.55   150.1214  0
+row12   = DNA      B        C    0.8333  0.62   1.55   110.0964  0
+row13   = DNA      B_curved P    0.8333  0.62   1.55    94.9696 -0.6
+row14   = DNA      B_curved S    0.8333  0.62   1.55    83.1104  0
+row15   = DNA      B_curved A    0.8333  0.62   1.55   134.122   0
+row16   = DNA      B_curved T    0.8333  0.62   1.55   125.1078  0
+row17   = DNA      B_curved G    0.8333  0.62   1.55   150.1214  0
+row18   = DNA      B_curved C    0.8333  0.62   1.55   110.0964  0
+row19   = Protein  Protein  CA   0.8333  0.62   1.55    12       0    #    15
+row20   = Protein  Protein  N    0       0.62   1.55     0       0    #    16
+row21   = Protein  Protein  O    0.8333  0.62   1.55    16       0    #    17
+row22   = Protein  Protein  CB   0.8333  0.62   1.55    12       1    #    18    19
+row23   = Protein  Protein  H    0       0.62   1.55     0       0
+row24   = Protein  Protein  C    0       0.62   1.55     0       0
 
 [Sequence dependent electrostatics]
-name	=	resname	charge	cutoff
-row1	=	A	0	40
-row2	=	C	0	40
-row3	=	D	-1	40
-row4	=	E	-1	40
-row5	=	F	0	40
-row6	=	G	0	40
-row7	=	H	0	40
-row8	=	I	0	40
-row9	=	K	1	40
-row10	=	L	0	40
-row11	=	M	0	40
-row12	=	N	0	40
-row13	=	P	0	40
-row14	=	Q	0	40
-row15	=	R	1	40
-row16	=	S	0	40
-row17	=	T	0	40
-row18	=	V	0	40
-row19	=	W	0	40
-row20	=	Y	0	40
-
-
+name    = resname charge cutoff
+#unit   = nan     e      nm
+row1    = A        0     4
+row2    = C        0     4
+row3    = D       -1     4
+row4    = E       -1     4
+row5    = F        0     4
+row6    = G        0     4
+row7    = H        0     4
+row8    = I        0     4
+row9    = K        1     4
+row10   = L        0     4
+row11   = M        0     4
+row12   = N        0     4
+row13   = P        0     4
+row14   = Q        0     4
+row15   = R        1     4
+row16   = S        0     4
+row17   = T        0     4
+row18   = V        0     4
+row19   = W        0     4
+row20   = Y        0     4
\ No newline at end of file
diff --git a/open3SPN2/3SPN2.xml b/open3SPN2/3SPN2.xml
index d3aab72..6bfdf3b 100644
--- a/open3SPN2/3SPN2.xml
+++ b/open3SPN2/3SPN2.xml
@@ -2,45 +2,73 @@
   <AtomTypes>
     <!-- DNA particles -->
     <Type name="P" class="Phosphate" element="P" mass="94.9696"/>
-    <Type name="S" class="Sugar" element="H" mass="83.1104"/>
-    <Type name="A base" class="Base" element="N" mass="134.1220"/>
-    <Type name="T base" class="Base" element="S" mass="125.1078"/>
-    <Type name="G base" class="Base" element="C" mass="150.1214"/>
-    <Type name="C base" class="Base" element="O" mass="110.0964"/>
+    <Type name="S DNA" class="Sugar" element="H" mass="83.1104"/>
+    <Type name="S RNA" class="Sugar" element="D" mass="99.1094"/>
+    <Type name="A Base" class="Base" element="N" mass="134.1220"/>
+    <Type name="T Base" class="Base" element="S" mass="125.1078"/>
+    <Type name="U Base" class="Base" element="S" mass="111.0868"/>
+    <Type name="G Base" class="Base" element="C" mass="150.1214"/>
+    <Type name="C Base" class="Base" element="O" mass="110.0964"/>
     <!-- Ion particles -->
     <Type name="Na" class="Ion" element="Mg" mass="0"/>
     <Type name="Mg" class="Ion" element="Fe" mass="0"/>
     <Type name="Ca" class="Ion" element="C" mass="0"/>
   </AtomTypes>
   <Residues>
+    <!-- DNA residues -->
     <Residue name="DA">
-        <AllowPatch name="Initial"/>
-        <Atom name="P"  type="P"/>
-        <Atom name="S"  type="S"/>
-        <Atom name="A"  type="A base"/>
+      <AllowPatch name="Initial"/>
+      <Atom name="P"  type="P"/>
+      <Atom name="S"  type="S DNA"/>
+      <Atom name="A"  type="A Base"/>
     </Residue>
     <Residue name="DC">
-        <AllowPatch name="Initial"/>
-        <Atom name="P"  type="P"/>
-        <Atom name="S"  type="S"/>
-        <Atom name="C"  type="C base"/>
+      <AllowPatch name="Initial"/>
+      <Atom name="P"  type="P"/>
+      <Atom name="S"  type="S DNA"/>
+      <Atom name="C"  type="C Base"/>
     </Residue>
     <Residue name="DT">
-        <AllowPatch name="Initial"/>
-        <Atom name="P"  type="P"/>
-        <Atom name="S"  type="S"/>
-        <Atom name="T"  type="T base"/>
+      <AllowPatch name="Initial"/>
+      <Atom name="P"  type="P"/>
+      <Atom name="S"  type="S DNA"/>
+      <Atom name="T"  type="T Base"/>
     </Residue>
     <Residue name="DG">
-        <AllowPatch name="Initial"/>
-        <Atom name="P"  type="P"/>
-        <Atom name="S"  type="S"/>
-        <Atom name="G"  type="G base"/>
+      <AllowPatch name="Initial"/>
+      <Atom name="P"  type="P"/>
+      <Atom name="S"  type="S DNA"/>
+      <Atom name="G"  type="G Base"/>
+    </Residue>
+    <!-- RNA residues -->
+    <Residue name="A">
+      <AllowPatch name="Initial"/>
+      <Atom name="P"  type="P"/>
+      <Atom name="Sr"  type="S RNA"/>
+      <Atom name="A"  type="A Base"/>
+    </Residue>
+    <Residue name="C">
+      <AllowPatch name="Initial"/>
+      <Atom name="P"  type="P"/>
+      <Atom name="Sr"  type="S RNA"/>
+      <Atom name="C"  type="C Base"/>
+    </Residue>
+    <Residue name="U">
+      <AllowPatch name="Initial"/>
+      <Atom name="P"  type="P"/>
+      <Atom name="Sr"  type="S RNA"/>
+      <Atom name="U"  type="U Base"/>
+    </Residue>
+    <Residue name="G">
+      <AllowPatch name="Initial"/>
+      <Atom name="P"  type="P"/>
+      <Atom name="Sr"  type="S RNA"/>
+      <Atom name="G"  type="G Base"/>
     </Residue>
   </Residues>
   <Patches>
-      <Patch name="Initial">
-          <RemoveAtom name="P"/>
-      </Patch>
+    <Patch name="Initial">
+      <RemoveAtom name="P"/>
+    </Patch>
   </Patches>
 </ForceField>
diff --git a/open3SPN2/DNA_atomic.csv b/open3SPN2/DNA_atomic.csv
new file mode 100644
index 0000000..fd6d51c
--- /dev/null
+++ b/open3SPN2/DNA_atomic.csv
@@ -0,0 +1,132 @@
+name,base,x,y,z
+P,A,-0.299,9.399,-1.529
+O1P,A,-0.377,10.734,-2.162
+O2P,A,0.714,9.245,-0.460
+O5',A,-1.738,8.985,-0.968
+C5',A,-2.674,8.343,-1.855
+C4',A,-3.346,7.182,-1.148
+O4',A,-2.596,5.941,-1.284
+C3',A,-3.530,7.338,0.361
+O3',A,-4.771,6.752,0.737
+C2',A,-2.398,6.514,0.976
+C1',A,-2.479,5.346,0.000
+N9,A,-1.291,4.498,0.000
+C8,A,0.024,4.897,0.000
+N7,A,0.877,3.902,0.000
+C5,A,0.071,2.771,0.000
+C6,A,0.369,1.398,0.000
+N6,A,1.611,0.909,0.000
+N1,A,-0.668,0.532,0.000
+C2,A,-1.912,1.023,0.000
+N3,A,-2.320,2.290,0.000
+C4,A,-1.267,3.124,0.000
+P,C,-0.275,9.443,-1.528
+O1P,C,-0.347,10.779,-2.161
+O2P,C,0.736,9.284,-0.459
+O5',C,-1.717,9.036,-0.967
+C5',C,-2.657,8.399,-1.854
+C4',C,-3.334,7.242,-1.148
+O4',C,-2.590,5.997,-1.284
+C3',C,-3.517,7.399,0.362
+O3',C,-4.762,6.819,0.737
+C2',C,-2.390,6.569,0.977
+C1',C,-2.477,5.402,0.000
+N1,C,-1.285,4.542,0.000
+C2,C,-1.472,3.158,0.000
+O2,C,-2.628,2.709,0.001
+N3,C,-0.391,2.344,0.000
+C4,C,0.837,2.868,0.000
+N4,C,1.875,2.027,0.001
+C5,C,1.056,4.275,0.000
+C6,C,-0.023,5.068,0.000
+P,G,-0.296,9.452,-1.527
+O1P,G,-0.375,10.788,-2.160
+O2P,G,0.716,9.298,-0.458
+O5',G,-1.736,9.038,-0.967
+C5',G,-2.672,8.396,-1.854
+C4',G,-3.344,7.235,-1.148
+O4',G,-2.593,5.994,-1.284
+C3',G,-3.528,7.391,0.362
+O3',G,-4.769,6.805,0.737
+C2',G,-2.396,6.566,0.977
+C1',G,-2.477,5.399,0.000
+N9,G,-1.289,4.551,0.000
+C8,G,0.023,4.962,0.000
+N7,G,0.870,3.969,0.000
+C5,G,0.071,2.833,0.000
+C6,G,0.424,1.460,0.000
+O6,G,1.554,0.955,0.000
+N1,G,-0.700,0.641,0.000
+C2,G,-1.999,1.087,0.000
+N2,G,-2.949,0.139,-0.001
+N3,G,-2.342,2.364,0.001
+C4,G,-1.265,3.177,0.000
+P,I,-0.296,9.452,-1.527
+O1P,I,-0.375,10.788,-2.160
+O2P,I,0.716,9.298,-0.458
+O5',I,-1.736,9.038,-0.967
+C5',I,-2.672,8.396,-1.854
+C4',I,-3.344,7.235,-1.148
+O4',I,-2.593,5.994,-1.284
+C3',I,-3.528,7.391,0.362
+O3',I,-4.769,6.805,0.737
+C2',I,-2.396,6.566,0.977
+C1',I,-2.477,5.399,0.000
+N9,I,-1.289,4.551,0.000
+C8,I,0.023,4.962,0.000
+N7,I,0.870,3.969,0.000
+C5,I,0.071,2.833,0.000
+C6,I,0.424,1.460,0.000
+O6,I,1.554,0.955,0.000
+N1,I,-0.700,0.641,0.000
+C2,I,-1.999,1.087,0.000
+N3,I,-2.342,2.364,0.001
+C4,I,-1.265,3.177,0.000
+C1',PSU,-2.506,5.371,0.000
+N1,PSU,1.087,4.295,0.000
+C2,PSU,1.037,2.915,0.000
+O2,PSU,2.036,2.217,0.000
+N3,PSU,-0.229,2.383,0.000
+C4,PSU,-1.422,3.076,0.000
+O4,PSU,-2.485,2.453,0.000
+C5,PSU,-1.284,4.500,0.000
+C6,PSU,-0.064,5.048,0.000
+P,T,-0.300,9.407,-1.526
+O1P,T,-0.379,10.743,-2.159
+O2P,T,0.712,9.253,-0.458
+O5',T,-1.740,8.993,-0.966
+C5',T,-2.676,8.351,-1.853
+C4',T,-3.348,7.190,-1.147
+O4',T,-2.597,5.949,-1.283
+C3',T,-3.532,7.346,0.362
+O3',T,-4.773,6.760,0.737
+C2',T,-2.400,6.521,0.977
+C1',T,-2.481,5.354,0.000
+N1,T,-1.284,4.500,0.000
+C2,T,-1.462,3.135,0.000
+O2,T,-2.562,2.608,0.000
+N3,T,-0.298,2.407,0.000
+C4,T,0.994,2.897,0.000
+O4,T,1.944,2.119,0.000
+C5,T,1.106,4.338,0.000
+C5M,T,2.466,4.961,0.001
+C6,T,-0.024,5.057,0.000
+P,U,-0.300,9.408,-1.525
+O1P,U,-0.378,10.743,-2.158
+O2P,U,0.712,9.253,-0.456
+O5',U,-1.740,8.993,-0.965
+C5',U,-2.675,8.351,-1.853
+C4',U,-3.347,7.190,-1.147
+O4',U,-2.597,5.949,-1.283
+C3',U,-3.532,7.346,0.362
+O3',U,-4.773,6.759,0.737
+C2',U,-2.400,6.521,0.978
+C1',U,-2.481,5.354,0.000
+N1,U,-1.284,4.500,0.000
+C2,U,-1.462,3.131,0.000
+O2,U,-2.563,2.608,0.000
+N3,U,-0.302,2.397,0.000
+C4,U,0.989,2.884,0.000
+O4,U,1.935,2.094,-0.001
+C5,U,1.089,4.311,0.000
+C6,U,-0.024,5.053,0.000
diff --git a/open3SPN2/__init__.py b/open3SPN2/__init__.py
index 9bb1ff3..dea601c 100644
--- a/open3SPN2/__init__.py
+++ b/open3SPN2/__init__.py
@@ -1,1 +1,4 @@
 from open3SPN2.ff3SPN2 import *
+
+from . import _version
+__version__ = _version.get_versions()['version']
diff --git a/open3SPN2/_version.py b/open3SPN2/_version.py
index 6b288d4..8424c79 100644
--- a/open3SPN2/_version.py
+++ b/open3SPN2/_version.py
@@ -5,8 +5,9 @@
 # directories (produced by setup.py build) will contain a much shorter file
 # that just contains the computed version number.
 
-# This file is released into the public domain. Generated by
-# versioneer-0.18 (https://github.com/warner/python-versioneer)
+# This file is released into the public domain.
+# Generated by versioneer-0.29
+# https://github.com/python-versioneer/python-versioneer
 
 """Git implementation of _version.py."""
 
@@ -15,9 +16,11 @@
 import re
 import subprocess
 import sys
+from typing import Any, Callable, Dict, List, Optional, Tuple
+import functools
 
 
-def get_keywords():
+def get_keywords() -> Dict[str, str]:
     """Get the keywords needed to look up the version information."""
     # these strings will be replaced by git during git-archive.
     # setup.py/versioneer.py will grep for the variable names, so they must
@@ -33,8 +36,15 @@ def get_keywords():
 class VersioneerConfig:
     """Container for Versioneer configuration parameters."""
 
+    VCS: str
+    style: str
+    tag_prefix: str
+    parentdir_prefix: str
+    versionfile_source: str
+    verbose: bool
 
-def get_config():
+
+def get_config() -> VersioneerConfig:
     """Create, populate and return the VersioneerConfig() object."""
     # these strings are filled in when 'setup.py versioneer' creates
     # _version.py
@@ -43,7 +53,7 @@ def get_config():
     cfg.style = "pep440"
     cfg.tag_prefix = ""
     cfg.parentdir_prefix = "None"
-    cfg.versionfile_source = "project/_version.py"
+    cfg.versionfile_source = "open3SPN2/_version.py"
     cfg.verbose = False
     return cfg
 
@@ -52,13 +62,13 @@ class NotThisMethod(Exception):
     """Exception raised if a method is not valid for the current scenario."""
 
 
-LONG_VERSION_PY = {}
-HANDLERS = {}
+LONG_VERSION_PY: Dict[str, str] = {}
+HANDLERS: Dict[str, Dict[str, Callable]] = {}
 
 
-def register_vcs_handler(vcs, method):  # decorator
-    """Decorator to mark a method as the handler for a particular VCS."""
-    def decorate(f):
+def register_vcs_handler(vcs: str, method: str) -> Callable:  # decorator
+    """Create decorator to mark a method as the handler of a VCS."""
+    def decorate(f: Callable) -> Callable:
         """Store f in HANDLERS[vcs][method]."""
         if vcs not in HANDLERS:
             HANDLERS[vcs] = {}
@@ -67,22 +77,35 @@ def decorate(f):
     return decorate
 
 
-def run_command(commands, args, cwd=None, verbose=False, hide_stderr=False,
-                env=None):
+def run_command(
+    commands: List[str],
+    args: List[str],
+    cwd: Optional[str] = None,
+    verbose: bool = False,
+    hide_stderr: bool = False,
+    env: Optional[Dict[str, str]] = None,
+) -> Tuple[Optional[str], Optional[int]]:
     """Call the given command(s)."""
     assert isinstance(commands, list)
-    p = None
-    for c in commands:
+    process = None
+
+    popen_kwargs: Dict[str, Any] = {}
+    if sys.platform == "win32":
+        # This hides the console window if pythonw.exe is used
+        startupinfo = subprocess.STARTUPINFO()
+        startupinfo.dwFlags |= subprocess.STARTF_USESHOWWINDOW
+        popen_kwargs["startupinfo"] = startupinfo
+
+    for command in commands:
         try:
-            dispcmd = str([c] + args)
+            dispcmd = str([command] + args)
             # remember shell=False, so use git.cmd on windows, not just git
-            p = subprocess.Popen([c] + args, cwd=cwd, env=env,
-                                 stdout=subprocess.PIPE,
-                                 stderr=(subprocess.PIPE if hide_stderr
-                                         else None))
+            process = subprocess.Popen([command] + args, cwd=cwd, env=env,
+                                       stdout=subprocess.PIPE,
+                                       stderr=(subprocess.PIPE if hide_stderr
+                                               else None), **popen_kwargs)
             break
-        except EnvironmentError:
-            e = sys.exc_info()[1]
+        except OSError as e:
             if e.errno == errno.ENOENT:
                 continue
             if verbose:
@@ -93,18 +116,20 @@ def run_command(commands, args, cwd=None, verbose=False, hide_stderr=False,
         if verbose:
             print("unable to find command, tried %s" % (commands,))
         return None, None
-    stdout = p.communicate()[0].strip()
-    if sys.version_info[0] >= 3:
-        stdout = stdout.decode()
-    if p.returncode != 0:
+    stdout = process.communicate()[0].strip().decode()
+    if process.returncode != 0:
         if verbose:
             print("unable to run %s (error)" % dispcmd)
             print("stdout was %s" % stdout)
-        return None, p.returncode
-    return stdout, p.returncode
+        return None, process.returncode
+    return stdout, process.returncode
 
 
-def versions_from_parentdir(parentdir_prefix, root, verbose):
+def versions_from_parentdir(
+    parentdir_prefix: str,
+    root: str,
+    verbose: bool,
+) -> Dict[str, Any]:
     """Try to determine the version from the parent directory name.
 
     Source tarballs conventionally unpack into a directory that includes both
@@ -113,15 +138,14 @@ def versions_from_parentdir(parentdir_prefix, root, verbose):
     """
     rootdirs = []
 
-    for i in range(3):
+    for _ in range(3):
         dirname = os.path.basename(root)
         if dirname.startswith(parentdir_prefix):
             return {"version": dirname[len(parentdir_prefix):],
                     "full-revisionid": None,
                     "dirty": False, "error": None, "date": None}
-        else:
-            rootdirs.append(root)
-            root = os.path.dirname(root)  # up a level
+        rootdirs.append(root)
+        root = os.path.dirname(root)  # up a level
 
     if verbose:
         print("Tried directories %s but none started with prefix %s" %
@@ -130,41 +154,48 @@ def versions_from_parentdir(parentdir_prefix, root, verbose):
 
 
 @register_vcs_handler("git", "get_keywords")
-def git_get_keywords(versionfile_abs):
+def git_get_keywords(versionfile_abs: str) -> Dict[str, str]:
     """Extract version information from the given file."""
     # the code embedded in _version.py can just fetch the value of these
     # keywords. When used from setup.py, we don't want to import _version.py,
     # so we do it with a regexp instead. This function is not used from
     # _version.py.
-    keywords = {}
+    keywords: Dict[str, str] = {}
     try:
-        f = open(versionfile_abs, "r")
-        for line in f.readlines():
-            if line.strip().startswith("git_refnames ="):
-                mo = re.search(r'=\s*"(.*)"', line)
-                if mo:
-                    keywords["refnames"] = mo.group(1)
-            if line.strip().startswith("git_full ="):
-                mo = re.search(r'=\s*"(.*)"', line)
-                if mo:
-                    keywords["full"] = mo.group(1)
-            if line.strip().startswith("git_date ="):
-                mo = re.search(r'=\s*"(.*)"', line)
-                if mo:
-                    keywords["date"] = mo.group(1)
-        f.close()
-    except EnvironmentError:
+        with open(versionfile_abs, "r") as fobj:
+            for line in fobj:
+                if line.strip().startswith("git_refnames ="):
+                    mo = re.search(r'=\s*"(.*)"', line)
+                    if mo:
+                        keywords["refnames"] = mo.group(1)
+                if line.strip().startswith("git_full ="):
+                    mo = re.search(r'=\s*"(.*)"', line)
+                    if mo:
+                        keywords["full"] = mo.group(1)
+                if line.strip().startswith("git_date ="):
+                    mo = re.search(r'=\s*"(.*)"', line)
+                    if mo:
+                        keywords["date"] = mo.group(1)
+    except OSError:
         pass
     return keywords
 
 
 @register_vcs_handler("git", "keywords")
-def git_versions_from_keywords(keywords, tag_prefix, verbose):
+def git_versions_from_keywords(
+    keywords: Dict[str, str],
+    tag_prefix: str,
+    verbose: bool,
+) -> Dict[str, Any]:
     """Get version information from git keywords."""
-    if not keywords:
-        raise NotThisMethod("no keywords at all, weird")
+    if "refnames" not in keywords:
+        raise NotThisMethod("Short version file found")
     date = keywords.get("date")
     if date is not None:
+        # Use only the last line.  Previous lines may contain GPG signature
+        # information.
+        date = date.splitlines()[-1]
+
         # git-2.2.0 added "%cI", which expands to an ISO-8601 -compliant
         # datestamp. However we prefer "%ci" (which expands to an "ISO-8601
         # -like" string, which we must then edit to make compliant), because
@@ -177,11 +208,11 @@ def git_versions_from_keywords(keywords, tag_prefix, verbose):
         if verbose:
             print("keywords are unexpanded, not using")
         raise NotThisMethod("unexpanded keywords, not a git-archive tarball")
-    refs = set([r.strip() for r in refnames.strip("()").split(",")])
+    refs = {r.strip() for r in refnames.strip("()").split(",")}
     # starting in git-1.8.3, tags are listed as "tag: foo-1.0" instead of
     # just "foo-1.0". If we see a "tag: " prefix, prefer those.
     TAG = "tag: "
-    tags = set([r[len(TAG):] for r in refs if r.startswith(TAG)])
+    tags = {r[len(TAG):] for r in refs if r.startswith(TAG)}
     if not tags:
         # Either we're using git < 1.8.3, or there really are no tags. We use
         # a heuristic: assume all version tags have a digit. The old git %d
@@ -190,7 +221,7 @@ def git_versions_from_keywords(keywords, tag_prefix, verbose):
         # between branches and tags. By ignoring refnames without digits, we
         # filter out many common branch names like "release" and
         # "stabilization", as well as "HEAD" and "master".
-        tags = set([r for r in refs if re.search(r'\d', r)])
+        tags = {r for r in refs if re.search(r'\d', r)}
         if verbose:
             print("discarding '%s', no digits" % ",".join(refs - tags))
     if verbose:
@@ -199,6 +230,11 @@ def git_versions_from_keywords(keywords, tag_prefix, verbose):
         # sorting will prefer e.g. "2.0" over "2.0rc1"
         if ref.startswith(tag_prefix):
             r = ref[len(tag_prefix):]
+            # Filter out refs that exactly match prefix or that don't start
+            # with a number once the prefix is stripped (mostly a concern
+            # when prefix is '')
+            if not re.match(r'\d', r):
+                continue
             if verbose:
                 print("picking %s" % r)
             return {"version": r,
@@ -214,7 +250,12 @@ def git_versions_from_keywords(keywords, tag_prefix, verbose):
 
 
 @register_vcs_handler("git", "pieces_from_vcs")
-def git_pieces_from_vcs(tag_prefix, root, verbose, run_command=run_command):
+def git_pieces_from_vcs(
+    tag_prefix: str,
+    root: str,
+    verbose: bool,
+    runner: Callable = run_command
+) -> Dict[str, Any]:
     """Get version from 'git describe' in the root of the source tree.
 
     This only gets called if the git-archive 'subst' keywords were *not*
@@ -225,8 +266,15 @@ def git_pieces_from_vcs(tag_prefix, root, verbose, run_command=run_command):
     if sys.platform == "win32":
         GITS = ["git.cmd", "git.exe"]
 
-    out, rc = run_command(GITS, ["rev-parse", "--git-dir"], cwd=root,
-                          hide_stderr=True)
+    # GIT_DIR can interfere with correct operation of Versioneer.
+    # It may be intended to be passed to the Versioneer-versioned project,
+    # but that should not change where we get our version from.
+    env = os.environ.copy()
+    env.pop("GIT_DIR", None)
+    runner = functools.partial(runner, env=env)
+
+    _, rc = runner(GITS, ["rev-parse", "--git-dir"], cwd=root,
+                   hide_stderr=not verbose)
     if rc != 0:
         if verbose:
             print("Directory %s not under git control" % root)
@@ -234,24 +282,57 @@ def git_pieces_from_vcs(tag_prefix, root, verbose, run_command=run_command):
 
     # if there is a tag matching tag_prefix, this yields TAG-NUM-gHEX[-dirty]
     # if there isn't one, this yields HEX[-dirty] (no NUM)
-    describe_out, rc = run_command(GITS, ["describe", "--tags", "--dirty",
-                                          "--always", "--long",
-                                          "--match", "%s*" % tag_prefix],
-                                   cwd=root)
+    describe_out, rc = runner(GITS, [
+        "describe", "--tags", "--dirty", "--always", "--long",
+        "--match", f"{tag_prefix}[[:digit:]]*"
+    ], cwd=root)
     # --long was added in git-1.5.5
     if describe_out is None:
         raise NotThisMethod("'git describe' failed")
     describe_out = describe_out.strip()
-    full_out, rc = run_command(GITS, ["rev-parse", "HEAD"], cwd=root)
+    full_out, rc = runner(GITS, ["rev-parse", "HEAD"], cwd=root)
     if full_out is None:
         raise NotThisMethod("'git rev-parse' failed")
     full_out = full_out.strip()
 
-    pieces = {}
+    pieces: Dict[str, Any] = {}
     pieces["long"] = full_out
     pieces["short"] = full_out[:7]  # maybe improved later
     pieces["error"] = None
 
+    branch_name, rc = runner(GITS, ["rev-parse", "--abbrev-ref", "HEAD"],
+                             cwd=root)
+    # --abbrev-ref was added in git-1.6.3
+    if rc != 0 or branch_name is None:
+        raise NotThisMethod("'git rev-parse --abbrev-ref' returned error")
+    branch_name = branch_name.strip()
+
+    if branch_name == "HEAD":
+        # If we aren't exactly on a branch, pick a branch which represents
+        # the current commit. If all else fails, we are on a branchless
+        # commit.
+        branches, rc = runner(GITS, ["branch", "--contains"], cwd=root)
+        # --contains was added in git-1.5.4
+        if rc != 0 or branches is None:
+            raise NotThisMethod("'git branch --contains' returned error")
+        branches = branches.split("\n")
+
+        # Remove the first line if we're running detached
+        if "(" in branches[0]:
+            branches.pop(0)
+
+        # Strip off the leading "* " from the list of branches.
+        branches = [branch[2:] for branch in branches]
+        if "master" in branches:
+            branch_name = "master"
+        elif not branches:
+            branch_name = None
+        else:
+            # Pick the first branch that is returned. Good or bad.
+            branch_name = branches[0]
+
+    pieces["branch"] = branch_name
+
     # parse describe_out. It will be like TAG-NUM-gHEX[-dirty] or HEX[-dirty]
     # TAG might have hyphens.
     git_describe = describe_out
@@ -268,7 +349,7 @@ def git_pieces_from_vcs(tag_prefix, root, verbose, run_command=run_command):
         # TAG-NUM-gHEX
         mo = re.search(r'^(.+)-(\d+)-g([0-9a-f]+)$', git_describe)
         if not mo:
-            # unparseable. Maybe git-describe is misbehaving?
+            # unparsable. Maybe git-describe is misbehaving?
             pieces["error"] = ("unable to parse git-describe output: '%s'"
                                % describe_out)
             return pieces
@@ -293,26 +374,27 @@ def git_pieces_from_vcs(tag_prefix, root, verbose, run_command=run_command):
     else:
         # HEX: no tags
         pieces["closest-tag"] = None
-        count_out, rc = run_command(GITS, ["rev-list", "HEAD", "--count"],
-                                    cwd=root)
-        pieces["distance"] = int(count_out)  # total number of commits
+        out, rc = runner(GITS, ["rev-list", "HEAD", "--left-right"], cwd=root)
+        pieces["distance"] = len(out.split())  # total number of commits
 
     # commit date: see ISO-8601 comment in git_versions_from_keywords()
-    date = run_command(GITS, ["show", "-s", "--format=%ci", "HEAD"],
-                       cwd=root)[0].strip()
+    date = runner(GITS, ["show", "-s", "--format=%ci", "HEAD"], cwd=root)[0].strip()
+    # Use only the last line.  Previous lines may contain GPG signature
+    # information.
+    date = date.splitlines()[-1]
     pieces["date"] = date.strip().replace(" ", "T", 1).replace(" ", "", 1)
 
     return pieces
 
 
-def plus_or_dot(pieces):
+def plus_or_dot(pieces: Dict[str, Any]) -> str:
     """Return a + if we don't already have one, else return a ."""
     if "+" in pieces.get("closest-tag", ""):
         return "."
     return "+"
 
 
-def render_pep440(pieces):
+def render_pep440(pieces: Dict[str, Any]) -> str:
     """Build up version string, with post-release "local version identifier".
 
     Our goal: TAG[+DISTANCE.gHEX[.dirty]] . Note that if you
@@ -337,23 +419,71 @@ def render_pep440(pieces):
     return rendered
 
 
-def render_pep440_pre(pieces):
-    """TAG[.post.devDISTANCE] -- No -dirty.
+def render_pep440_branch(pieces: Dict[str, Any]) -> str:
+    """TAG[[.dev0]+DISTANCE.gHEX[.dirty]] .
+
+    The ".dev0" means not master branch. Note that .dev0 sorts backwards
+    (a feature branch will appear "older" than the master branch).
 
     Exceptions:
-    1: no tags. 0.post.devDISTANCE
+    1: no tags. 0[.dev0]+untagged.DISTANCE.gHEX[.dirty]
     """
     if pieces["closest-tag"]:
         rendered = pieces["closest-tag"]
+        if pieces["distance"] or pieces["dirty"]:
+            if pieces["branch"] != "master":
+                rendered += ".dev0"
+            rendered += plus_or_dot(pieces)
+            rendered += "%d.g%s" % (pieces["distance"], pieces["short"])
+            if pieces["dirty"]:
+                rendered += ".dirty"
+    else:
+        # exception #1
+        rendered = "0"
+        if pieces["branch"] != "master":
+            rendered += ".dev0"
+        rendered += "+untagged.%d.g%s" % (pieces["distance"],
+                                          pieces["short"])
+        if pieces["dirty"]:
+            rendered += ".dirty"
+    return rendered
+
+
+def pep440_split_post(ver: str) -> Tuple[str, Optional[int]]:
+    """Split pep440 version string at the post-release segment.
+
+    Returns the release segments before the post-release and the
+    post-release version number (or -1 if no post-release segment is present).
+    """
+    vc = str.split(ver, ".post")
+    return vc[0], int(vc[1] or 0) if len(vc) == 2 else None
+
+
+def render_pep440_pre(pieces: Dict[str, Any]) -> str:
+    """TAG[.postN.devDISTANCE] -- No -dirty.
+
+    Exceptions:
+    1: no tags. 0.post0.devDISTANCE
+    """
+    if pieces["closest-tag"]:
         if pieces["distance"]:
-            rendered += ".post.dev%d" % pieces["distance"]
+            # update the post release segment
+            tag_version, post_version = pep440_split_post(pieces["closest-tag"])
+            rendered = tag_version
+            if post_version is not None:
+                rendered += ".post%d.dev%d" % (post_version + 1, pieces["distance"])
+            else:
+                rendered += ".post0.dev%d" % (pieces["distance"])
+        else:
+            # no commits, use the tag as the version
+            rendered = pieces["closest-tag"]
     else:
         # exception #1
-        rendered = "0.post.dev%d" % pieces["distance"]
+        rendered = "0.post0.dev%d" % pieces["distance"]
     return rendered
 
 
-def render_pep440_post(pieces):
+def render_pep440_post(pieces: Dict[str, Any]) -> str:
     """TAG[.postDISTANCE[.dev0]+gHEX] .
 
     The ".dev0" means dirty. Note that .dev0 sorts backwards
@@ -380,12 +510,41 @@ def render_pep440_post(pieces):
     return rendered
 
 
-def render_pep440_old(pieces):
+def render_pep440_post_branch(pieces: Dict[str, Any]) -> str:
+    """TAG[.postDISTANCE[.dev0]+gHEX[.dirty]] .
+
+    The ".dev0" means not master branch.
+
+    Exceptions:
+    1: no tags. 0.postDISTANCE[.dev0]+gHEX[.dirty]
+    """
+    if pieces["closest-tag"]:
+        rendered = pieces["closest-tag"]
+        if pieces["distance"] or pieces["dirty"]:
+            rendered += ".post%d" % pieces["distance"]
+            if pieces["branch"] != "master":
+                rendered += ".dev0"
+            rendered += plus_or_dot(pieces)
+            rendered += "g%s" % pieces["short"]
+            if pieces["dirty"]:
+                rendered += ".dirty"
+    else:
+        # exception #1
+        rendered = "0.post%d" % pieces["distance"]
+        if pieces["branch"] != "master":
+            rendered += ".dev0"
+        rendered += "+g%s" % pieces["short"]
+        if pieces["dirty"]:
+            rendered += ".dirty"
+    return rendered
+
+
+def render_pep440_old(pieces: Dict[str, Any]) -> str:
     """TAG[.postDISTANCE[.dev0]] .
 
     The ".dev0" means dirty.
 
-    Eexceptions:
+    Exceptions:
     1: no tags. 0.postDISTANCE[.dev0]
     """
     if pieces["closest-tag"]:
@@ -402,7 +561,7 @@ def render_pep440_old(pieces):
     return rendered
 
 
-def render_git_describe(pieces):
+def render_git_describe(pieces: Dict[str, Any]) -> str:
     """TAG[-DISTANCE-gHEX][-dirty].
 
     Like 'git describe --tags --dirty --always'.
@@ -422,7 +581,7 @@ def render_git_describe(pieces):
     return rendered
 
 
-def render_git_describe_long(pieces):
+def render_git_describe_long(pieces: Dict[str, Any]) -> str:
     """TAG-DISTANCE-gHEX[-dirty].
 
     Like 'git describe --tags --dirty --always -long'.
@@ -442,7 +601,7 @@ def render_git_describe_long(pieces):
     return rendered
 
 
-def render(pieces, style):
+def render(pieces: Dict[str, Any], style: str) -> Dict[str, Any]:
     """Render the given version pieces into the requested style."""
     if pieces["error"]:
         return {"version": "unknown",
@@ -456,10 +615,14 @@ def render(pieces, style):
 
     if style == "pep440":
         rendered = render_pep440(pieces)
+    elif style == "pep440-branch":
+        rendered = render_pep440_branch(pieces)
     elif style == "pep440-pre":
         rendered = render_pep440_pre(pieces)
     elif style == "pep440-post":
         rendered = render_pep440_post(pieces)
+    elif style == "pep440-post-branch":
+        rendered = render_pep440_post_branch(pieces)
     elif style == "pep440-old":
         rendered = render_pep440_old(pieces)
     elif style == "git-describe":
@@ -474,7 +637,7 @@ def render(pieces, style):
             "date": pieces.get("date")}
 
 
-def get_versions():
+def get_versions() -> Dict[str, Any]:
     """Get version information or return default if unable to do so."""
     # I am in _version.py, which lives at ROOT/VERSIONFILE_SOURCE. If we have
     # __file__, we can work backwards from there to the root. Some
@@ -495,7 +658,7 @@ def get_versions():
         # versionfile_source is the relative path from the top of the source
         # tree (where the .git directory might live) to this file. Invert
         # this to find the root from __file__.
-        for i in cfg.versionfile_source.split('/'):
+        for _ in cfg.versionfile_source.split('/'):
             root = os.path.dirname(root)
     except NameError:
         return {"version": "0+unknown", "full-revisionid": None,
diff --git a/open3SPN2/ff3SPN2.py b/open3SPN2/ff3SPN2.py
index 75f6376..89499ba 100644
--- a/open3SPN2/ff3SPN2.py
+++ b/open3SPN2/ff3SPN2.py
@@ -8,29 +8,31 @@
 
 
 __author__ = 'Carlos Bueno'
-__version__ = '0.2'
+__version__ = '0.3.2'
 
-import simtk.openmm.app
-import simtk.openmm
-import simtk.unit as unit
+import os
 import configparser
-import numpy as np
+import subprocess
 import itertools
+import numpy as np
+
+import openmm.app
+import openmm
+import openmm.unit as unit
 import scipy.spatial.distance as sdist
-import os
 import pdbfixer
 import pandas
-import subprocess
-import nose
+from .geometry import Transform
+from pathlib import Path
 
-__location__ = os.path.realpath(os.path.join(os.getcwd(), os.path.dirname(__file__)))
+__location__ = Path(__file__).parent.resolve()
 _ef = 1 * unit.kilocalorie / unit.kilojoule  # energy scaling factor
 _df = 1 * unit.angstrom / unit.nanometer  # distance scaling factor
 _af = 1 * unit.degree / unit.radian  # angle scaling factor
 _complement = {'A': 'T', 'T': 'A', 'G': 'C', 'C': 'G'}
 _dnaResidues = ['DA', 'DC', 'DT', 'DG']
 _proteinResidues = ['IPR', 'IGL', 'NGP']
-xml = f'{__location__}/3SPN2.xml'
+xml = __location__/'3SPN2.xml'
 
 
 def parseConfigTable(config_section):
@@ -141,6 +143,26 @@ def pdb2table(pdb):
     atom_list.index = atom_list['serial']
     return atom_list
 
+def writePDB(atoms, pdb_file='clean.pdb'):
+    """ Writes a minimal version of the pdb file needed for openmm """
+    # Compute chain field
+    if type(atoms['chainID'].iloc[0]) is not str:
+        chain_ix = 'ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz'
+        atoms['chainID'] = [chain_ix[i - 1] for i in atoms['chainID']]
+
+    #Complete missing fields
+    for field, value in {'occupancy':0.0, 'beta':0.0, 'segment':'', 'charge':'', 'insertion':''}.items():
+        if field not in atoms:
+            atoms[field]=value
+
+    # Write pdb file
+    with open(pdb_file, 'w+') as pdb:
+        for i, atom in atoms.iterrows():
+            pdb_line = f'ATOM  {i + 1:>5} {atom["name"]:^4} {atom.resname:<3} {atom.chainID:1}{atom.resSeq:>4}{atom.insertion:1}   {atom.x:>8.3f}{atom.y:>8.3f}{atom.z:>8.3f}{atom.occupancy:>6.2f}{atom.beta:>6.2f}      {atom.segment:4}{atom.element:2}{atom.charge:2}'
+            assert len(pdb_line) == 80, 'An item in the atom table is longer than expected'
+            pdb.write(pdb_line + '\n')
+    return pdb_file
+
 class DNA(object):
     """ A Coarse Grained DNA object."""
     def __init__(self, periodic=True):
@@ -151,7 +173,7 @@ def __repr__(self):
         return f'3SPN2 DNA object ({len(self.atoms)} atoms)'
         # print the sequence and the identity of the DNA object
 
-    def parseConfigurationFile(self, configuration_file=f'{__location__}/3SPN2.conf'):
+    def parseConfigurationFile(self, configuration_file=__location__ /'3SPN2.conf'):
         """Reads the configuration file for the forcefield. The default configuration file is 3SPN2.conf
         and it contains most of the parameters used in the simulation."""
         self.configuration_file = configuration_file
@@ -206,40 +228,101 @@ def computeGeometry(self, sequence=None, temp_name='temp'):
             data += [pandas.concat([pandas.Series([f'{s}-{_complement[s]}'], index=['Sequence']), pair_s, step_s])]
             _s = s
         data = pandas.concat(data, axis=1).T
-        try:
-            location_x3dna = os.environ["X3DNA"]
-        except KeyError as ex:
-            raise X3DNAnotFound from ex
-
-        with open(f'{temp_name}_parameters.par', 'w+') as par:
-            par.write(f' {len(data)} # Number of base pairs\n')
-            par.write(f' 0 # local base-pair & step parameters\n')
-            par.write('#')
-            par.write(data.to_csv(sep=' ', index=False))
-
-        # Attempting to call rebuild multiple times
-        # This function fails sometimes when called by multiple because a file with the same name is created
-        attempt = 0
-        max_attempts = 10
-        while attempt < max_attempts:
-            try:
-                subprocess.check_output([f'{location_x3dna}/bin/x3dna_utils',
-                                         'cp_std', 'BDNA'])
-                subprocess.check_output([f'{location_x3dna}/bin/rebuild',
-                                         '-atomic', f'{temp_name}_parameters.par',
-                                         f'{temp_name}_template.pdb'])
-                break
-            except subprocess.CalledProcessError as e:
-                attempt += 1
-                if attempt == max_attempts:
-                    print(f"subprocess.CalledProcessError failed {max_attempts} times {e.args[0]}: {e.args[1]}")
-        template_dna = self.fromPDB(f'{temp_name}_template.pdb', output_pdb=f'{temp_name}_temp.pdb', compute_topology=False)
-        template = template_dna.atoms.copy()
+        # try:
+        #     location_x3dna = os.environ["X3DNA"]
+        # except KeyError as ex:
+        #     raise X3DNAnotFound from ex
+
+        # with open(f'{temp_name}_parameters.par', 'w+') as par:
+        #     par.write(f' {len(data)} # Number of base pairs\n')
+        #     par.write(f' 0 # local base-pair & step parameters\n')
+        #     par.write('#')
+        #     par.write(data.to_csv(sep=' ', index=False))
+
+        # # Attempting to call rebuild multiple times
+        # # This function fails sometimes when called by multiple because a file with the same name is created
+        # attempt = 0
+        # max_attempts = 10
+        # while attempt < max_attempts:
+        #     try:
+        #         subprocess.check_output([f'{location_x3dna}/bin/x3dna_utils',
+        #                                  'cp_std', 'BDNA'])
+        #         subprocess.check_output([f'{location_x3dna}/bin/rebuild',
+        #                                  '-atomic', f'{temp_name}_parameters.par',
+        #                                  f'{temp_name}_template.pdb'])
+        #         break
+        #     except subprocess.CalledProcessError as e:
+        #         attempt += 1
+        #         if attempt == max_attempts:
+        #             print(f"subprocess.CalledProcessError failed {max_attempts} times {e.args[0]}: {e.args[1]}")
+        # template_dna = self.fromPDB(f'{temp_name}_template.pdb', output_pdb=f'{temp_name}_temp.pdb', compute_topology=False)
+        # template = template_dna.atoms.copy()
+
+        new_orien = np.eye(3)
+        new_pos = np.zeros(3)
+
+        forward_strand=[]
+        reverse_strand=[]
+        xyz = pandas.read_csv(__location__/'DNA_atomic.csv')
+
+        for row in data.itertuples():
+            base1 = row.Sequence[0]
+            base2 = row.Sequence[-1]
+            base_parameters = np.array(row[2:8])
+            step_parameters = np.array(row[8:16])
+            
+            base_parameters[3:]*=np.pi/180
+            step_parameters[3:]*=np.pi/180
+            
+            t=Transform(*base_parameters)
+
+            b1 = xyz[xyz['base']==base1].copy()
+            b2 = xyz[xyz['base']==base2].copy()
+            
+            xyz1=b1[['x','y','z']].copy()
+            xyz1 = np.dot(b1[['x','y','z']],t.full_rotation.T)+t.full_displacement
+            xyz1 = np.dot((xyz1-t.half_displacement),t.half_rotation)
+            
+            xyz2 = b2[['x','y','z']].copy()
+            xyz2[['y','z']]=-xyz2[['y','z']]
+            xyz2 = np.dot((xyz2-t.half_displacement),t.half_rotation)
+            
+            # Step parameters
+            t=Transform(*step_parameters)
+            new_pos+=np.dot(t.full_displacement,new_orien)
+            new_orien=np.dot(t.full_rotation.T,new_orien)
+            
+            xyz1 = np.dot(xyz1,new_orien)+new_pos
+            #xyz2 = np.dot(xyz2,new_orien)+new_pos
+            
+            b1[['x','y','z']] = xyz1
+            b2[['x','y','z']] = xyz2
+            
+            for b in [b1,b2]:
+                b['recname']='ATOM'
+                b['resname']='D'+b['base']
+                b['resSeq']=row[0]+1
+                b['occupancy']=1.00
+                b['tempFactor']=1.00
+                b['element']=b['name'].str[0]
+            
+            b1['chainID']='A'
+            b2['chainID']='B'
+
+            forward_strand=forward_strand + [b1]
+            reverse_strand=[b2] + reverse_strand
+            
+
+        template=pandas.concat(forward_strand+reverse_strand)
+        template.reindex()
+        template['serial']=template.index+1
+        writePDB(template,pdb_file='temp_template.pdb')
+
         try:
             self.atoms
         except AttributeError:
             return template
-
+        
         template = template[template['chainID'] == 'A']
         original = self.atoms.copy()
         original.index = original['chainID'].astype(str) + '_' + original['resSeq'].astype(str) + '_' + original['name']
@@ -433,24 +516,12 @@ def computeTopology(self, template_from_X3DNA=True, temp_name='temp'):
             self.dihedrals['t0'] = -d - 180
 
     def writePDB(self, pdb_file='clean.pdb'):
-        """ Writes a minimal version of the pdb file needed for openmm """
-        # Compute chain field
-        if type(self.atoms['chainID'].iloc[0]) is not str:
-            chain_ix = 'ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz'
-            self.atoms['chainID'] = [chain_ix[i - 1] for i in self.atoms['chainID']]
-
         # Compute element fields
         element_ix = {'P': 'P', 'S': 'H', 'A': 'N', 'T': 'S', 'C': 'O', 'G': 'C'}  # Elements choosen to keep VMD colors
         self.atoms.loc[:, 'element'] = [element_ix[atomType] for atomType in self.atoms['name']]
-
-        # Write pdb file
-        with open(pdb_file, 'w+') as pdb:
-            for i, atom in self.atoms.iterrows():
-                pdb_line = f'ATOM  {i + 1:>5} {atom["name"]:^4} {atom.resname:<3} {atom.chainID}{atom.resSeq:>4}    {atom.x:>8.3f}{atom.y:>8.3f}{atom.z:>8.3f}' + ' ' * 22 + f'{atom.element:2}' + ' ' * 2
-                assert len(pdb_line) == 80, 'An item in the atom table is longer than expected'
-                pdb.write(pdb_line + '\n')
-        self.pdb_file = pdb_file
-        return pdb_file
+        
+        self.pdb_file = writePDB(self.atoms,pdb_file)
+        return self.pdb_file
 
     @classmethod
     def fromCoarsePDB(cls, pdb_file, dna_type='B_curved', template_from_X3DNA=True, temp_name='temp',
@@ -486,7 +557,7 @@ def fromPDB(cls, pdb_file, dna_type='B_curved', template_from_X3DNA=True, output
         return self
 
     @staticmethod
-    def CoarseGrain(pdb_table):
+    def CoarseGrain(pdb_table, rna=False):
         """ Selects DNA atoms from a pdb table and returns a table containing only the coarse-grained atoms for 3SPN2"""
         masses = {"H": 1.00794, "C": 12.0107, "N": 14.0067, "O": 15.9994, "P": 30.973762, }
         CG = {"O5\'": 'P', "C5\'": 'S', "C4\'": 'S', "O4\'": 'S', "C3\'": 'S', "O3\'": 'P',
@@ -499,6 +570,8 @@ def CoarseGrain(pdb_table):
               "H1'": 'S', "H8": 'B', "H61": 'B', "H62": 'B', 'H2': 'B', 'H1': 'B', 'H21': 'B',
               'H22': 'B', 'H3': 'B', 'H71': 'B', 'H72': 'B', 'H73': 'B', 'H6': 'B', 'H41': 'B',
               'H42': 'B', 'H5': 'B', "HO3'": 'P'}
+        if rna:
+            CG.update({"HO2\'": 'S', "O2\'": 'S'})
         cols = ['recname', 'serial', 'name', 'altLoc',
                 'resname', 'chainID', 'resSeq', 'iCode',
                 'x', 'y', 'z', 'occupancy', 'tempFactor',
@@ -506,7 +579,10 @@ def CoarseGrain(pdb_table):
         temp = pdb_table.copy()
 
         # Select DNA residues
-        temp = temp[temp['resname'].isin(['DA', 'DT', 'DG', 'DC'])]
+        if rna:
+            temp = temp[temp['resname'].isin(['A', 'U', 'G', 'C'])]
+        else:
+            temp = temp[temp['resname'].isin(['DA', 'DT', 'DG', 'DC'])]
 
         # Group the atoms by sugar, phosphate or base
         temp['group'] = temp['name'].replace(CG)
@@ -522,7 +598,7 @@ def CoarseGrain(pdb_table):
         temp = temp[temp['resname'] != 'remove']
 
         # Calculate center of mass
-        temp['element']=temp['element'].str.strip()
+        temp['element'] = temp['element'].str.strip()
         temp['mass'] = temp.element.replace(masses).astype(float)
         temp[['x', 'y', 'z']] = (temp[['x', 'y', 'z']].T * temp['mass']).T[['x', 'y', 'z']]
         temp = temp[temp['element'] != 'H']  # Exclude hydrogens
@@ -540,13 +616,17 @@ def CoarseGrain(pdb_table):
         Coarse.loc[mask, 'name'] = Coarse[mask].resname.str[-1]  # takes last letter from the residue name
         Coarse['type'] = Coarse['name']
         # Set element (depends on base)
-        Coarse['element'] = Coarse['name'].replace({'P': 'P', 'S': 'H', 'A': 'N', 'T': 'S', 'G': 'C', 'C': 'O'})
-        # Remove P from the beggining
+        if rna:
+            Coarse['name'] = Coarse['name'].replace({'S', 'Sr'})
+        Coarse['element'] = Coarse['name'].replace({'P': 'P', 'S': 'H', 'Sr': 'D',
+                                                    'A': 'N', 'T': 'S', 'U': 'S', 'G': 'C', 'C': 'O'})
+        # Remove P from the beginning
         drop_list = []
         for chain in Coarse.chainID.unique():
             sel = Coarse[Coarse.chainID == chain]
             drop_list += list(sel[(sel.resSeq == sel.resSeq.min()) & sel['name'].isin(['P'])].index)
         Coarse = Coarse.drop(drop_list)
+
         # Renumber
         Coarse.index = range(len(Coarse))
         Coarse['serial'] = Coarse.index
@@ -627,19 +707,19 @@ def fromXYZ(cls, xyz_file, dnatype='B_curved', template_from_X3DNA=True, output_
         return self
 
 
-class System(simtk.openmm.System):
+class System(openmm.System):
     """ Wrapper of openmm system class, adds some openmm simulation attributes"""
 
-    def __init__(self, dna, forcefieldFiles=[f'{__location__}/3SPN2.xml'], periodicBox=None):
+    def __init__(self, dna, forcefieldFiles=[__location__/'3SPN2.xml'], periodicBox=None):
         self.dna = dna
-        self.top = simtk.openmm.app.PDBFile(dna.pdb_file).getTopology()
+        self.top = openmm.app.PDBFile(dna.pdb_file).getTopology()
         if self.top.getUnitCellDimensions() is None:
             x = dna.atoms[['x', 'y', 'z']]
             d = np.round((x.max() - x.min()) * 2 + 5, -1)
             self.top.setUnitCellDimensions(d)
 
-        self.coord = simtk.openmm.app.PDBFile(dna.pdb_file)
-        self.forcefield = simtk.openmm.app.ForceField(*forcefieldFiles)
+        self.coord = openmm.app.PDBFile(dna.pdb_file)
+        self.forcefield = openmm.app.ForceField(*forcefieldFiles)
         self._wrapped_system = self.forcefield.createSystem(self.top)
         self.periodicBox = periodicBox
         if periodicBox is not None:
@@ -660,7 +740,7 @@ def clearForces(self, keepCMMotionRemover=True):
         openMM usually adds a "CMMotionRemover" force to keep the center of mass of the system from drifting."""
         j = 0
         for i, f in enumerate(self.getForces()):
-            if keepCMMotionRemover and i == 0 and f.__class__ == simtk.openmm.CMMotionRemover:
+            if keepCMMotionRemover and i == 0 and f.__class__ == openmm.CMMotionRemover:
                 # print('Kept ', f.__class__)
                 j += 1
                 continue
@@ -695,9 +775,9 @@ def addProteinDNAforces(self, verbose=False):
 
     def initializeMD(self, temperature=300 * unit.kelvin, platform_name='Reference', damping=2/unit.picosecond, timestep=2*unit.femtoseconds):
         """Starts a simple simulation using the selected system"""
-        self.integrator = simtk.openmm.LangevinIntegrator(temperature, damping, timestep)
-        self.platform = simtk.openmm.Platform.getPlatformByName(platform_name)
-        self.simulation = simtk.openmm.app.Simulation(self.top, self._wrapped_system, self.integrator, self.platform)
+        self.integrator = openmm.LangevinIntegrator(temperature, damping, timestep)
+        self.platform = openmm.Platform.getPlatformByName(platform_name)
+        self.simulation = openmm.app.Simulation(self.top, self._wrapped_system, self.integrator, self.platform)
         self.simulation.context.setPositions(self.coord.positions)
         return self.simulation
 
@@ -795,7 +875,7 @@ def computeSingleEnergy(self, system, trajectory):
         # return a table with the energy
 
 
-class Bond(Force, simtk.openmm.CustomBondForce):
+class Bond(Force, openmm.CustomBondForce):
     def __init__(self, dna, force_group=6, OpenCLPatch=True):
         self.force_group = force_group
         super().__init__(dna, OpenCLPatch=OpenCLPatch)
@@ -810,7 +890,7 @@ def getParameterNames(self):
         return [self.perInteractionParameters, self.GlobalParameters]
 
     def reset(self):
-        bondForce = simtk.openmm.CustomBondForce("Kb2*(r-r0)^2+Kb3*(r-r0)^3+Kb4*(r-r0)^4")
+        bondForce = openmm.CustomBondForce("Kb2*(r-r0)^2+Kb3*(r-r0)^3+Kb4*(r-r0)^4")
         bondForce.addPerBondParameter('r0')
         bondForce.addPerBondParameter('Kb2')
         bondForce.addPerBondParameter('Kb3')
@@ -829,14 +909,14 @@ def defineInteraction(self):
             self.force.addBond(int(b['aai']), int(b['aaj']), parameters)
 
 
-class Angle(Force, simtk.openmm.HarmonicAngleForce):
+class Angle(Force, openmm.HarmonicAngleForce):
 
     def __init__(self, dna, force_group=7, OpenCLPatch=True):
         self.force_group = force_group
         super().__init__(dna, OpenCLPatch=OpenCLPatch)
 
     def reset(self):
-        angleForce = simtk.openmm.HarmonicAngleForce()
+        angleForce = openmm.HarmonicAngleForce()
         angleForce.setUsesPeriodicBoundaryConditions(self.periodic)
         angleForce.setForceGroup(self.force_group)
         self.force = angleForce
@@ -848,13 +928,13 @@ def defineInteraction(self):
             self.force.addAngle(int(a['aai']), int(a['aaj']), int(a['aak']), *parameters)
 
 
-class Stacking(Force, simtk.openmm.CustomCompoundBondForce):
+class Stacking(Force, openmm.CustomCompoundBondForce):
     def __init__(self, dna, force_group=8, OpenCLPatch=True):
         self.force_group = force_group
         super().__init__(dna, OpenCLPatch=OpenCLPatch)
 
     def reset(self):
-        stackingForce = simtk.openmm.CustomCompoundBondForce(3, """energy;
+        stackingForce = openmm.CustomCompoundBondForce(3, """energy;
                         energy=rep+f2*attr;
                         rep=epsilon*(1-exp(-alpha*(dr)))^2*step(-dr);
                         attr=epsilon*(1-exp(-alpha*(dr)))^2*step(dr)-epsilon;
@@ -884,13 +964,13 @@ def defineInteraction(self):
             self.force.addBond([a['aai'], a['aaj'], a['aak']], parameters)
 
 
-class Dihedral(Force, simtk.openmm.CustomTorsionForce):
+class Dihedral(Force, openmm.CustomTorsionForce):
     def __init__(self, dna, force_group=9, OpenCLPatch=True):
         self.force_group = force_group
         super().__init__(dna, OpenCLPatch=OpenCLPatch)
 
     def reset(self):
-        dihedralForce = simtk.openmm.CustomTorsionForce("""energy;
+        dihedralForce = openmm.CustomTorsionForce("""energy;
                         energy = K_periodic*(1-cs)-K_gaussian*exp(-dt_periodic^2/2/sigma^2);
                         cs = cos(dt);
                         dt_periodic = dt-floor((dt+pi)/(2*pi))*(2*pi);
@@ -915,14 +995,14 @@ def defineInteraction(self):
             self.force.addTorsion(*particles, parameters)
 
 
-class BasePair(Force, simtk.openmm.CustomHbondForce):
+class BasePair(Force, openmm.CustomHbondForce):
     def __init__(self, dna, force_group=10, OpenCLPatch=True):
         self.force_group = force_group
         super().__init__(dna, OpenCLPatch=OpenCLPatch)
 
     def reset(self):
         def basePairForce():
-            pairForce = simtk.openmm.CustomHbondForce('''energy;
+            pairForce = openmm.CustomHbondForce('''energy;
                         energy=rep+1/2*(1+cos(dphi))*fdt1*fdt2*attr;
                         rep  = epsilon*(1-exp(-alpha*dr))^2*(1-step(dr));
                         attr = epsilon*(1-exp(-alpha*dr))^2*step(dr)-epsilon;
@@ -995,12 +1075,12 @@ def defineInteraction(self):
 
             # Define parameters
             parameters = [pair.torsion * _af,
-                          pair.sigma * _df,
+                          pair.sigma,
                           pair.t1 * _af,
                           pair.t2 * _af,
                           pair.rang,
-                          pair.epsilon * _ef,
-                          pair.alpha / _df]
+                          pair.epsilon,
+                          pair.alpha]
 
             # Add donors and acceptors
             # Here I am including the same atom twice,
@@ -1037,7 +1117,7 @@ def __init__(self, dna, force_group=11, OpenCLPatch=True):
 
     def reset(self):
         def crossStackingForce(parametersOnDonor=False):
-            crossForce = simtk.openmm.CustomHbondForce(f'''energy;
+            crossForce = openmm.CustomHbondForce(f'''energy;
                          energy   = fdt3*fdtCS*attr/2;
                          attr     = epsilon*(1-exp(-alpha*dr))^2*step(dr)-epsilon;
                          fdt3     = max(f1*pair0t3,pair1t3);
@@ -1132,9 +1212,9 @@ def defineInteraction(self):
                           param['T0CS_2'] * _af,
                           param['rng_cs2'],
                           param['rng_bp'],
-                          param['eps_cs2'] * _ef,
-                          param['alpha_cs2'] / _df,
-                          param['Sigma_2'] * _df]
+                          param['eps_cs2'],
+                          param['alpha_cs2'],
+                          param['Sigma_2']]
             c1.addDonor(d1, d2, d3)
             c2.addAcceptor(d1, d2, d3, parameters)
             # print("Donor", d1t, d1, d2, d3)
@@ -1152,9 +1232,9 @@ def defineInteraction(self):
                           param['T0CS_1'] * _af,
                           param['rng_cs1'],
                           param['rng_bp'],
-                          param['eps_cs1'] * _ef,
-                          param['alpha_cs1'] / _df,
-                          param['Sigma_1'] * _df]
+                          param['eps_cs1'],
+                          param['alpha_cs1'],
+                          param['Sigma_1']]
             c1.addAcceptor(a1, a2, a3, parameters)
             c2.addDonor(a1, a2, a3)
             # print("Aceptor", a1t, a1, a2, a3)
@@ -1192,33 +1272,49 @@ def getForceGroup(self):
             assert fg == c2.getForceGroup()
         return fg
 
-
 def addNonBondedExclusions(dna, force, OpenCLPatch=True):
     is_dna = dna.atoms['resname'].isin(_dnaResidues)
     atoms = dna.atoms.copy()
-    selection = atoms[is_dna]
-    for (i, atom_a), (j, atom_b) in itertools.combinations(selection.iterrows(), r=2):
-        if j < i:
-            i, j = j, i
-            atom_a, atom_b = atom_b, atom_a
-        # Neighboring residues
-        if atom_a.chainID == atom_b.chainID and (abs(atom_a.resSeq - atom_b.resSeq) <= 1):
-            force.addExclusion(i, j)
-            # print(i, j)
-        # Base-pair residues
-        elif OpenCLPatch and (atom_a['name'] in _complement.keys()) and (atom_b['name'] in _complement.keys()) and (
-                atom_a['name'] == _complement[atom_b['name']]):
-            force.addExclusion(i, j)
-            # print(i, j)
-
-
-class Exclusion(Force, simtk.openmm.CustomNonbondedForce):
+    selection = atoms[is_dna].sort_index()
+    selection['index'] = selection.index
+    selection['neighbor'] = selection['chainID'].astype(str) + '_' + (selection['resSeq'] - 1).astype(str)
+    selection.index = selection['chainID'].astype(str) + '_' + (selection['resSeq']).astype(str)
+
+    exclusions = []
+    for i, neighbor_res, self_res in zip(selection['index'], selection['neighbor'], selection.index):
+        # Add exclusions for the same residue
+        for j in selection.loc[self_res, 'index']:
+            if i > j:
+                exclusions += [(j, i)]
+        # Add exclusions with the neighboring residue on the same chain
+        try:
+            for j in selection.loc[neighbor_res, 'index']:
+                exclusions += [(j, i)]
+        except KeyError:
+            continue
+
+    exclusions_bp = []
+    if OpenCLPatch:
+        # Add basepair exclusions
+        for key in _complement:
+            selection_N = atoms[is_dna & (atoms['name'] == key)]
+            selection_C = atoms[is_dna & (atoms['name'] == _complement[key])]
+            for i, j in itertools.product(selection_N.index, selection_C.index):
+                if i > j:
+                    exclusions_bp += [(j, i)]
+
+    exclusions = list(set(exclusions+exclusions_bp))
+    for i, j in exclusions:
+        force.addExclusion(i, j)
+
+
+class Exclusion(Force, openmm.CustomNonbondedForce):
     def __init__(self, dna, force_group = 12, OpenCLPatch=True):
         self.force_group = force_group
         super().__init__(dna, OpenCLPatch=OpenCLPatch)
 
     def reset(self):
-        exclusionForce = simtk.openmm.CustomNonbondedForce("""energy;
+        exclusionForce = openmm.CustomNonbondedForce("""energy;
                          energy=(epsilon*((sigma/r)^12-2*(sigma/r)^6)+epsilon)*step(sigma-r);
                          sigma=0.5*(sigma1+sigma2);
                          epsilon=sqrt(epsilon1*epsilon2)""")
@@ -1238,7 +1334,7 @@ def defineInteraction(self):
         particle_definition.index = particle_definition.name
 
         # Reduces or increases the cutoff to the maximum particle radius
-        self.force.setCutoffDistance(particle_definition.radius.max() * _df)
+        self.force.setCutoffDistance(particle_definition.radius.max())
 
         # Select only dna atoms
         is_dna = self.dna.atoms['resname'].isin(_dnaResidues)
@@ -1247,8 +1343,8 @@ def defineInteraction(self):
         for i, atom in atoms.iterrows():
             if atom.is_dna:
                 param = particle_definition.loc[atom['name']]
-                parameters = [param.epsilon * _ef,
-                              param.radius * _df]
+                parameters = [param.epsilon,
+                              param.radius]
             else:
                 parameters = [0, .1]  # Null energy and some radius)
             # print(i, parameters)
@@ -1258,7 +1354,7 @@ def defineInteraction(self):
         addNonBondedExclusions(self.dna, self.force)
 
 
-class Electrostatics(Force, simtk.openmm.CustomNonbondedForce):
+class Electrostatics(Force, openmm.CustomNonbondedForce):
     def __init__(self, dna, force_group=13, temperature=300*unit.kelvin, salt_concentration=100*unit.millimolar, OpenCLPatch=True):
         self.force_group = force_group
         self.T = temperature
@@ -1284,7 +1380,7 @@ def reset(self):
         denominator = denominator.in_units_of(unit.kilocalorie_per_mole**-1 * unit.nanometer**-1)
         #print(ldby, denominator)
 
-        electrostaticForce = simtk.openmm.CustomNonbondedForce("""energy;
+        electrostaticForce = openmm.CustomNonbondedForce("""energy;
                                                                 energy=q1*q2*exp(-r/dh_length)/denominator/r;""")
         electrostaticForce.addPerParticleParameter('q')
         electrostaticForce.addGlobalParameter('dh_length', ldby)
@@ -1335,13 +1431,13 @@ def __init__(self, dna, protein, k=1, force_group=14):
         super().__init__(dna, protein)
 
     def reset(self):
-        k = self.k
-        exclusionForce = simtk.openmm.CustomNonbondedForce(f"""{k}*energy;
+        exclusionForce = openmm.CustomNonbondedForce(f"""k_exclusion_protein_DNA*energy;
                          energy=(4*epsilon*((sigma/r)^12-(sigma/r)^6)-offset)*step(cutoff-r);
                          offset=4*epsilon*((sigma/cutoff)^12-(sigma/cutoff)^6);
                          sigma=0.5*(sigma1+sigma2); 
                          epsilon=sqrt(epsilon1*epsilon2);
                          cutoff=sqrt(cutoff1*cutoff2)""")
+        exclusionForce.addGlobalParameter('k_exclusion_protein_DNA', self.k)
         exclusionForce.addPerParticleParameter('epsilon')
         exclusionForce.addPerParticleParameter('sigma')
         exclusionForce.addPerParticleParameter('cutoff')
@@ -1379,15 +1475,15 @@ def defineInteraction(self):
         for i, atom in atoms.iterrows():
             if atom.is_dna:
                 param = particle_definition.loc['DNA' + atom['name']]
-                parameters = [param.epsilon * _ef,
-                              param.radius * _df,
-                              param.cutoff * _df]
+                parameters = [param.epsilon,
+                              param.radius,
+                              param.cutoff]
                 DNA_list += [i]
             elif atom.is_protein:
                 param = particle_definition.loc['Protein' + atom['name']]
-                parameters = [param.epsilon * _ef,
-                              param.radius * _df,
-                              param.cutoff * _df]
+                parameters = [param.epsilon,
+                              param.radius,
+                              param.cutoff]
                 protein_list += [i]
             else:
                 print(f'Residue {i} not included in protein-DNA interactions')
@@ -1421,7 +1517,7 @@ def reset(self):
         denominator = denominator.in_units_of(unit.kilocalorie_per_mole**-1 * unit.nanometer**-1)
         #print(ldby, denominator)
         k = self.k
-        electrostaticForce = simtk.openmm.CustomNonbondedForce(f"""k_electro_protein_DNA*energy;
+        electrostaticForce = openmm.CustomNonbondedForce(f"""k_electro_protein_DNA*energy;
                              energy=q1*q2*exp(-r/inter_dh_length)/inter_denominator/r;""")
         electrostaticForce.addPerParticleParameter('q')
         electrostaticForce.addGlobalParameter('k_electro_protein_DNA', k)
@@ -1488,6 +1584,7 @@ def defineInteraction(self):
         # addExclusions
         addNonBondedExclusions(self.dna, self.force)
 
+
 class AMHgoProteinDNA(ProteinDNAForce):
     """ Protein-DNA amhgo potential"""
     def __init__(self, dna, protein, chain_protein='A', chain_DNA='B', k_amhgo_PD=1*unit.kilocalorie_per_mole, sigma_sq=0.05*unit.nanometers**2, aaweight=False, globalct=True, cutoff=1.8, force_group=16):
@@ -1505,9 +1602,9 @@ def reset(self):
         cutoff = self.cutoff
         k_3spn2 = self.k_3spn2
         if self.globalct:
-                amhgoForce = simtk.openmm.CustomBondForce(f"-k_amhgo_PD*gamma_ij*exp(-(r-r_ijN)^2/(2*sigma_sq))*step({cutoff}-r)")
+                amhgoForce = openmm.CustomBondForce(f"-k_amhgo_PD*gamma_ij*exp(-(r-r_ijN)^2/(2*sigma_sq))*step({cutoff}-r)")
         else:
-                amhgoForce = simtk.openmm.CustomBondForce(f"-k_amhgo_PD*gamma_ij*exp(-(r-r_ijN)^2/(2*sigma_sq))*step(r_ijN+{cutoff}-r)")
+                amhgoForce = openmm.CustomBondForce(f"-k_amhgo_PD*gamma_ij*exp(-(r-r_ijN)^2/(2*sigma_sq))*step(r_ijN+{cutoff}-r)")
         amhgoForce.addGlobalParameter("k_amhgo_PD", k_3spn2*self.k_amhgo_PD)
         amhgoForce.addGlobalParameter("sigma_sq", self.sigma_sq)
         amhgoForce.addPerBondParameter("gamma_ij")
@@ -1523,7 +1620,10 @@ def defineInteraction(self):
 
         contact_list = np.loadtxt("contact_protein_DNA.dat")
         for i in range(len(contact_list)):
-            gamma_ij = contact_list[i][3] if self.aaweight else 1.0
+            if self.aaweight:
+                gamma_ij = contact_list[i][3]
+            else:
+                gamma_ij = 1.0
             if (self.chain_protein, int(contact_list[i][0]), 'CB') in atoms.index:
                  CB_protein = atoms[(atoms['chainID'] == self.chain_protein) & (atoms['resSeq'] == int(contact_list[i][0])) & (atoms['name'] == 'CB') & atoms['resname'].isin(_proteinResidues)].copy()
             else:
@@ -1531,10 +1631,10 @@ def defineInteraction(self):
             base_DNA = atoms[(atoms['chainID'] == self.chain_DNA) & (atoms['resSeq'] == int(contact_list[i][1])) & (atoms['name'].isin(['A', 'T', 'G', 'C'])) & atoms['resname'].isin(_dnaResidues)].copy()
             r_ijN = contact_list[i][2]/10.0*unit.nanometers
             self.force.addBond(int(CB_protein['index'].values[0]), int(base_DNA['index'].values[0]), [gamma_ij, r_ijN])
-            #print(int(CB_protein['index'].values[0]), int(base_DNA['index'].values[0]), [gamma_ij, r_ijN])
+            print(int(CB_protein['index'].values[0]), int(base_DNA['index'].values[0]), [gamma_ij, r_ijN])
 
 
-# class AMHgoProteinDNA(ProteinDNAForce):
+#class AMHgoProteinDNA(ProteinDNAForce):
 #    """ Protein-DNA amhgo potential (Xinyu)"""
 #    def __init__(self, dna, protein, chain_protein='A', chain_DNA='B', k_amhgo_PD=1*unit.kilocalorie_per_mole,
 #                 sigma_sq=0.05*unit.nanometers**2, aaweight=False, cutoff=1.8, force_group=16):
@@ -1549,7 +1649,7 @@ def defineInteraction(self):
 #
 #    def reset(self):
 #        cutoff = self.cutoff
-#        amhgoForce = simtk.openmm.CustomBondForce(f"-k_amhgo_PD*gamma_ij*exp(-(r-r_ijN)^2/(2*sigma_sq))*step({cutoff}-r)")
+#        amhgoForce = openmm.CustomBondForce(f"-k_amhgo_PD*gamma_ij*exp(-(r-r_ijN)^2/(2*sigma_sq))*step({cutoff}-r)")
 #        amhgoForce.addGlobalParameter("k_amhgo_PD", self.k_amhgo_PD)
 #        amhgoForce.addGlobalParameter("sigma_sq", self.sigma_sq)
 #        amhgoForce.addPerBondParameter("gamma_ij")
@@ -1573,8 +1673,7 @@ def defineInteraction(self):
 #            else:
 #                 CB_protein = atoms[(atoms['chainID'] == self.chain_protein) & (atoms['resSeq'] == int(contact_list[i][0])) &
 #                                    (atoms['name'] == 'CA') & atoms['resname'].isin(_proteinResidues)].copy()
-#            base_DNA = atoms[(atoms['chainID'] == self.chain_DNA) & (atoms['resSeq'] == int(contact_list[i][1])) &
-#                             (atoms['name'].isin(['A', 'T', 'G', 'C'])) & atoms['resname'].isin(_dnaResidues)].copy()
+#            base_DNA = atoms[(atoms['chainID'] == self.chain_DNA) & (atoms['resSeq'] == int(contact_list[i][1])) & (atoms['name'].isin(['A', 'T', 'G', 'C'])) & atoms['resname'].isin(_dnaResidues)].copy()
 #            r_ijN = contact_list[i][2]/10.0*unit.nanometers
 #            self.force.addBond(int(CB_protein['index'].values[0]), int(base_DNA['index'].values[0]), [gamma_ij, r_ijN])
 #            print(int(CB_protein['index'].values[0]), int(base_DNA['index'].values[0]), [gamma_ij, r_ijN])
@@ -1593,7 +1692,7 @@ def __init__(self, dna, protein, r0, chain_protein='A', chain_DNA='B', k_string_
     def reset(self):
         r0=self.r0
         k_string_PD=self.k_string_PD
-        stringForce = simtk.openmm.CustomCentroidBondForce(2, f"0.5*{k_string_PD}*(distance(g1,g2)-{r0})^2")
+        stringForce = openmm.CustomCentroidBondForce(2, f"0.5*{k_string_PD}*(distance(g1,g2)-{r0})^2")
         self.force = stringForce
         print("String_PD bias on: r0, k_string = ", r0, k_string_PD)
 
@@ -1624,7 +1723,7 @@ def __init__(self, dna, protein, chain_protein='A', chain_DNA='B', protein_seg=F
         super().__init__(dna, protein)
 
     def reset(self):
-        length = simtk.openmm.CustomCentroidBondForce(2, "distance(g1,g2)")
+        length = openmm.CustomCentroidBondForce(2, "distance(g1,g2)")
         length.setForceGroup(self.force_group)
         self.force = length
 
@@ -1657,77 +1756,6 @@ def defineInteraction(self):
                         ElectrostaticsProteinDNA=ElectrostaticsProteinDNA)
 
 
-# Unit testing
-def test_DNA_from_pdb():
-    """ Test correct DNA initialization from PDB"""
-    mol = DNA.fromPDB("Tests/1svc/1svc.pdb", template_from_X3DNA=False)
-
-
-def test_DNA_from_gro():
-    """ Test correct DNA initialization from gromacs files"""
-    pass
-
-
-def test_DNA_from_seq():
-    """ Test correct DNA initialization from sequence files"""
-    return True #Needs X3DNA
-    seq = 'ATACAAAGGTGCGAGGTTTCTATGCTCCCACG'
-    dna = DNA.fromSequence(seq, dna_type='B_curved')
-
-    # Compute the topology for the DNA structure.
-    # Since the dna was generated from the sequence using X3DNA,
-    # it is not necesary to recompute the geometry.
-
-    dna.computeTopology(template_from_X3DNA=False)
-
-    # Create the system.
-    # To set periodic boundary conditions (periodicBox=[50,50,50]).
-    # The periodic box size is in nanometers.
-    dna.periodic = False
-    s = System(dna, periodicBox=None)
-
-    # Add 3SPN2 forces
-    s.add3SPN2forces(verbose=True)
-
-    import simtk.openmm
-    import simtk.openmm.app
-    import simtk.unit
-    import sys
-    import numpy as np
-
-    # Initialize Molecular Dynamics simulations
-    s.initializeMD(temperature=300 * simtk.unit.kelvin, platform_name='OpenCL')
-    simulation = s.simulation
-
-    # Set initial positions
-    simulation.context.setPositions(s.coord.getPositions())
-
-    energy_unit = simtk.openmm.unit.kilojoule_per_mole
-    # Total energy
-    state = simulation.context.getState(getEnergy=True)
-    energy = state.getPotentialEnergy().value_in_unit(energy_unit)
-    print('TotalEnergy', round(energy, 6), energy_unit.get_symbol())
-
-    # Detailed energy
-    energies = {}
-    for force_name, force in s.forces.items():
-        group = force.getForceGroup()
-        state = simulation.context.getState(getEnergy=True, groups=2 ** group)
-        energies[force_name] = state.getPotentialEnergy().value_in_unit(energy_unit)
-
-    for force_name in s.forces.keys():
-        print(force_name, round(energies[force_name], 6), energy_unit.get_symbol())
-
-
-def test_DNA_from_xyz():
-    """Tests the correct parsing from an xyz file"""
-    mol = DNA.fromXYZ('Tests/adna/in00_conf.xyz', template_from_X3DNA=False)
-    assert mol.atoms.at[8, 'name'] == 'P'
-    assert round(mol.atoms.at[188, 'y'], 6) == -8.779343
-
-
-# Functional testing
-
 def parse_xyz(filename=''):
     columns = ['N', 'timestep', 'id', 'name', 'x', 'y', 'z']
     data = []
@@ -1781,125 +1809,6 @@ def parse_log(filename=''):
     return log_data
 
 
-def test_parse_xyz():
-    """Tests the example trajectory parsing"""
-    xyz_data = parse_xyz('Tests/adna/traj.xyz')
-    assert xyz_data.at[1, 'name'] == 7
-    assert xyz_data.at[1, 'x'] == 4.34621
-
-
-def test_parse_log():
-    """Tests the example log parsing"""
-    log_data = parse_log('Tests/adna/sim.log')
-    assert log_data.at[1, 'Step'] == 2000
-    assert log_data.at[1, 'eexcl'] == 0.45734636
-
-
-class TestEnergies:
-    """Tests that the energies are the same as the example outputs from lammps"""
-
-    def _test_energy(self,
-                     log_energy='E_bond',
-                     log_file='Tests/adna/sim.log',
-                     traj_file='Tests/adna/traj.xyz',
-                     force='Bond', periodic_size=94.2,
-                     platform_name='Reference', dna=None, system=None):
-        self.dna = dna
-        self.system = system
-        self.system.clearForces()
-        self.system.setDefaultPeriodicBoxVectors(*np.diag([2 * periodic_size / 10] * 3))
-
-        log = parse_log(log_file)
-
-        self.system.clearForces()
-        f = forces[force]
-        tempforce = f(self.dna)
-        try:
-            tempforce.addForce(self.system)
-        except AttributeError:
-            self.system.addForce(tempforce)
-        energies = self.system.recomputeEnergy(traj_file, platform_name=platform_name)
-        d = (energies / _ef - log[log_energy])
-        diff = np.sqrt((d ** 2).sum() / len(energies))
-        print(f'The difference in the energy term {log_energy} is {diff} Kcal/mol')
-        print(f'The DNA type of the system is {self.dna.DNAtype}')
-        data = np.array([energies / _ef, np.array(log[log_energy]), np.array(d)])
-        results = pandas.DataFrame(data.T, columns=['Openmm energy', 'Lammps energy', 'Difference'])
-        # print(data)
-        print(pandas.DataFrame(data.T, columns=['Openmm energy', 'Lammps energy', 'Difference']))
-        # The maximum error seems to be bonds in curved BDNA (0.002)
-        assert diff < 2E-3, diff
-        assert len(results.dropna()) == len(results), results
-
-    def test_energies(self):
-        test_sets = pandas.read_csv('Tests/test_cases.csv', comment='#')
-        for i, tests in test_sets.groupby(['Folder', 'DNA type']):
-            folder = i[0]
-            dna_type = i[1]
-            self.dna = DNA.fromXYZ(f'{folder}/in00_conf.xyz', dna_type, template_from_X3DNA=False)
-            self.system = System(self.dna)
-            for j, test in tests.iterrows():
-                print(j)
-                yield self._test_energy, test['Energy term'], f'{folder}/{test.Log}', f'{folder}/{test.Trajectory}', \
-                      test['Name'], test['periodic size'], test['Platform'], self.dna, self.system
-
-    def _test_force(self,
-                    log_energy='E_bond',
-                    log_file='Tests/adna/sim.log',
-                    traj_file='Tests/adna/traj.xyz',
-                    force='Bond', periodic_size=94.2,
-                    platform_name='Reference', dna=None, system=None):
-        self.dna = dna
-        self.system = system
-        self.system.clearForces()
-        self.system.setDefaultPeriodicBoxVectors(*np.diag([2 * periodic_size / 10] * 3))
-        log = parse_log(log_file)
-        self.system.clearForces()
-        f = forces[force]
-        tempforce = f(self.dna)
-        try:
-            tempforce.addForce(self.system)
-        except AttributeError:
-            self.system.addForce(tempforce)
-        temperature = 300 * simtk.openmm.unit.kelvin
-        integrator = simtk.openmm.LangevinIntegrator(temperature, 1 / simtk.openmm.unit.picosecond,
-                                                     2 * simtk.openmm.unit.femtoseconds)
-        platform = simtk.openmm.Platform.getPlatformByName(platform_name)
-        simulation = simtk.openmm.app.Simulation(self.system.top, self.system, integrator, platform)
-        simulation.context.setPositions(self.system.coord.getPositions())
-        energy_unit = simtk.openmm.unit.kilojoule_per_mole
-        nan_force_particles = 0
-        for i in range(10):
-            state = simulation.context.getState(getForces=True)
-            ff = state.getForces()
-            sf = (np.array(ff) ** 2).sum(axis=1) ** .5
-            for j, f in enumerate(sf):
-                if np.isnan(f.value_in_unit(unit.kilojoule_per_mole / unit.nanometer)):
-                    print(f"Particle {j + 1}/{len(sf)} has force {f} at step {i}")
-                    nan_force_particles += 1
-            simulation.step(1)
-        assert nan_force_particles == 0, "At least one particle has undefined force"
-
-    def test_forces(self):
-        test_sets = pandas.read_csv('Tests/test_cases.csv', comment='#')
-        for i, tests in test_sets.groupby(['Folder', 'DNA type']):
-            folder = i[0]
-            dna_type = i[1]
-            self.dna = DNA.fromXYZ(f'{folder}/in00_conf.xyz', dna_type, template_from_X3DNA=False)
-            self.system = System(self.dna)
-            for j, test in tests.iterrows():
-                print(j)
-                yield self._test_force, test['Energy term'], f'{folder}/{test.Log}', f'{folder}/{test.Trajectory}', \
-                      test['Name'], test['periodic size'], test['Platform'], self.dna, self.system
-
-    def _test_energies_slow(self):
-        test_sets = pandas.read_csv('Tests/test_cases.csv', comment='#')
-        for i, test in test_sets.iterrows():
-            dna_type = test['DNA type']
-            folder = test['Folder']
-            yield self._test_energy, test['Energy term'], f'{folder}/{test.Log}', f'{folder}/{test.Trajectory}', \
-                  test['Name'], folder, dna_type, test['periodic size'], test['Platform']
-
 # Some typical errors
 class BaseError(Exception):
     pass
diff --git a/open3SPN2/function_terms.py b/open3SPN2/function_terms.py
index 174cad2..333c8eb 100644
--- a/open3SPN2/function_terms.py
+++ b/open3SPN2/function_terms.py
@@ -1,6 +1,6 @@
 
-from simtk.openmm.app import *
-from simtk.openmm import *
+from openmm.app import *
+from openmm import *
 from simtk.unit import *
 import numpy as np
 
diff --git a/open3SPN2/geometry.py b/open3SPN2/geometry.py
new file mode 100644
index 0000000..6df0fcd
--- /dev/null
+++ b/open3SPN2/geometry.py
@@ -0,0 +1,107 @@
+import numpy as np
+
+class Transform:
+    def __init__(self, dx=0, dy=0, dz=0, rx=0, ry=0, rz=0):
+        """
+        Initialize the Transform object with displacements and rotations.
+        
+        Parameters:
+        - dx, dy, dz: Displacements in the x, y, z directions respectively (Angstroms).
+        - rx, ry, rz: Rotations around the x, y, z axes respectively (Radians).
+        """
+        self.dx = dx
+        self.dy = dy
+        self.dz = dz
+        self.rx = rx
+        self.ry = ry
+        self.rz = rz
+        
+        # Angle needed to rotate around the z axis to achieve a single rotation around the y axis
+        self.phi = np.arctan2(self.rx, self.ry)
+        
+        # Magnitude of the combined x and y rotation
+        self.yx_ang = np.sqrt(self.rx**2 + self.ry**2)
+
+    @property
+    def full_rotation(self):
+        """
+        Compute the full rotation transformation matrix that transitions from one state to another.
+        
+        Returns:
+        - A 3x3 rotation matrix.
+        """
+        rotz1 = self.rotz(0.5 * self.rz - self.phi)
+        rotyx = self.roty(self.yx_ang)
+        rotz2 = self.rotz(0.5 * self.rz + self.phi)
+        return rotz1 @ rotyx @ rotz2
+
+    @property
+    def half_rotation(self):
+        """
+        Compute the rotation matrix for half the rotation.
+        
+        This aims to rotate two objects equally relative to a reference frame.
+        
+        Returns:
+        - A 3x3 rotation matrix.
+        """
+        rotz1 = self.rotz(0.5 * self.rz - self.phi)
+        rotyx = self.roty(0.5 * self.yx_ang)
+        rotz2 = self.rotz(self.phi)
+        return rotz1 @ rotyx @ rotz2
+
+    @property
+    def full_displacement(self):
+        """
+        Compute the displacement vector for transitioning from one state to another.
+        
+        Returns:
+        - A 1x3 displacement vector.
+        """
+        displacement = np.array([self.dx, self.dy, self.dz])
+        return np.dot(displacement, self.half_rotation.T)
+
+    @property
+    def half_displacement(self):
+        """
+        Compute the midpoint of the displacement vector.
+        
+        Returns:
+        - A 1x3 displacement vector.
+        """
+        return 0.5 * self.full_displacement
+    
+    @staticmethod
+    def rotz(theta):
+        """
+        Compute the rotation matrix around the Z-axis.
+        
+        Parameters:
+        - theta: Rotation angle in radians.
+        
+        Returns:
+        - A 3x3 rotation matrix.
+        """
+        c, s = np.cos(theta), np.sin(theta)
+        return np.array([[c, -s, 0], [s, c, 0], [0, 0, 1]])
+
+    @staticmethod
+    def roty(theta):
+        """
+        Compute the rotation matrix around the Y-axis.
+        
+        Parameters:
+        - theta: Rotation angle in radians.
+        
+        Returns:
+        - A 3x3 rotation matrix.
+        """
+        c, s = np.cos(theta), np.sin(theta)
+        return np.array([[c, 0, s], [0, 1, 0], [-s, 0, c]])
+
+
+
+
+if __name__=='__main__':
+    test_default_values()
+    test_given_values()
\ No newline at end of file
diff --git a/setup.py b/setup.py
index 52f8ab5..e700c5e 100644
--- a/setup.py
+++ b/setup.py
@@ -29,7 +29,7 @@
       url = 'https://github.com/cabb99/open3spn2', 
       download_url = 'https://github.com/cabb99/open3spn2/archive/0.2.0.tar.gz',
       keywords = ['dna', 'forcefield', 'openmm'],
-      package_data={'open3SPN2': ['3SPN2.conf','3SPN2.xml']},
+      package_data={'open3SPN2': ['3SPN2.conf','3SPN2.xml','DNA_atomic.csv']},
       include_package_data=True,
       install_requires=[
               'biopython',
diff --git a/versioneer.py b/versioneer.py
index 64fea1c..1e3753e 100644
--- a/versioneer.py
+++ b/versioneer.py
@@ -1,5 +1,5 @@
 
-# Version: 0.18
+# Version: 0.29
 
 """The Versioneer - like a rocketeer, but for versions.
 
@@ -7,18 +7,14 @@
 ==============
 
 * like a rocketeer, but for versions!
-* https://github.com/warner/python-versioneer
+* https://github.com/python-versioneer/python-versioneer
 * Brian Warner
-* License: Public Domain
-* Compatible With: python2.6, 2.7, 3.2, 3.3, 3.4, 3.5, 3.6, and pypy
-* [![Latest Version]
-(https://pypip.in/version/versioneer/badge.svg?style=flat)
-](https://pypi.python.org/pypi/versioneer/)
-* [![Build Status]
-(https://travis-ci.org/warner/python-versioneer.png?branch=master)
-](https://travis-ci.org/warner/python-versioneer)
-
-This is a tool for managing a recorded version number in distutils-based
+* License: Public Domain (Unlicense)
+* Compatible with: Python 3.7, 3.8, 3.9, 3.10, 3.11 and pypy3
+* [![Latest Version][pypi-image]][pypi-url]
+* [![Build Status][travis-image]][travis-url]
+
+This is a tool for managing a recorded version number in setuptools-based
 python projects. The goal is to remove the tedious and error-prone "update
 the embedded version string" step from your release process. Making a new
 release should be as easy as recording a new tag in your version-control
@@ -27,9 +23,38 @@
 
 ## Quick Install
 
-* `pip install versioneer` to somewhere to your $PATH
-* add a `[versioneer]` section to your setup.cfg (see below)
-* run `versioneer install` in your source tree, commit the results
+Versioneer provides two installation modes. The "classic" vendored mode installs
+a copy of versioneer into your repository. The experimental build-time dependency mode
+is intended to allow you to skip this step and simplify the process of upgrading.
+
+### Vendored mode
+
+* `pip install versioneer` to somewhere in your $PATH
+   * A [conda-forge recipe](https://github.com/conda-forge/versioneer-feedstock) is
+     available, so you can also use `conda install -c conda-forge versioneer`
+* add a `[tool.versioneer]` section to your `pyproject.toml` or a
+  `[versioneer]` section to your `setup.cfg` (see [Install](INSTALL.md))
+   * Note that you will need to add `tomli; python_version < "3.11"` to your
+     build-time dependencies if you use `pyproject.toml`
+* run `versioneer install --vendor` in your source tree, commit the results
+* verify version information with `python setup.py version`
+
+### Build-time dependency mode
+
+* `pip install versioneer` to somewhere in your $PATH
+   * A [conda-forge recipe](https://github.com/conda-forge/versioneer-feedstock) is
+     available, so you can also use `conda install -c conda-forge versioneer`
+* add a `[tool.versioneer]` section to your `pyproject.toml` or a
+  `[versioneer]` section to your `setup.cfg` (see [Install](INSTALL.md))
+* add `versioneer` (with `[toml]` extra, if configuring in `pyproject.toml`)
+  to the `requires` key of the `build-system` table in `pyproject.toml`:
+  ```toml
+  [build-system]
+  requires = ["setuptools", "versioneer[toml]"]
+  build-backend = "setuptools.build_meta"
+  ```
+* run `versioneer install --no-vendor` in your source tree, commit the results
+* verify version information with `python setup.py version`
 
 ## Version Identifiers
 
@@ -61,7 +86,7 @@
 for example `git describe --tags --dirty --always` reports things like
 "0.7-1-g574ab98-dirty" to indicate that the checkout is one revision past the
 0.7 tag, has a unique revision id of "574ab98", and is "dirty" (it has
-uncommitted changes.
+uncommitted changes).
 
 The version identifier is used for multiple purposes:
 
@@ -166,7 +191,7 @@
 
 Some situations are known to cause problems for Versioneer. This details the
 most significant ones. More can be found on Github
-[issues page](https://github.com/warner/python-versioneer/issues).
+[issues page](https://github.com/python-versioneer/python-versioneer/issues).
 
 ### Subprojects
 
@@ -180,7 +205,7 @@
   `setup.cfg`, and `tox.ini`. Projects like these produce multiple PyPI
   distributions (and upload multiple independently-installable tarballs).
 * Source trees whose main purpose is to contain a C library, but which also
-  provide bindings to Python (and perhaps other langauges) in subdirectories.
+  provide bindings to Python (and perhaps other languages) in subdirectories.
 
 Versioneer will look for `.git` in parent directories, and most operations
 should get the right version string. However `pip` and `setuptools` have bugs
@@ -194,9 +219,9 @@
 Pip-8.1.1 is known to have this problem, but hopefully it will get fixed in
 some later version.
 
-[Bug #38](https://github.com/warner/python-versioneer/issues/38) is tracking
+[Bug #38](https://github.com/python-versioneer/python-versioneer/issues/38) is tracking
 this issue. The discussion in
-[PR #61](https://github.com/warner/python-versioneer/pull/61) describes the
+[PR #61](https://github.com/python-versioneer/python-versioneer/pull/61) describes the
 issue from the Versioneer side in more detail.
 [pip PR#3176](https://github.com/pypa/pip/pull/3176) and
 [pip PR#3615](https://github.com/pypa/pip/pull/3615) contain work to improve
@@ -224,31 +249,20 @@
 cause egg_info to be rebuilt (including `sdist`, `wheel`, and installing into
 a different virtualenv), so this can be surprising.
 
-[Bug #83](https://github.com/warner/python-versioneer/issues/83) describes
+[Bug #83](https://github.com/python-versioneer/python-versioneer/issues/83) describes
 this one, but upgrading to a newer version of setuptools should probably
 resolve it.
 
-### Unicode version strings
-
-While Versioneer works (and is continually tested) with both Python 2 and
-Python 3, it is not entirely consistent with bytes-vs-unicode distinctions.
-Newer releases probably generate unicode version strings on py2. It's not
-clear that this is wrong, but it may be surprising for applications when then
-write these strings to a network connection or include them in bytes-oriented
-APIs like cryptographic checksums.
-
-[Bug #71](https://github.com/warner/python-versioneer/issues/71) investigates
-this question.
-
 
 ## Updating Versioneer
 
 To upgrade your project to a new release of Versioneer, do the following:
 
 * install the new Versioneer (`pip install -U versioneer` or equivalent)
-* edit `setup.cfg`, if necessary, to include any new configuration settings
-  indicated by the release notes. See [UPGRADING](./UPGRADING.md) for details.
-* re-run `versioneer install` in your source tree, to replace
+* edit `setup.cfg` and `pyproject.toml`, if necessary,
+  to include any new configuration settings indicated by the release notes.
+  See [UPGRADING](./UPGRADING.md) for details.
+* re-run `versioneer install --[no-]vendor` in your source tree, to replace
   `SRC/_version.py`
 * commit any changed files
 
@@ -265,35 +279,70 @@
 direction and include code from all supported VCS systems, reducing the
 number of intermediate scripts.
 
+## Similar projects
+
+* [setuptools_scm](https://github.com/pypa/setuptools_scm/) - a non-vendored build-time
+  dependency
+* [minver](https://github.com/jbweston/miniver) - a lightweight reimplementation of
+  versioneer
+* [versioningit](https://github.com/jwodder/versioningit) - a PEP 518-based setuptools
+  plugin
 
 ## License
 
 To make Versioneer easier to embed, all its code is dedicated to the public
 domain. The `_version.py` that it creates is also in the public domain.
-Specifically, both are released under the Creative Commons "Public Domain
-Dedication" license (CC0-1.0), as described in
-https://creativecommons.org/publicdomain/zero/1.0/ .
+Specifically, both are released under the "Unlicense", as described in
+https://unlicense.org/.
+
+[pypi-image]: https://img.shields.io/pypi/v/versioneer.svg
+[pypi-url]: https://pypi.python.org/pypi/versioneer/
+[travis-image]:
+https://img.shields.io/travis/com/python-versioneer/python-versioneer.svg
+[travis-url]: https://travis-ci.com/github/python-versioneer/python-versioneer
 
 """
+# pylint:disable=invalid-name,import-outside-toplevel,missing-function-docstring
+# pylint:disable=missing-class-docstring,too-many-branches,too-many-statements
+# pylint:disable=raise-missing-from,too-many-lines,too-many-locals,import-error
+# pylint:disable=too-few-public-methods,redefined-outer-name,consider-using-with
+# pylint:disable=attribute-defined-outside-init,too-many-arguments
 
-from __future__ import print_function
-try:
-    import configparser
-except ImportError:
-    import ConfigParser as configparser
+import configparser
 import errno
 import json
 import os
 import re
 import subprocess
 import sys
+from pathlib import Path
+from typing import Any, Callable, cast, Dict, List, Optional, Tuple, Union
+from typing import NoReturn
+import functools
+
+have_tomllib = True
+if sys.version_info >= (3, 11):
+    import tomllib
+else:
+    try:
+        import tomli as tomllib
+    except ImportError:
+        have_tomllib = False
 
 
 class VersioneerConfig:
     """Container for Versioneer configuration parameters."""
 
+    VCS: str
+    style: str
+    tag_prefix: str
+    versionfile_source: str
+    versionfile_build: Optional[str]
+    parentdir_prefix: Optional[str]
+    verbose: Optional[bool]
+
 
-def get_root():
+def get_root() -> str:
     """Get the project root directory.
 
     We require that all commands are run from the project root, i.e. the
@@ -301,13 +350,23 @@ def get_root():
     """
     root = os.path.realpath(os.path.abspath(os.getcwd()))
     setup_py = os.path.join(root, "setup.py")
+    pyproject_toml = os.path.join(root, "pyproject.toml")
     versioneer_py = os.path.join(root, "versioneer.py")
-    if not (os.path.exists(setup_py) or os.path.exists(versioneer_py)):
+    if not (
+        os.path.exists(setup_py)
+        or os.path.exists(pyproject_toml)
+        or os.path.exists(versioneer_py)
+    ):
         # allow 'python path/to/setup.py COMMAND'
         root = os.path.dirname(os.path.realpath(os.path.abspath(sys.argv[0])))
         setup_py = os.path.join(root, "setup.py")
+        pyproject_toml = os.path.join(root, "pyproject.toml")
         versioneer_py = os.path.join(root, "versioneer.py")
-    if not (os.path.exists(setup_py) or os.path.exists(versioneer_py)):
+    if not (
+        os.path.exists(setup_py)
+        or os.path.exists(pyproject_toml)
+        or os.path.exists(versioneer_py)
+    ):
         err = ("Versioneer was unable to run the project root directory. "
                "Versioneer requires setup.py to be executed from "
                "its immediate directory (like 'python setup.py COMMAND'), "
@@ -321,43 +380,62 @@ def get_root():
         # module-import table will cache the first one. So we can't use
         # os.path.dirname(__file__), as that will find whichever
         # versioneer.py was first imported, even in later projects.
-        me = os.path.realpath(os.path.abspath(__file__))
-        me_dir = os.path.normcase(os.path.splitext(me)[0])
+        my_path = os.path.realpath(os.path.abspath(__file__))
+        me_dir = os.path.normcase(os.path.splitext(my_path)[0])
         vsr_dir = os.path.normcase(os.path.splitext(versioneer_py)[0])
-        if me_dir != vsr_dir:
+        if me_dir != vsr_dir and "VERSIONEER_PEP518" not in globals():
             print("Warning: build in %s is using versioneer.py from %s"
-                  % (os.path.dirname(me), versioneer_py))
+                  % (os.path.dirname(my_path), versioneer_py))
     except NameError:
         pass
     return root
 
 
-def get_config_from_root(root):
+def get_config_from_root(root: str) -> VersioneerConfig:
     """Read the project setup.cfg file to determine Versioneer config."""
-    # This might raise EnvironmentError (if setup.cfg is missing), or
+    # This might raise OSError (if setup.cfg is missing), or
     # configparser.NoSectionError (if it lacks a [versioneer] section), or
     # configparser.NoOptionError (if it lacks "VCS="). See the docstring at
     # the top of versioneer.py for instructions on writing your setup.cfg .
-    setup_cfg = os.path.join(root, "setup.cfg")
-    parser = configparser.SafeConfigParser()
-    with open(setup_cfg, "r") as f:
-        parser.readfp(f)
-    VCS = parser.get("versioneer", "VCS")  # mandatory
-
-    def get(parser, name):
-        if parser.has_option("versioneer", name):
-            return parser.get("versioneer", name)
-        return None
+    root_pth = Path(root)
+    pyproject_toml = root_pth / "pyproject.toml"
+    setup_cfg = root_pth / "setup.cfg"
+    section: Union[Dict[str, Any], configparser.SectionProxy, None] = None
+    if pyproject_toml.exists() and have_tomllib:
+        try:
+            with open(pyproject_toml, 'rb') as fobj:
+                pp = tomllib.load(fobj)
+            section = pp['tool']['versioneer']
+        except (tomllib.TOMLDecodeError, KeyError) as e:
+            print(f"Failed to load config from {pyproject_toml}: {e}")
+            print("Try to load it from setup.cfg")
+    if not section:
+        parser = configparser.ConfigParser()
+        with open(setup_cfg) as cfg_file:
+            parser.read_file(cfg_file)
+        parser.get("versioneer", "VCS")  # raise error if missing
+
+        section = parser["versioneer"]
+
+    # `cast`` really shouldn't be used, but its simplest for the
+    # common VersioneerConfig users at the moment. We verify against
+    # `None` values elsewhere where it matters
+
     cfg = VersioneerConfig()
-    cfg.VCS = VCS
-    cfg.style = get(parser, "style") or ""
-    cfg.versionfile_source = get(parser, "versionfile_source")
-    cfg.versionfile_build = get(parser, "versionfile_build")
-    cfg.tag_prefix = get(parser, "tag_prefix")
-    if cfg.tag_prefix in ("''", '""'):
+    cfg.VCS = section['VCS']
+    cfg.style = section.get("style", "")
+    cfg.versionfile_source = cast(str, section.get("versionfile_source"))
+    cfg.versionfile_build = section.get("versionfile_build")
+    cfg.tag_prefix = cast(str, section.get("tag_prefix"))
+    if cfg.tag_prefix in ("''", '""', None):
         cfg.tag_prefix = ""
-    cfg.parentdir_prefix = get(parser, "parentdir_prefix")
-    cfg.verbose = get(parser, "verbose")
+    cfg.parentdir_prefix = section.get("parentdir_prefix")
+    if isinstance(section, configparser.SectionProxy):
+        # Make sure configparser translates to bool
+        cfg.verbose = section.getboolean("verbose")
+    else:
+        cfg.verbose = section.get("verbose")
+
     return cfg
 
 
@@ -366,37 +444,48 @@ class NotThisMethod(Exception):
 
 
 # these dictionaries contain VCS-specific tools
-LONG_VERSION_PY = {}
-HANDLERS = {}
+LONG_VERSION_PY: Dict[str, str] = {}
+HANDLERS: Dict[str, Dict[str, Callable]] = {}
 
 
-def register_vcs_handler(vcs, method):  # decorator
-    """Decorator to mark a method as the handler for a particular VCS."""
-    def decorate(f):
+def register_vcs_handler(vcs: str, method: str) -> Callable:  # decorator
+    """Create decorator to mark a method as the handler of a VCS."""
+    def decorate(f: Callable) -> Callable:
         """Store f in HANDLERS[vcs][method]."""
-        if vcs not in HANDLERS:
-            HANDLERS[vcs] = {}
-        HANDLERS[vcs][method] = f
+        HANDLERS.setdefault(vcs, {})[method] = f
         return f
     return decorate
 
 
-def run_command(commands, args, cwd=None, verbose=False, hide_stderr=False,
-                env=None):
+def run_command(
+    commands: List[str],
+    args: List[str],
+    cwd: Optional[str] = None,
+    verbose: bool = False,
+    hide_stderr: bool = False,
+    env: Optional[Dict[str, str]] = None,
+) -> Tuple[Optional[str], Optional[int]]:
     """Call the given command(s)."""
     assert isinstance(commands, list)
-    p = None
-    for c in commands:
+    process = None
+
+    popen_kwargs: Dict[str, Any] = {}
+    if sys.platform == "win32":
+        # This hides the console window if pythonw.exe is used
+        startupinfo = subprocess.STARTUPINFO()
+        startupinfo.dwFlags |= subprocess.STARTF_USESHOWWINDOW
+        popen_kwargs["startupinfo"] = startupinfo
+
+    for command in commands:
         try:
-            dispcmd = str([c] + args)
+            dispcmd = str([command] + args)
             # remember shell=False, so use git.cmd on windows, not just git
-            p = subprocess.Popen([c] + args, cwd=cwd, env=env,
-                                 stdout=subprocess.PIPE,
-                                 stderr=(subprocess.PIPE if hide_stderr
-                                         else None))
+            process = subprocess.Popen([command] + args, cwd=cwd, env=env,
+                                       stdout=subprocess.PIPE,
+                                       stderr=(subprocess.PIPE if hide_stderr
+                                               else None), **popen_kwargs)
             break
-        except EnvironmentError:
-            e = sys.exc_info()[1]
+        except OSError as e:
             if e.errno == errno.ENOENT:
                 continue
             if verbose:
@@ -407,26 +496,25 @@ def run_command(commands, args, cwd=None, verbose=False, hide_stderr=False,
         if verbose:
             print("unable to find command, tried %s" % (commands,))
         return None, None
-    stdout = p.communicate()[0].strip()
-    if sys.version_info[0] >= 3:
-        stdout = stdout.decode()
-    if p.returncode != 0:
+    stdout = process.communicate()[0].strip().decode()
+    if process.returncode != 0:
         if verbose:
             print("unable to run %s (error)" % dispcmd)
             print("stdout was %s" % stdout)
-        return None, p.returncode
-    return stdout, p.returncode
+        return None, process.returncode
+    return stdout, process.returncode
 
 
-LONG_VERSION_PY['git'] = '''
+LONG_VERSION_PY['git'] = r'''
 # This file helps to compute a version number in source trees obtained from
 # git-archive tarball (such as those provided by githubs download-from-tag
 # feature). Distribution tarballs (built by setup.py sdist) and build
 # directories (produced by setup.py build) will contain a much shorter file
 # that just contains the computed version number.
 
-# This file is released into the public domain. Generated by
-# versioneer-0.18 (https://github.com/warner/python-versioneer)
+# This file is released into the public domain.
+# Generated by versioneer-0.29
+# https://github.com/python-versioneer/python-versioneer
 
 """Git implementation of _version.py."""
 
@@ -435,9 +523,11 @@ def run_command(commands, args, cwd=None, verbose=False, hide_stderr=False,
 import re
 import subprocess
 import sys
+from typing import Any, Callable, Dict, List, Optional, Tuple
+import functools
 
 
-def get_keywords():
+def get_keywords() -> Dict[str, str]:
     """Get the keywords needed to look up the version information."""
     # these strings will be replaced by git during git-archive.
     # setup.py/versioneer.py will grep for the variable names, so they must
@@ -453,8 +543,15 @@ def get_keywords():
 class VersioneerConfig:
     """Container for Versioneer configuration parameters."""
 
+    VCS: str
+    style: str
+    tag_prefix: str
+    parentdir_prefix: str
+    versionfile_source: str
+    verbose: bool
+
 
-def get_config():
+def get_config() -> VersioneerConfig:
     """Create, populate and return the VersioneerConfig() object."""
     # these strings are filled in when 'setup.py versioneer' creates
     # _version.py
@@ -472,13 +569,13 @@ class NotThisMethod(Exception):
     """Exception raised if a method is not valid for the current scenario."""
 
 
-LONG_VERSION_PY = {}
-HANDLERS = {}
+LONG_VERSION_PY: Dict[str, str] = {}
+HANDLERS: Dict[str, Dict[str, Callable]] = {}
 
 
-def register_vcs_handler(vcs, method):  # decorator
-    """Decorator to mark a method as the handler for a particular VCS."""
-    def decorate(f):
+def register_vcs_handler(vcs: str, method: str) -> Callable:  # decorator
+    """Create decorator to mark a method as the handler of a VCS."""
+    def decorate(f: Callable) -> Callable:
         """Store f in HANDLERS[vcs][method]."""
         if vcs not in HANDLERS:
             HANDLERS[vcs] = {}
@@ -487,22 +584,35 @@ def decorate(f):
     return decorate
 
 
-def run_command(commands, args, cwd=None, verbose=False, hide_stderr=False,
-                env=None):
+def run_command(
+    commands: List[str],
+    args: List[str],
+    cwd: Optional[str] = None,
+    verbose: bool = False,
+    hide_stderr: bool = False,
+    env: Optional[Dict[str, str]] = None,
+) -> Tuple[Optional[str], Optional[int]]:
     """Call the given command(s)."""
     assert isinstance(commands, list)
-    p = None
-    for c in commands:
+    process = None
+
+    popen_kwargs: Dict[str, Any] = {}
+    if sys.platform == "win32":
+        # This hides the console window if pythonw.exe is used
+        startupinfo = subprocess.STARTUPINFO()
+        startupinfo.dwFlags |= subprocess.STARTF_USESHOWWINDOW
+        popen_kwargs["startupinfo"] = startupinfo
+
+    for command in commands:
         try:
-            dispcmd = str([c] + args)
+            dispcmd = str([command] + args)
             # remember shell=False, so use git.cmd on windows, not just git
-            p = subprocess.Popen([c] + args, cwd=cwd, env=env,
-                                 stdout=subprocess.PIPE,
-                                 stderr=(subprocess.PIPE if hide_stderr
-                                         else None))
+            process = subprocess.Popen([command] + args, cwd=cwd, env=env,
+                                       stdout=subprocess.PIPE,
+                                       stderr=(subprocess.PIPE if hide_stderr
+                                               else None), **popen_kwargs)
             break
-        except EnvironmentError:
-            e = sys.exc_info()[1]
+        except OSError as e:
             if e.errno == errno.ENOENT:
                 continue
             if verbose:
@@ -513,18 +623,20 @@ def run_command(commands, args, cwd=None, verbose=False, hide_stderr=False,
         if verbose:
             print("unable to find command, tried %%s" %% (commands,))
         return None, None
-    stdout = p.communicate()[0].strip()
-    if sys.version_info[0] >= 3:
-        stdout = stdout.decode()
-    if p.returncode != 0:
+    stdout = process.communicate()[0].strip().decode()
+    if process.returncode != 0:
         if verbose:
             print("unable to run %%s (error)" %% dispcmd)
             print("stdout was %%s" %% stdout)
-        return None, p.returncode
-    return stdout, p.returncode
+        return None, process.returncode
+    return stdout, process.returncode
 
 
-def versions_from_parentdir(parentdir_prefix, root, verbose):
+def versions_from_parentdir(
+    parentdir_prefix: str,
+    root: str,
+    verbose: bool,
+) -> Dict[str, Any]:
     """Try to determine the version from the parent directory name.
 
     Source tarballs conventionally unpack into a directory that includes both
@@ -533,15 +645,14 @@ def versions_from_parentdir(parentdir_prefix, root, verbose):
     """
     rootdirs = []
 
-    for i in range(3):
+    for _ in range(3):
         dirname = os.path.basename(root)
         if dirname.startswith(parentdir_prefix):
             return {"version": dirname[len(parentdir_prefix):],
                     "full-revisionid": None,
                     "dirty": False, "error": None, "date": None}
-        else:
-            rootdirs.append(root)
-            root = os.path.dirname(root)  # up a level
+        rootdirs.append(root)
+        root = os.path.dirname(root)  # up a level
 
     if verbose:
         print("Tried directories %%s but none started with prefix %%s" %%
@@ -550,41 +661,48 @@ def versions_from_parentdir(parentdir_prefix, root, verbose):
 
 
 @register_vcs_handler("git", "get_keywords")
-def git_get_keywords(versionfile_abs):
+def git_get_keywords(versionfile_abs: str) -> Dict[str, str]:
     """Extract version information from the given file."""
     # the code embedded in _version.py can just fetch the value of these
     # keywords. When used from setup.py, we don't want to import _version.py,
     # so we do it with a regexp instead. This function is not used from
     # _version.py.
-    keywords = {}
+    keywords: Dict[str, str] = {}
     try:
-        f = open(versionfile_abs, "r")
-        for line in f.readlines():
-            if line.strip().startswith("git_refnames ="):
-                mo = re.search(r'=\s*"(.*)"', line)
-                if mo:
-                    keywords["refnames"] = mo.group(1)
-            if line.strip().startswith("git_full ="):
-                mo = re.search(r'=\s*"(.*)"', line)
-                if mo:
-                    keywords["full"] = mo.group(1)
-            if line.strip().startswith("git_date ="):
-                mo = re.search(r'=\s*"(.*)"', line)
-                if mo:
-                    keywords["date"] = mo.group(1)
-        f.close()
-    except EnvironmentError:
+        with open(versionfile_abs, "r") as fobj:
+            for line in fobj:
+                if line.strip().startswith("git_refnames ="):
+                    mo = re.search(r'=\s*"(.*)"', line)
+                    if mo:
+                        keywords["refnames"] = mo.group(1)
+                if line.strip().startswith("git_full ="):
+                    mo = re.search(r'=\s*"(.*)"', line)
+                    if mo:
+                        keywords["full"] = mo.group(1)
+                if line.strip().startswith("git_date ="):
+                    mo = re.search(r'=\s*"(.*)"', line)
+                    if mo:
+                        keywords["date"] = mo.group(1)
+    except OSError:
         pass
     return keywords
 
 
 @register_vcs_handler("git", "keywords")
-def git_versions_from_keywords(keywords, tag_prefix, verbose):
+def git_versions_from_keywords(
+    keywords: Dict[str, str],
+    tag_prefix: str,
+    verbose: bool,
+) -> Dict[str, Any]:
     """Get version information from git keywords."""
-    if not keywords:
-        raise NotThisMethod("no keywords at all, weird")
+    if "refnames" not in keywords:
+        raise NotThisMethod("Short version file found")
     date = keywords.get("date")
     if date is not None:
+        # Use only the last line.  Previous lines may contain GPG signature
+        # information.
+        date = date.splitlines()[-1]
+
         # git-2.2.0 added "%%cI", which expands to an ISO-8601 -compliant
         # datestamp. However we prefer "%%ci" (which expands to an "ISO-8601
         # -like" string, which we must then edit to make compliant), because
@@ -597,11 +715,11 @@ def git_versions_from_keywords(keywords, tag_prefix, verbose):
         if verbose:
             print("keywords are unexpanded, not using")
         raise NotThisMethod("unexpanded keywords, not a git-archive tarball")
-    refs = set([r.strip() for r in refnames.strip("()").split(",")])
+    refs = {r.strip() for r in refnames.strip("()").split(",")}
     # starting in git-1.8.3, tags are listed as "tag: foo-1.0" instead of
     # just "foo-1.0". If we see a "tag: " prefix, prefer those.
     TAG = "tag: "
-    tags = set([r[len(TAG):] for r in refs if r.startswith(TAG)])
+    tags = {r[len(TAG):] for r in refs if r.startswith(TAG)}
     if not tags:
         # Either we're using git < 1.8.3, or there really are no tags. We use
         # a heuristic: assume all version tags have a digit. The old git %%d
@@ -610,7 +728,7 @@ def git_versions_from_keywords(keywords, tag_prefix, verbose):
         # between branches and tags. By ignoring refnames without digits, we
         # filter out many common branch names like "release" and
         # "stabilization", as well as "HEAD" and "master".
-        tags = set([r for r in refs if re.search(r'\d', r)])
+        tags = {r for r in refs if re.search(r'\d', r)}
         if verbose:
             print("discarding '%%s', no digits" %% ",".join(refs - tags))
     if verbose:
@@ -619,6 +737,11 @@ def git_versions_from_keywords(keywords, tag_prefix, verbose):
         # sorting will prefer e.g. "2.0" over "2.0rc1"
         if ref.startswith(tag_prefix):
             r = ref[len(tag_prefix):]
+            # Filter out refs that exactly match prefix or that don't start
+            # with a number once the prefix is stripped (mostly a concern
+            # when prefix is '')
+            if not re.match(r'\d', r):
+                continue
             if verbose:
                 print("picking %%s" %% r)
             return {"version": r,
@@ -634,7 +757,12 @@ def git_versions_from_keywords(keywords, tag_prefix, verbose):
 
 
 @register_vcs_handler("git", "pieces_from_vcs")
-def git_pieces_from_vcs(tag_prefix, root, verbose, run_command=run_command):
+def git_pieces_from_vcs(
+    tag_prefix: str,
+    root: str,
+    verbose: bool,
+    runner: Callable = run_command
+) -> Dict[str, Any]:
     """Get version from 'git describe' in the root of the source tree.
 
     This only gets called if the git-archive 'subst' keywords were *not*
@@ -645,8 +773,15 @@ def git_pieces_from_vcs(tag_prefix, root, verbose, run_command=run_command):
     if sys.platform == "win32":
         GITS = ["git.cmd", "git.exe"]
 
-    out, rc = run_command(GITS, ["rev-parse", "--git-dir"], cwd=root,
-                          hide_stderr=True)
+    # GIT_DIR can interfere with correct operation of Versioneer.
+    # It may be intended to be passed to the Versioneer-versioned project,
+    # but that should not change where we get our version from.
+    env = os.environ.copy()
+    env.pop("GIT_DIR", None)
+    runner = functools.partial(runner, env=env)
+
+    _, rc = runner(GITS, ["rev-parse", "--git-dir"], cwd=root,
+                   hide_stderr=not verbose)
     if rc != 0:
         if verbose:
             print("Directory %%s not under git control" %% root)
@@ -654,24 +789,57 @@ def git_pieces_from_vcs(tag_prefix, root, verbose, run_command=run_command):
 
     # if there is a tag matching tag_prefix, this yields TAG-NUM-gHEX[-dirty]
     # if there isn't one, this yields HEX[-dirty] (no NUM)
-    describe_out, rc = run_command(GITS, ["describe", "--tags", "--dirty",
-                                          "--always", "--long",
-                                          "--match", "%%s*" %% tag_prefix],
-                                   cwd=root)
+    describe_out, rc = runner(GITS, [
+        "describe", "--tags", "--dirty", "--always", "--long",
+        "--match", f"{tag_prefix}[[:digit:]]*"
+    ], cwd=root)
     # --long was added in git-1.5.5
     if describe_out is None:
         raise NotThisMethod("'git describe' failed")
     describe_out = describe_out.strip()
-    full_out, rc = run_command(GITS, ["rev-parse", "HEAD"], cwd=root)
+    full_out, rc = runner(GITS, ["rev-parse", "HEAD"], cwd=root)
     if full_out is None:
         raise NotThisMethod("'git rev-parse' failed")
     full_out = full_out.strip()
 
-    pieces = {}
+    pieces: Dict[str, Any] = {}
     pieces["long"] = full_out
     pieces["short"] = full_out[:7]  # maybe improved later
     pieces["error"] = None
 
+    branch_name, rc = runner(GITS, ["rev-parse", "--abbrev-ref", "HEAD"],
+                             cwd=root)
+    # --abbrev-ref was added in git-1.6.3
+    if rc != 0 or branch_name is None:
+        raise NotThisMethod("'git rev-parse --abbrev-ref' returned error")
+    branch_name = branch_name.strip()
+
+    if branch_name == "HEAD":
+        # If we aren't exactly on a branch, pick a branch which represents
+        # the current commit. If all else fails, we are on a branchless
+        # commit.
+        branches, rc = runner(GITS, ["branch", "--contains"], cwd=root)
+        # --contains was added in git-1.5.4
+        if rc != 0 or branches is None:
+            raise NotThisMethod("'git branch --contains' returned error")
+        branches = branches.split("\n")
+
+        # Remove the first line if we're running detached
+        if "(" in branches[0]:
+            branches.pop(0)
+
+        # Strip off the leading "* " from the list of branches.
+        branches = [branch[2:] for branch in branches]
+        if "master" in branches:
+            branch_name = "master"
+        elif not branches:
+            branch_name = None
+        else:
+            # Pick the first branch that is returned. Good or bad.
+            branch_name = branches[0]
+
+    pieces["branch"] = branch_name
+
     # parse describe_out. It will be like TAG-NUM-gHEX[-dirty] or HEX[-dirty]
     # TAG might have hyphens.
     git_describe = describe_out
@@ -688,7 +856,7 @@ def git_pieces_from_vcs(tag_prefix, root, verbose, run_command=run_command):
         # TAG-NUM-gHEX
         mo = re.search(r'^(.+)-(\d+)-g([0-9a-f]+)$', git_describe)
         if not mo:
-            # unparseable. Maybe git-describe is misbehaving?
+            # unparsable. Maybe git-describe is misbehaving?
             pieces["error"] = ("unable to parse git-describe output: '%%s'"
                                %% describe_out)
             return pieces
@@ -713,26 +881,27 @@ def git_pieces_from_vcs(tag_prefix, root, verbose, run_command=run_command):
     else:
         # HEX: no tags
         pieces["closest-tag"] = None
-        count_out, rc = run_command(GITS, ["rev-list", "HEAD", "--count"],
-                                    cwd=root)
-        pieces["distance"] = int(count_out)  # total number of commits
+        out, rc = runner(GITS, ["rev-list", "HEAD", "--left-right"], cwd=root)
+        pieces["distance"] = len(out.split())  # total number of commits
 
     # commit date: see ISO-8601 comment in git_versions_from_keywords()
-    date = run_command(GITS, ["show", "-s", "--format=%%ci", "HEAD"],
-                       cwd=root)[0].strip()
+    date = runner(GITS, ["show", "-s", "--format=%%ci", "HEAD"], cwd=root)[0].strip()
+    # Use only the last line.  Previous lines may contain GPG signature
+    # information.
+    date = date.splitlines()[-1]
     pieces["date"] = date.strip().replace(" ", "T", 1).replace(" ", "", 1)
 
     return pieces
 
 
-def plus_or_dot(pieces):
+def plus_or_dot(pieces: Dict[str, Any]) -> str:
     """Return a + if we don't already have one, else return a ."""
     if "+" in pieces.get("closest-tag", ""):
         return "."
     return "+"
 
 
-def render_pep440(pieces):
+def render_pep440(pieces: Dict[str, Any]) -> str:
     """Build up version string, with post-release "local version identifier".
 
     Our goal: TAG[+DISTANCE.gHEX[.dirty]] . Note that if you
@@ -757,23 +926,71 @@ def render_pep440(pieces):
     return rendered
 
 
-def render_pep440_pre(pieces):
-    """TAG[.post.devDISTANCE] -- No -dirty.
+def render_pep440_branch(pieces: Dict[str, Any]) -> str:
+    """TAG[[.dev0]+DISTANCE.gHEX[.dirty]] .
+
+    The ".dev0" means not master branch. Note that .dev0 sorts backwards
+    (a feature branch will appear "older" than the master branch).
 
     Exceptions:
-    1: no tags. 0.post.devDISTANCE
+    1: no tags. 0[.dev0]+untagged.DISTANCE.gHEX[.dirty]
     """
     if pieces["closest-tag"]:
         rendered = pieces["closest-tag"]
+        if pieces["distance"] or pieces["dirty"]:
+            if pieces["branch"] != "master":
+                rendered += ".dev0"
+            rendered += plus_or_dot(pieces)
+            rendered += "%%d.g%%s" %% (pieces["distance"], pieces["short"])
+            if pieces["dirty"]:
+                rendered += ".dirty"
+    else:
+        # exception #1
+        rendered = "0"
+        if pieces["branch"] != "master":
+            rendered += ".dev0"
+        rendered += "+untagged.%%d.g%%s" %% (pieces["distance"],
+                                          pieces["short"])
+        if pieces["dirty"]:
+            rendered += ".dirty"
+    return rendered
+
+
+def pep440_split_post(ver: str) -> Tuple[str, Optional[int]]:
+    """Split pep440 version string at the post-release segment.
+
+    Returns the release segments before the post-release and the
+    post-release version number (or -1 if no post-release segment is present).
+    """
+    vc = str.split(ver, ".post")
+    return vc[0], int(vc[1] or 0) if len(vc) == 2 else None
+
+
+def render_pep440_pre(pieces: Dict[str, Any]) -> str:
+    """TAG[.postN.devDISTANCE] -- No -dirty.
+
+    Exceptions:
+    1: no tags. 0.post0.devDISTANCE
+    """
+    if pieces["closest-tag"]:
         if pieces["distance"]:
-            rendered += ".post.dev%%d" %% pieces["distance"]
+            # update the post release segment
+            tag_version, post_version = pep440_split_post(pieces["closest-tag"])
+            rendered = tag_version
+            if post_version is not None:
+                rendered += ".post%%d.dev%%d" %% (post_version + 1, pieces["distance"])
+            else:
+                rendered += ".post0.dev%%d" %% (pieces["distance"])
+        else:
+            # no commits, use the tag as the version
+            rendered = pieces["closest-tag"]
     else:
         # exception #1
-        rendered = "0.post.dev%%d" %% pieces["distance"]
+        rendered = "0.post0.dev%%d" %% pieces["distance"]
     return rendered
 
 
-def render_pep440_post(pieces):
+def render_pep440_post(pieces: Dict[str, Any]) -> str:
     """TAG[.postDISTANCE[.dev0]+gHEX] .
 
     The ".dev0" means dirty. Note that .dev0 sorts backwards
@@ -800,12 +1017,41 @@ def render_pep440_post(pieces):
     return rendered
 
 
-def render_pep440_old(pieces):
+def render_pep440_post_branch(pieces: Dict[str, Any]) -> str:
+    """TAG[.postDISTANCE[.dev0]+gHEX[.dirty]] .
+
+    The ".dev0" means not master branch.
+
+    Exceptions:
+    1: no tags. 0.postDISTANCE[.dev0]+gHEX[.dirty]
+    """
+    if pieces["closest-tag"]:
+        rendered = pieces["closest-tag"]
+        if pieces["distance"] or pieces["dirty"]:
+            rendered += ".post%%d" %% pieces["distance"]
+            if pieces["branch"] != "master":
+                rendered += ".dev0"
+            rendered += plus_or_dot(pieces)
+            rendered += "g%%s" %% pieces["short"]
+            if pieces["dirty"]:
+                rendered += ".dirty"
+    else:
+        # exception #1
+        rendered = "0.post%%d" %% pieces["distance"]
+        if pieces["branch"] != "master":
+            rendered += ".dev0"
+        rendered += "+g%%s" %% pieces["short"]
+        if pieces["dirty"]:
+            rendered += ".dirty"
+    return rendered
+
+
+def render_pep440_old(pieces: Dict[str, Any]) -> str:
     """TAG[.postDISTANCE[.dev0]] .
 
     The ".dev0" means dirty.
 
-    Eexceptions:
+    Exceptions:
     1: no tags. 0.postDISTANCE[.dev0]
     """
     if pieces["closest-tag"]:
@@ -822,7 +1068,7 @@ def render_pep440_old(pieces):
     return rendered
 
 
-def render_git_describe(pieces):
+def render_git_describe(pieces: Dict[str, Any]) -> str:
     """TAG[-DISTANCE-gHEX][-dirty].
 
     Like 'git describe --tags --dirty --always'.
@@ -842,7 +1088,7 @@ def render_git_describe(pieces):
     return rendered
 
 
-def render_git_describe_long(pieces):
+def render_git_describe_long(pieces: Dict[str, Any]) -> str:
     """TAG-DISTANCE-gHEX[-dirty].
 
     Like 'git describe --tags --dirty --always -long'.
@@ -862,7 +1108,7 @@ def render_git_describe_long(pieces):
     return rendered
 
 
-def render(pieces, style):
+def render(pieces: Dict[str, Any], style: str) -> Dict[str, Any]:
     """Render the given version pieces into the requested style."""
     if pieces["error"]:
         return {"version": "unknown",
@@ -876,10 +1122,14 @@ def render(pieces, style):
 
     if style == "pep440":
         rendered = render_pep440(pieces)
+    elif style == "pep440-branch":
+        rendered = render_pep440_branch(pieces)
     elif style == "pep440-pre":
         rendered = render_pep440_pre(pieces)
     elif style == "pep440-post":
         rendered = render_pep440_post(pieces)
+    elif style == "pep440-post-branch":
+        rendered = render_pep440_post_branch(pieces)
     elif style == "pep440-old":
         rendered = render_pep440_old(pieces)
     elif style == "git-describe":
@@ -894,7 +1144,7 @@ def render(pieces, style):
             "date": pieces.get("date")}
 
 
-def get_versions():
+def get_versions() -> Dict[str, Any]:
     """Get version information or return default if unable to do so."""
     # I am in _version.py, which lives at ROOT/VERSIONFILE_SOURCE. If we have
     # __file__, we can work backwards from there to the root. Some
@@ -915,7 +1165,7 @@ def get_versions():
         # versionfile_source is the relative path from the top of the source
         # tree (where the .git directory might live) to this file. Invert
         # this to find the root from __file__.
-        for i in cfg.versionfile_source.split('/'):
+        for _ in cfg.versionfile_source.split('/'):
             root = os.path.dirname(root)
     except NameError:
         return {"version": "0+unknown", "full-revisionid": None,
@@ -942,41 +1192,48 @@ def get_versions():
 
 
 @register_vcs_handler("git", "get_keywords")
-def git_get_keywords(versionfile_abs):
+def git_get_keywords(versionfile_abs: str) -> Dict[str, str]:
     """Extract version information from the given file."""
     # the code embedded in _version.py can just fetch the value of these
     # keywords. When used from setup.py, we don't want to import _version.py,
     # so we do it with a regexp instead. This function is not used from
     # _version.py.
-    keywords = {}
+    keywords: Dict[str, str] = {}
     try:
-        f = open(versionfile_abs, "r")
-        for line in f.readlines():
-            if line.strip().startswith("git_refnames ="):
-                mo = re.search(r'=\s*"(.*)"', line)
-                if mo:
-                    keywords["refnames"] = mo.group(1)
-            if line.strip().startswith("git_full ="):
-                mo = re.search(r'=\s*"(.*)"', line)
-                if mo:
-                    keywords["full"] = mo.group(1)
-            if line.strip().startswith("git_date ="):
-                mo = re.search(r'=\s*"(.*)"', line)
-                if mo:
-                    keywords["date"] = mo.group(1)
-        f.close()
-    except EnvironmentError:
+        with open(versionfile_abs, "r") as fobj:
+            for line in fobj:
+                if line.strip().startswith("git_refnames ="):
+                    mo = re.search(r'=\s*"(.*)"', line)
+                    if mo:
+                        keywords["refnames"] = mo.group(1)
+                if line.strip().startswith("git_full ="):
+                    mo = re.search(r'=\s*"(.*)"', line)
+                    if mo:
+                        keywords["full"] = mo.group(1)
+                if line.strip().startswith("git_date ="):
+                    mo = re.search(r'=\s*"(.*)"', line)
+                    if mo:
+                        keywords["date"] = mo.group(1)
+    except OSError:
         pass
     return keywords
 
 
 @register_vcs_handler("git", "keywords")
-def git_versions_from_keywords(keywords, tag_prefix, verbose):
+def git_versions_from_keywords(
+    keywords: Dict[str, str],
+    tag_prefix: str,
+    verbose: bool,
+) -> Dict[str, Any]:
     """Get version information from git keywords."""
-    if not keywords:
-        raise NotThisMethod("no keywords at all, weird")
+    if "refnames" not in keywords:
+        raise NotThisMethod("Short version file found")
     date = keywords.get("date")
     if date is not None:
+        # Use only the last line.  Previous lines may contain GPG signature
+        # information.
+        date = date.splitlines()[-1]
+
         # git-2.2.0 added "%cI", which expands to an ISO-8601 -compliant
         # datestamp. However we prefer "%ci" (which expands to an "ISO-8601
         # -like" string, which we must then edit to make compliant), because
@@ -989,11 +1246,11 @@ def git_versions_from_keywords(keywords, tag_prefix, verbose):
         if verbose:
             print("keywords are unexpanded, not using")
         raise NotThisMethod("unexpanded keywords, not a git-archive tarball")
-    refs = set([r.strip() for r in refnames.strip("()").split(",")])
+    refs = {r.strip() for r in refnames.strip("()").split(",")}
     # starting in git-1.8.3, tags are listed as "tag: foo-1.0" instead of
     # just "foo-1.0". If we see a "tag: " prefix, prefer those.
     TAG = "tag: "
-    tags = set([r[len(TAG):] for r in refs if r.startswith(TAG)])
+    tags = {r[len(TAG):] for r in refs if r.startswith(TAG)}
     if not tags:
         # Either we're using git < 1.8.3, or there really are no tags. We use
         # a heuristic: assume all version tags have a digit. The old git %d
@@ -1002,7 +1259,7 @@ def git_versions_from_keywords(keywords, tag_prefix, verbose):
         # between branches and tags. By ignoring refnames without digits, we
         # filter out many common branch names like "release" and
         # "stabilization", as well as "HEAD" and "master".
-        tags = set([r for r in refs if re.search(r'\d', r)])
+        tags = {r for r in refs if re.search(r'\d', r)}
         if verbose:
             print("discarding '%s', no digits" % ",".join(refs - tags))
     if verbose:
@@ -1011,6 +1268,11 @@ def git_versions_from_keywords(keywords, tag_prefix, verbose):
         # sorting will prefer e.g. "2.0" over "2.0rc1"
         if ref.startswith(tag_prefix):
             r = ref[len(tag_prefix):]
+            # Filter out refs that exactly match prefix or that don't start
+            # with a number once the prefix is stripped (mostly a concern
+            # when prefix is '')
+            if not re.match(r'\d', r):
+                continue
             if verbose:
                 print("picking %s" % r)
             return {"version": r,
@@ -1026,7 +1288,12 @@ def git_versions_from_keywords(keywords, tag_prefix, verbose):
 
 
 @register_vcs_handler("git", "pieces_from_vcs")
-def git_pieces_from_vcs(tag_prefix, root, verbose, run_command=run_command):
+def git_pieces_from_vcs(
+    tag_prefix: str,
+    root: str,
+    verbose: bool,
+    runner: Callable = run_command
+) -> Dict[str, Any]:
     """Get version from 'git describe' in the root of the source tree.
 
     This only gets called if the git-archive 'subst' keywords were *not*
@@ -1037,8 +1304,15 @@ def git_pieces_from_vcs(tag_prefix, root, verbose, run_command=run_command):
     if sys.platform == "win32":
         GITS = ["git.cmd", "git.exe"]
 
-    out, rc = run_command(GITS, ["rev-parse", "--git-dir"], cwd=root,
-                          hide_stderr=True)
+    # GIT_DIR can interfere with correct operation of Versioneer.
+    # It may be intended to be passed to the Versioneer-versioned project,
+    # but that should not change where we get our version from.
+    env = os.environ.copy()
+    env.pop("GIT_DIR", None)
+    runner = functools.partial(runner, env=env)
+
+    _, rc = runner(GITS, ["rev-parse", "--git-dir"], cwd=root,
+                   hide_stderr=not verbose)
     if rc != 0:
         if verbose:
             print("Directory %s not under git control" % root)
@@ -1046,24 +1320,57 @@ def git_pieces_from_vcs(tag_prefix, root, verbose, run_command=run_command):
 
     # if there is a tag matching tag_prefix, this yields TAG-NUM-gHEX[-dirty]
     # if there isn't one, this yields HEX[-dirty] (no NUM)
-    describe_out, rc = run_command(GITS, ["describe", "--tags", "--dirty",
-                                          "--always", "--long",
-                                          "--match", "%s*" % tag_prefix],
-                                   cwd=root)
+    describe_out, rc = runner(GITS, [
+        "describe", "--tags", "--dirty", "--always", "--long",
+        "--match", f"{tag_prefix}[[:digit:]]*"
+    ], cwd=root)
     # --long was added in git-1.5.5
     if describe_out is None:
         raise NotThisMethod("'git describe' failed")
     describe_out = describe_out.strip()
-    full_out, rc = run_command(GITS, ["rev-parse", "HEAD"], cwd=root)
+    full_out, rc = runner(GITS, ["rev-parse", "HEAD"], cwd=root)
     if full_out is None:
         raise NotThisMethod("'git rev-parse' failed")
     full_out = full_out.strip()
 
-    pieces = {}
+    pieces: Dict[str, Any] = {}
     pieces["long"] = full_out
     pieces["short"] = full_out[:7]  # maybe improved later
     pieces["error"] = None
 
+    branch_name, rc = runner(GITS, ["rev-parse", "--abbrev-ref", "HEAD"],
+                             cwd=root)
+    # --abbrev-ref was added in git-1.6.3
+    if rc != 0 or branch_name is None:
+        raise NotThisMethod("'git rev-parse --abbrev-ref' returned error")
+    branch_name = branch_name.strip()
+
+    if branch_name == "HEAD":
+        # If we aren't exactly on a branch, pick a branch which represents
+        # the current commit. If all else fails, we are on a branchless
+        # commit.
+        branches, rc = runner(GITS, ["branch", "--contains"], cwd=root)
+        # --contains was added in git-1.5.4
+        if rc != 0 or branches is None:
+            raise NotThisMethod("'git branch --contains' returned error")
+        branches = branches.split("\n")
+
+        # Remove the first line if we're running detached
+        if "(" in branches[0]:
+            branches.pop(0)
+
+        # Strip off the leading "* " from the list of branches.
+        branches = [branch[2:] for branch in branches]
+        if "master" in branches:
+            branch_name = "master"
+        elif not branches:
+            branch_name = None
+        else:
+            # Pick the first branch that is returned. Good or bad.
+            branch_name = branches[0]
+
+    pieces["branch"] = branch_name
+
     # parse describe_out. It will be like TAG-NUM-gHEX[-dirty] or HEX[-dirty]
     # TAG might have hyphens.
     git_describe = describe_out
@@ -1080,7 +1387,7 @@ def git_pieces_from_vcs(tag_prefix, root, verbose, run_command=run_command):
         # TAG-NUM-gHEX
         mo = re.search(r'^(.+)-(\d+)-g([0-9a-f]+)$', git_describe)
         if not mo:
-            # unparseable. Maybe git-describe is misbehaving?
+            # unparsable. Maybe git-describe is misbehaving?
             pieces["error"] = ("unable to parse git-describe output: '%s'"
                                % describe_out)
             return pieces
@@ -1105,19 +1412,20 @@ def git_pieces_from_vcs(tag_prefix, root, verbose, run_command=run_command):
     else:
         # HEX: no tags
         pieces["closest-tag"] = None
-        count_out, rc = run_command(GITS, ["rev-list", "HEAD", "--count"],
-                                    cwd=root)
-        pieces["distance"] = int(count_out)  # total number of commits
+        out, rc = runner(GITS, ["rev-list", "HEAD", "--left-right"], cwd=root)
+        pieces["distance"] = len(out.split())  # total number of commits
 
     # commit date: see ISO-8601 comment in git_versions_from_keywords()
-    date = run_command(GITS, ["show", "-s", "--format=%ci", "HEAD"],
-                       cwd=root)[0].strip()
+    date = runner(GITS, ["show", "-s", "--format=%ci", "HEAD"], cwd=root)[0].strip()
+    # Use only the last line.  Previous lines may contain GPG signature
+    # information.
+    date = date.splitlines()[-1]
     pieces["date"] = date.strip().replace(" ", "T", 1).replace(" ", "", 1)
 
     return pieces
 
 
-def do_vcs_install(manifest_in, versionfile_source, ipy):
+def do_vcs_install(versionfile_source: str, ipy: Optional[str]) -> None:
     """Git-specific installation logic for Versioneer.
 
     For Git, this means creating/changing .gitattributes to mark _version.py
@@ -1126,36 +1434,40 @@ def do_vcs_install(manifest_in, versionfile_source, ipy):
     GITS = ["git"]
     if sys.platform == "win32":
         GITS = ["git.cmd", "git.exe"]
-    files = [manifest_in, versionfile_source]
+    files = [versionfile_source]
     if ipy:
         files.append(ipy)
-    try:
-        me = __file__
-        if me.endswith(".pyc") or me.endswith(".pyo"):
-            me = os.path.splitext(me)[0] + ".py"
-        versioneer_file = os.path.relpath(me)
-    except NameError:
-        versioneer_file = "versioneer.py"
-    files.append(versioneer_file)
+    if "VERSIONEER_PEP518" not in globals():
+        try:
+            my_path = __file__
+            if my_path.endswith((".pyc", ".pyo")):
+                my_path = os.path.splitext(my_path)[0] + ".py"
+            versioneer_file = os.path.relpath(my_path)
+        except NameError:
+            versioneer_file = "versioneer.py"
+        files.append(versioneer_file)
     present = False
     try:
-        f = open(".gitattributes", "r")
-        for line in f.readlines():
-            if line.strip().startswith(versionfile_source):
-                if "export-subst" in line.strip().split()[1:]:
-                    present = True
-        f.close()
-    except EnvironmentError:
+        with open(".gitattributes", "r") as fobj:
+            for line in fobj:
+                if line.strip().startswith(versionfile_source):
+                    if "export-subst" in line.strip().split()[1:]:
+                        present = True
+                        break
+    except OSError:
         pass
     if not present:
-        f = open(".gitattributes", "a+")
-        f.write("%s export-subst\n" % versionfile_source)
-        f.close()
+        with open(".gitattributes", "a+") as fobj:
+            fobj.write(f"{versionfile_source} export-subst\n")
         files.append(".gitattributes")
     run_command(GITS, ["add", "--"] + files)
 
 
-def versions_from_parentdir(parentdir_prefix, root, verbose):
+def versions_from_parentdir(
+    parentdir_prefix: str,
+    root: str,
+    verbose: bool,
+) -> Dict[str, Any]:
     """Try to determine the version from the parent directory name.
 
     Source tarballs conventionally unpack into a directory that includes both
@@ -1164,15 +1476,14 @@ def versions_from_parentdir(parentdir_prefix, root, verbose):
     """
     rootdirs = []
 
-    for i in range(3):
+    for _ in range(3):
         dirname = os.path.basename(root)
         if dirname.startswith(parentdir_prefix):
             return {"version": dirname[len(parentdir_prefix):],
                     "full-revisionid": None,
                     "dirty": False, "error": None, "date": None}
-        else:
-            rootdirs.append(root)
-            root = os.path.dirname(root)  # up a level
+        rootdirs.append(root)
+        root = os.path.dirname(root)  # up a level
 
     if verbose:
         print("Tried directories %s but none started with prefix %s" %
@@ -1181,7 +1492,7 @@ def versions_from_parentdir(parentdir_prefix, root, verbose):
 
 
 SHORT_VERSION_PY = """
-# This file was generated by 'versioneer.py' (0.18) from
+# This file was generated by 'versioneer.py' (0.29) from
 # revision-control system data, or from the parent directory name of an
 # unpacked source archive. Distribution tarballs contain a pre-generated copy
 # of this file.
@@ -1198,12 +1509,12 @@ def get_versions():
 """
 
 
-def versions_from_file(filename):
+def versions_from_file(filename: str) -> Dict[str, Any]:
     """Try to determine the version from _version.py if present."""
     try:
         with open(filename) as f:
             contents = f.read()
-    except EnvironmentError:
+    except OSError:
         raise NotThisMethod("unable to read _version.py")
     mo = re.search(r"version_json = '''\n(.*)'''  # END VERSION_JSON",
                    contents, re.M | re.S)
@@ -1215,9 +1526,8 @@ def versions_from_file(filename):
     return json.loads(mo.group(1))
 
 
-def write_to_version_file(filename, versions):
+def write_to_version_file(filename: str, versions: Dict[str, Any]) -> None:
     """Write the given version number to the given _version.py file."""
-    os.unlink(filename)
     contents = json.dumps(versions, sort_keys=True,
                           indent=1, separators=(",", ": "))
     with open(filename, "w") as f:
@@ -1226,14 +1536,14 @@ def write_to_version_file(filename, versions):
     print("set %s to '%s'" % (filename, versions["version"]))
 
 
-def plus_or_dot(pieces):
+def plus_or_dot(pieces: Dict[str, Any]) -> str:
     """Return a + if we don't already have one, else return a ."""
     if "+" in pieces.get("closest-tag", ""):
         return "."
     return "+"
 
 
-def render_pep440(pieces):
+def render_pep440(pieces: Dict[str, Any]) -> str:
     """Build up version string, with post-release "local version identifier".
 
     Our goal: TAG[+DISTANCE.gHEX[.dirty]] . Note that if you
@@ -1258,23 +1568,71 @@ def render_pep440(pieces):
     return rendered
 
 
-def render_pep440_pre(pieces):
-    """TAG[.post.devDISTANCE] -- No -dirty.
+def render_pep440_branch(pieces: Dict[str, Any]) -> str:
+    """TAG[[.dev0]+DISTANCE.gHEX[.dirty]] .
+
+    The ".dev0" means not master branch. Note that .dev0 sorts backwards
+    (a feature branch will appear "older" than the master branch).
 
     Exceptions:
-    1: no tags. 0.post.devDISTANCE
+    1: no tags. 0[.dev0]+untagged.DISTANCE.gHEX[.dirty]
     """
     if pieces["closest-tag"]:
         rendered = pieces["closest-tag"]
+        if pieces["distance"] or pieces["dirty"]:
+            if pieces["branch"] != "master":
+                rendered += ".dev0"
+            rendered += plus_or_dot(pieces)
+            rendered += "%d.g%s" % (pieces["distance"], pieces["short"])
+            if pieces["dirty"]:
+                rendered += ".dirty"
+    else:
+        # exception #1
+        rendered = "0"
+        if pieces["branch"] != "master":
+            rendered += ".dev0"
+        rendered += "+untagged.%d.g%s" % (pieces["distance"],
+                                          pieces["short"])
+        if pieces["dirty"]:
+            rendered += ".dirty"
+    return rendered
+
+
+def pep440_split_post(ver: str) -> Tuple[str, Optional[int]]:
+    """Split pep440 version string at the post-release segment.
+
+    Returns the release segments before the post-release and the
+    post-release version number (or -1 if no post-release segment is present).
+    """
+    vc = str.split(ver, ".post")
+    return vc[0], int(vc[1] or 0) if len(vc) == 2 else None
+
+
+def render_pep440_pre(pieces: Dict[str, Any]) -> str:
+    """TAG[.postN.devDISTANCE] -- No -dirty.
+
+    Exceptions:
+    1: no tags. 0.post0.devDISTANCE
+    """
+    if pieces["closest-tag"]:
         if pieces["distance"]:
-            rendered += ".post.dev%d" % pieces["distance"]
+            # update the post release segment
+            tag_version, post_version = pep440_split_post(pieces["closest-tag"])
+            rendered = tag_version
+            if post_version is not None:
+                rendered += ".post%d.dev%d" % (post_version + 1, pieces["distance"])
+            else:
+                rendered += ".post0.dev%d" % (pieces["distance"])
+        else:
+            # no commits, use the tag as the version
+            rendered = pieces["closest-tag"]
     else:
         # exception #1
-        rendered = "0.post.dev%d" % pieces["distance"]
+        rendered = "0.post0.dev%d" % pieces["distance"]
     return rendered
 
 
-def render_pep440_post(pieces):
+def render_pep440_post(pieces: Dict[str, Any]) -> str:
     """TAG[.postDISTANCE[.dev0]+gHEX] .
 
     The ".dev0" means dirty. Note that .dev0 sorts backwards
@@ -1301,12 +1659,41 @@ def render_pep440_post(pieces):
     return rendered
 
 
-def render_pep440_old(pieces):
+def render_pep440_post_branch(pieces: Dict[str, Any]) -> str:
+    """TAG[.postDISTANCE[.dev0]+gHEX[.dirty]] .
+
+    The ".dev0" means not master branch.
+
+    Exceptions:
+    1: no tags. 0.postDISTANCE[.dev0]+gHEX[.dirty]
+    """
+    if pieces["closest-tag"]:
+        rendered = pieces["closest-tag"]
+        if pieces["distance"] or pieces["dirty"]:
+            rendered += ".post%d" % pieces["distance"]
+            if pieces["branch"] != "master":
+                rendered += ".dev0"
+            rendered += plus_or_dot(pieces)
+            rendered += "g%s" % pieces["short"]
+            if pieces["dirty"]:
+                rendered += ".dirty"
+    else:
+        # exception #1
+        rendered = "0.post%d" % pieces["distance"]
+        if pieces["branch"] != "master":
+            rendered += ".dev0"
+        rendered += "+g%s" % pieces["short"]
+        if pieces["dirty"]:
+            rendered += ".dirty"
+    return rendered
+
+
+def render_pep440_old(pieces: Dict[str, Any]) -> str:
     """TAG[.postDISTANCE[.dev0]] .
 
     The ".dev0" means dirty.
 
-    Eexceptions:
+    Exceptions:
     1: no tags. 0.postDISTANCE[.dev0]
     """
     if pieces["closest-tag"]:
@@ -1323,7 +1710,7 @@ def render_pep440_old(pieces):
     return rendered
 
 
-def render_git_describe(pieces):
+def render_git_describe(pieces: Dict[str, Any]) -> str:
     """TAG[-DISTANCE-gHEX][-dirty].
 
     Like 'git describe --tags --dirty --always'.
@@ -1343,7 +1730,7 @@ def render_git_describe(pieces):
     return rendered
 
 
-def render_git_describe_long(pieces):
+def render_git_describe_long(pieces: Dict[str, Any]) -> str:
     """TAG-DISTANCE-gHEX[-dirty].
 
     Like 'git describe --tags --dirty --always -long'.
@@ -1363,7 +1750,7 @@ def render_git_describe_long(pieces):
     return rendered
 
 
-def render(pieces, style):
+def render(pieces: Dict[str, Any], style: str) -> Dict[str, Any]:
     """Render the given version pieces into the requested style."""
     if pieces["error"]:
         return {"version": "unknown",
@@ -1377,10 +1764,14 @@ def render(pieces, style):
 
     if style == "pep440":
         rendered = render_pep440(pieces)
+    elif style == "pep440-branch":
+        rendered = render_pep440_branch(pieces)
     elif style == "pep440-pre":
         rendered = render_pep440_pre(pieces)
     elif style == "pep440-post":
         rendered = render_pep440_post(pieces)
+    elif style == "pep440-post-branch":
+        rendered = render_pep440_post_branch(pieces)
     elif style == "pep440-old":
         rendered = render_pep440_old(pieces)
     elif style == "git-describe":
@@ -1399,7 +1790,7 @@ class VersioneerBadRootError(Exception):
     """The project root directory is unknown or missing key files."""
 
 
-def get_versions(verbose=False):
+def get_versions(verbose: bool = False) -> Dict[str, Any]:
     """Get the project version from whatever source is available.
 
     Returns dict with two keys: 'version' and 'full'.
@@ -1414,7 +1805,7 @@ def get_versions(verbose=False):
     assert cfg.VCS is not None, "please set [versioneer]VCS= in setup.cfg"
     handlers = HANDLERS.get(cfg.VCS)
     assert handlers, "unrecognized VCS '%s'" % cfg.VCS
-    verbose = verbose or cfg.verbose
+    verbose = verbose or bool(cfg.verbose)  # `bool()` used to avoid `None`
     assert cfg.versionfile_source is not None, \
         "please set versioneer.versionfile_source"
     assert cfg.tag_prefix is not None, "please set versioneer.tag_prefix"
@@ -1475,13 +1866,17 @@ def get_versions(verbose=False):
             "date": None}
 
 
-def get_version():
+def get_version() -> str:
     """Get the short version string for this project."""
     return get_versions()["version"]
 
 
-def get_cmdclass():
-    """Get the custom setuptools/distutils subclasses used by Versioneer."""
+def get_cmdclass(cmdclass: Optional[Dict[str, Any]] = None):
+    """Get the custom setuptools subclasses used by Versioneer.
+
+    If the package uses a different cmdclass (e.g. one from numpy), it
+    should be provide as an argument.
+    """
     if "versioneer" in sys.modules:
         del sys.modules["versioneer"]
         # this fixes the "python setup.py develop" case (also 'install' and
@@ -1495,25 +1890,25 @@ def get_cmdclass():
         # parent is protected against the child's "import versioneer". By
         # removing ourselves from sys.modules here, before the child build
         # happens, we protect the child from the parent's versioneer too.
-        # Also see https://github.com/warner/python-versioneer/issues/52
+        # Also see https://github.com/python-versioneer/python-versioneer/issues/52
 
-    cmds = {}
+    cmds = {} if cmdclass is None else cmdclass.copy()
 
-    # we add "version" to both distutils and setuptools
-    from distutils.core import Command
+    # we add "version" to setuptools
+    from setuptools import Command
 
     class cmd_version(Command):
         description = "report generated version string"
-        user_options = []
-        boolean_options = []
+        user_options: List[Tuple[str, str, str]] = []
+        boolean_options: List[str] = []
 
-        def initialize_options(self):
+        def initialize_options(self) -> None:
             pass
 
-        def finalize_options(self):
+        def finalize_options(self) -> None:
             pass
 
-        def run(self):
+        def run(self) -> None:
             vers = get_versions(verbose=True)
             print("Version: %s" % vers["version"])
             print(" full-revisionid: %s" % vers.get("full-revisionid"))
@@ -1523,7 +1918,7 @@ def run(self):
                 print(" error: %s" % vers["error"])
     cmds["version"] = cmd_version
 
-    # we override "build_py" in both distutils and setuptools
+    # we override "build_py" in setuptools
     #
     # most invocation pathways end up running build_py:
     #  distutils/build -> build_py
@@ -1538,18 +1933,25 @@ def run(self):
     #   then does setup.py bdist_wheel, or sometimes setup.py install
     #  setup.py egg_info -> ?
 
+    # pip install -e . and setuptool/editable_wheel will invoke build_py
+    # but the build_py command is not expected to copy any files.
+
     # we override different "build_py" commands for both environments
-    if "setuptools" in sys.modules:
-        from setuptools.command.build_py import build_py as _build_py
+    if 'build_py' in cmds:
+        _build_py: Any = cmds['build_py']
     else:
-        from distutils.command.build_py import build_py as _build_py
+        from setuptools.command.build_py import build_py as _build_py
 
     class cmd_build_py(_build_py):
-        def run(self):
+        def run(self) -> None:
             root = get_root()
             cfg = get_config_from_root(root)
             versions = get_versions()
             _build_py.run(self)
+            if getattr(self, "editable_mode", False):
+                # During editable installs `.py` and data files are
+                # not copied to build_lib
+                return
             # now locate _version.py in the new build/ directory and replace
             # it with an updated value
             if cfg.versionfile_build:
@@ -1559,8 +1961,40 @@ def run(self):
                 write_to_version_file(target_versionfile, versions)
     cmds["build_py"] = cmd_build_py
 
+    if 'build_ext' in cmds:
+        _build_ext: Any = cmds['build_ext']
+    else:
+        from setuptools.command.build_ext import build_ext as _build_ext
+
+    class cmd_build_ext(_build_ext):
+        def run(self) -> None:
+            root = get_root()
+            cfg = get_config_from_root(root)
+            versions = get_versions()
+            _build_ext.run(self)
+            if self.inplace:
+                # build_ext --inplace will only build extensions in
+                # build/lib<..> dir with no _version.py to write to.
+                # As in place builds will already have a _version.py
+                # in the module dir, we do not need to write one.
+                return
+            # now locate _version.py in the new build/ directory and replace
+            # it with an updated value
+            if not cfg.versionfile_build:
+                return
+            target_versionfile = os.path.join(self.build_lib,
+                                              cfg.versionfile_build)
+            if not os.path.exists(target_versionfile):
+                print(f"Warning: {target_versionfile} does not exist, skipping "
+                      "version update. This can happen if you are running build_ext "
+                      "without first running build_py.")
+                return
+            print("UPDATING %s" % target_versionfile)
+            write_to_version_file(target_versionfile, versions)
+    cmds["build_ext"] = cmd_build_ext
+
     if "cx_Freeze" in sys.modules:  # cx_freeze enabled?
-        from cx_Freeze.dist import build_exe as _build_exe
+        from cx_Freeze.dist import build_exe as _build_exe  # type: ignore
         # nczeczulin reports that py2exe won't like the pep440-style string
         # as FILEVERSION, but it can be used for PRODUCTVERSION, e.g.
         # setup(console=[{
@@ -1569,7 +2003,7 @@ def run(self):
         #   ...
 
         class cmd_build_exe(_build_exe):
-            def run(self):
+            def run(self) -> None:
                 root = get_root()
                 cfg = get_config_from_root(root)
                 versions = get_versions()
@@ -1593,12 +2027,12 @@ def run(self):
 
     if 'py2exe' in sys.modules:  # py2exe enabled?
         try:
-            from py2exe.distutils_buildexe import py2exe as _py2exe  # py3
+            from py2exe.setuptools_buildexe import py2exe as _py2exe  # type: ignore
         except ImportError:
-            from py2exe.build_exe import py2exe as _py2exe  # py2
+            from py2exe.distutils_buildexe import py2exe as _py2exe  # type: ignore
 
         class cmd_py2exe(_py2exe):
-            def run(self):
+            def run(self) -> None:
                 root = get_root()
                 cfg = get_config_from_root(root)
                 versions = get_versions()
@@ -1619,14 +2053,51 @@ def run(self):
                              })
         cmds["py2exe"] = cmd_py2exe
 
+    # sdist farms its file list building out to egg_info
+    if 'egg_info' in cmds:
+        _egg_info: Any = cmds['egg_info']
+    else:
+        from setuptools.command.egg_info import egg_info as _egg_info
+
+    class cmd_egg_info(_egg_info):
+        def find_sources(self) -> None:
+            # egg_info.find_sources builds the manifest list and writes it
+            # in one shot
+            super().find_sources()
+
+            # Modify the filelist and normalize it
+            root = get_root()
+            cfg = get_config_from_root(root)
+            self.filelist.append('versioneer.py')
+            if cfg.versionfile_source:
+                # There are rare cases where versionfile_source might not be
+                # included by default, so we must be explicit
+                self.filelist.append(cfg.versionfile_source)
+            self.filelist.sort()
+            self.filelist.remove_duplicates()
+
+            # The write method is hidden in the manifest_maker instance that
+            # generated the filelist and was thrown away
+            # We will instead replicate their final normalization (to unicode,
+            # and POSIX-style paths)
+            from setuptools import unicode_utils
+            normalized = [unicode_utils.filesys_decode(f).replace(os.sep, '/')
+                          for f in self.filelist.files]
+
+            manifest_filename = os.path.join(self.egg_info, 'SOURCES.txt')
+            with open(manifest_filename, 'w') as fobj:
+                fobj.write('\n'.join(normalized))
+
+    cmds['egg_info'] = cmd_egg_info
+
     # we override different "sdist" commands for both environments
-    if "setuptools" in sys.modules:
-        from setuptools.command.sdist import sdist as _sdist
+    if 'sdist' in cmds:
+        _sdist: Any = cmds['sdist']
     else:
-        from distutils.command.sdist import sdist as _sdist
+        from setuptools.command.sdist import sdist as _sdist
 
     class cmd_sdist(_sdist):
-        def run(self):
+        def run(self) -> None:
             versions = get_versions()
             self._versioneer_generated_versions = versions
             # unless we update this, the command will keep using the old
@@ -1634,7 +2105,7 @@ def run(self):
             self.distribution.metadata.version = versions["version"]
             return _sdist.run(self)
 
-        def make_release_tree(self, base_dir, files):
+        def make_release_tree(self, base_dir: str, files: List[str]) -> None:
             root = get_root()
             cfg = get_config_from_root(root)
             _sdist.make_release_tree(self, base_dir, files)
@@ -1687,21 +2158,26 @@ def make_release_tree(self, base_dir, files):
 
 """
 
-INIT_PY_SNIPPET = """
+OLD_SNIPPET = """
 from ._version import get_versions
 __version__ = get_versions()['version']
 del get_versions
 """
 
+INIT_PY_SNIPPET = """
+from . import {0}
+__version__ = {0}.get_versions()['version']
+"""
 
-def do_setup():
-    """Main VCS-independent setup function for installing Versioneer."""
+
+def do_setup() -> int:
+    """Do main VCS-independent setup function for installing Versioneer."""
     root = get_root()
     try:
         cfg = get_config_from_root(root)
-    except (EnvironmentError, configparser.NoSectionError,
+    except (OSError, configparser.NoSectionError,
             configparser.NoOptionError) as e:
-        if isinstance(e, (EnvironmentError, configparser.NoSectionError)):
+        if isinstance(e, (OSError, configparser.NoSectionError)):
             print("Adding sample versioneer config to setup.cfg",
                   file=sys.stderr)
             with open(os.path.join(root, "setup.cfg"), "a") as f:
@@ -1721,62 +2197,37 @@ def do_setup():
 
     ipy = os.path.join(os.path.dirname(cfg.versionfile_source),
                        "__init__.py")
+    maybe_ipy: Optional[str] = ipy
     if os.path.exists(ipy):
         try:
             with open(ipy, "r") as f:
                 old = f.read()
-        except EnvironmentError:
+        except OSError:
             old = ""
-        if INIT_PY_SNIPPET not in old:
+        module = os.path.splitext(os.path.basename(cfg.versionfile_source))[0]
+        snippet = INIT_PY_SNIPPET.format(module)
+        if OLD_SNIPPET in old:
+            print(" replacing boilerplate in %s" % ipy)
+            with open(ipy, "w") as f:
+                f.write(old.replace(OLD_SNIPPET, snippet))
+        elif snippet not in old:
             print(" appending to %s" % ipy)
             with open(ipy, "a") as f:
-                f.write(INIT_PY_SNIPPET)
+                f.write(snippet)
         else:
             print(" %s unmodified" % ipy)
     else:
         print(" %s doesn't exist, ok" % ipy)
-        ipy = None
-
-    # Make sure both the top-level "versioneer.py" and versionfile_source
-    # (PKG/_version.py, used by runtime code) are in MANIFEST.in, so
-    # they'll be copied into source distributions. Pip won't be able to
-    # install the package without this.
-    manifest_in = os.path.join(root, "MANIFEST.in")
-    simple_includes = set()
-    try:
-        with open(manifest_in, "r") as f:
-            for line in f:
-                if line.startswith("include "):
-                    for include in line.split()[1:]:
-                        simple_includes.add(include)
-    except EnvironmentError:
-        pass
-    # That doesn't cover everything MANIFEST.in can do
-    # (http://docs.python.org/2/distutils/sourcedist.html#commands), so
-    # it might give some false negatives. Appending redundant 'include'
-    # lines is safe, though.
-    if "versioneer.py" not in simple_includes:
-        print(" appending 'versioneer.py' to MANIFEST.in")
-        with open(manifest_in, "a") as f:
-            f.write("include versioneer.py\n")
-    else:
-        print(" 'versioneer.py' already in MANIFEST.in")
-    if cfg.versionfile_source not in simple_includes:
-        print(" appending versionfile_source ('%s') to MANIFEST.in" %
-              cfg.versionfile_source)
-        with open(manifest_in, "a") as f:
-            f.write("include %s\n" % cfg.versionfile_source)
-    else:
-        print(" versionfile_source already in MANIFEST.in")
+        maybe_ipy = None
 
     # Make VCS-specific changes. For git, this means creating/changing
     # .gitattributes to mark _version.py for export-subst keyword
     # substitution.
-    do_vcs_install(manifest_in, cfg.versionfile_source, ipy)
+    do_vcs_install(cfg.versionfile_source, maybe_ipy)
     return 0
 
 
-def scan_setup_py():
+def scan_setup_py() -> int:
     """Validate the contents of setup.py against Versioneer's expectations."""
     found = set()
     setters = False
@@ -1813,10 +2264,14 @@ def scan_setup_py():
     return errors
 
 
+def setup_command() -> NoReturn:
+    """Set up Versioneer and exit with appropriate error code."""
+    errors = do_setup()
+    errors += scan_setup_py()
+    sys.exit(1 if errors else 0)
+
+
 if __name__ == "__main__":
     cmd = sys.argv[1]
     if cmd == "setup":
-        errors = do_setup()
-        errors += scan_setup_py()
-        if errors:
-            sys.exit(1)
+        setup_command()


Test patch:
diff --git a/Tests/conftest.py b/Tests/conftest.py
new file mode 100644
index 0000000..98e64ae
--- /dev/null
+++ b/Tests/conftest.py
@@ -0,0 +1,5 @@
+def pytest_addoption(parser):
+    parser.addoption("--skip-platform", action="store", default=None, help="Skip tests for a specific platform")
+
+def pytest_configure(config):
+    config.addinivalue_line("markers", "skip_platform: mark test to be skipped for a specific platform")
\ No newline at end of file
diff --git a/Tests/test_IO.py b/Tests/test_IO.py
new file mode 100644
index 0000000..15da72f
--- /dev/null
+++ b/Tests/test_IO.py
@@ -0,0 +1,106 @@
+from open3SPN2 import *
+
+# Unit testing
+def test_DNA_from_pdb():
+    """ Test correct DNA initialization from PDB"""
+    mol = DNA.fromPDB("Tests/1svc/1svc.pdb", template_from_X3DNA=False)
+
+
+def test_DNA_from_gro():
+    """ Test correct DNA initialization from gromacs files"""
+    pass
+
+
+def test_DNA_from_seq():
+    """ Test correct DNA initialization from sequence files"""
+    #return True #Needs X3DNA
+    seq = 'ATACAAAGGTGCGAGGTTTCTATGCTCCCACG'
+    dna = DNA.fromSequence(seq, dna_type='B_curved')
+
+    # Compute the topology for the DNA structure.
+    # Since the dna was generated from the sequence using X3DNA,
+    # it is not necesary to recompute the geometry.
+
+    dna.computeTopology(template_from_X3DNA=True)
+
+    # Create the system.
+    # To set periodic boundary conditions (periodicBox=[50,50,50]).
+    # The periodic box size is in nanometers.
+    dna.periodic = False
+    s = System(dna, periodicBox=None)
+
+    # Add 3SPN2 forces
+    s.add3SPN2forces(verbose=True)
+
+    # Initialize Molecular Dynamics simulations
+    s.initializeMD(temperature=300 * openmm.unit.kelvin, platform_name='CPU')
+    simulation = s.simulation
+
+    # Set initial positions
+    simulation.context.setPositions(s.coord.getPositions())
+
+    energy_unit = openmm.unit.kilojoule_per_mole
+    # Total energy
+    state = simulation.context.getState(getEnergy=True)
+    energy = state.getPotentialEnergy().value_in_unit(energy_unit)
+    print('TotalEnergy', round(energy, 6), energy_unit.get_symbol())
+
+    # Detailed energy
+    energies = {}
+    for force_name, force in s.forces.items():
+        group = force.getForceGroup()
+        state = simulation.context.getState(getEnergy=True, groups=2 ** group)
+        energies[force_name] = state.getPotentialEnergy().value_in_unit(energy_unit)
+
+    for force_name in s.forces.keys():
+        print(force_name, round(energies[force_name], 6), energy_unit.get_symbol())
+
+
+def test_DNA_from_xyz():
+    """Tests the correct parsing from an xyz file"""
+    mol = DNA.fromXYZ('Tests/adna/in00_conf.xyz', template_from_X3DNA=False)
+    assert mol.atoms.at[8, 'name'] == 'P'
+    assert round(mol.atoms.at[188, 'y'], 6) == -8.779343
+
+
+# Functional testing
+
+def test_parse_xyz():
+    """Tests the example trajectory parsing"""
+    xyz_data = parse_xyz('Tests/adna/traj.xyz')
+    assert xyz_data.at[1, 'name'] == 7
+    assert xyz_data.at[1, 'x'] == 4.34621
+
+
+def test_parse_log():
+    """Tests the example log parsing"""
+    log_data = parse_log('Tests/adna/sim.log')
+    assert log_data.at[1, 'Step'] == 2000
+    assert log_data.at[1, 'eexcl'] == 0.45734636
+
+
+def test_transform_zero():
+    transform = Transform()
+    
+    np.testing.assert_almost_equal(transform.full_rotation, np.eye(3), decimal=7)
+    np.testing.assert_almost_equal(transform.half_rotation, np.eye(3), decimal=7)
+    np.testing.assert_almost_equal(transform.full_displacement, [0, 0, 0], decimal=7)
+    np.testing.assert_almost_equal(transform.half_displacement, [0, 0, 0], decimal=7)
+
+def test_transform_known():
+    transform = Transform(0.07, -0.19, 0.07, 
+                          1.8/180*np.pi, -15.0/180*np.pi, 1.5/180*np.pi)
+    
+    expected_full_rotation = np.array([[0.965592, -0.029813, -0.258348], 
+                                       [0.021636,  0.999173, -0.034438], 
+                                       [0.259161,  0.027663,  0.965438]])
+    expected_half_rotation = np.array([[0.991374, -0.014114, -0.130305],
+                                       [0.011951,  0.999778, -0.017370],
+                                       [0.130521,  0.015662,  0.991322]])
+    expected_full_displacement = [0.062957, -0.190337, 0.075553]
+    expected_half_displacement = [0.031478, -0.095169, 0.037777]
+    
+    np.testing.assert_almost_equal(transform.full_rotation, expected_full_rotation, decimal=6)
+    np.testing.assert_almost_equal(transform.half_rotation, expected_half_rotation, decimal=6)
+    np.testing.assert_almost_equal(transform.full_displacement, expected_full_displacement, decimal=6)
+    np.testing.assert_almost_equal(transform.half_displacement, expected_half_displacement, decimal=6)
\ No newline at end of file
diff --git a/Tests/test_energies.py b/Tests/test_energies.py
new file mode 100644
index 0000000..a141ea5
--- /dev/null
+++ b/Tests/test_energies.py
@@ -0,0 +1,135 @@
+from open3SPN2 import DNA, System, parse_log, forces
+import pandas
+import numpy as np
+import openmm.app
+import openmm
+import openmm.unit as unit
+import pytest
+import pandas as pd
+from pathlib import Path
+
+# Define a fixture to load the test sets
+# Define a module-level variable to store the test sets
+_test_sets = pd.read_csv('Tests/test_cases.csv', comment='#')
+_ef = 1 * unit.kilocalorie / unit.kilojoule  # energy scaling factor
+
+@pytest.fixture
+def skip_platform(request):
+    return request.config.getoption("--skip-platform")
+
+id_list = [f"{row['Name']}_{row['Platform']}" for _, row in _test_sets.iterrows()]
+class TestEnergies:
+    """Tests that the energies are the same as the example outputs from lammps"""
+    @pytest.mark.parametrize(
+        "name, dna_type, folder, trajectory, log, energy_term, periodic_size, platform",
+        [(row['Name'], row['DNA type'], row['Folder'], row['Trajectory'], row['Log'],
+          row['Energy term'], row['periodic size'], row['Platform']) for _, row in _test_sets.iterrows()],
+        ids=id_list  # Use the function here
+    )
+    def test_energies(self, name, dna_type, folder, trajectory, log, energy_term, periodic_size, platform, skip_platform):
+        if skip_platform and skip_platform.lower() == platform.lower():
+            pytest.skip(f"Skipping tests for platform: {platform}")
+        folder = Path(folder)
+        dna = DNA.fromXYZ(f'{folder}/in00_conf.xyz', dna_type, template_from_X3DNA=False)
+        system = System(dna)
+        self._test_energy(energy_term, folder/log, folder/trajectory, name, periodic_size, platform, dna, system)
+
+    @pytest.mark.parametrize(
+        "name, dna_type, folder, trajectory, log, energy_term, periodic_size, platform",
+        [(row['Name'], row['DNA type'], row['Folder'], row['Trajectory'], row['Log'],
+          row['Energy term'], row['periodic size'], row['Platform']) for _, row in _test_sets.iterrows()],
+        ids=id_list  # Use the function here
+    )
+    def test_forces(self, name, dna_type, folder, trajectory, log, energy_term, periodic_size, platform, skip_platform):
+        if skip_platform and skip_platform.lower() == platform.lower():
+            pytest.skip(f"Skipping tests for platform: {platform}")
+        folder = Path(folder)
+        dna = DNA.fromXYZ(f'{folder}/in00_conf.xyz', dna_type, template_from_X3DNA=False)
+        system = System(dna)
+        self._test_force(energy_term, folder/log, folder/trajectory, name, periodic_size, platform, dna, system)
+
+  
+    def _test_energy(self,
+                     log_energy='E_bond',
+                     log_file='Tests/adna/sim.log',
+                     traj_file='Tests/adna/traj.xyz',
+                     force='Bond', periodic_size=94.2,
+                     platform_name='Reference', dna=None, system=None):
+        self.dna = dna
+        self.system = system
+        self.system.clearForces()
+        self.system.setDefaultPeriodicBoxVectors(*np.diag([2 * periodic_size / 10] * 3))
+
+        log = parse_log(log_file)
+
+        self.system.clearForces()
+        f = forces[force]
+        tempforce = f(self.dna)
+        try:
+            tempforce.addForce(self.system)
+        except AttributeError:
+            self.system.addForce(tempforce)
+        energies = self.system.recomputeEnergy(traj_file, platform_name=platform_name)
+        d = (energies / _ef - log[log_energy])
+        diff = np.sqrt((d ** 2).sum() / len(energies))
+        print(f'The difference in the energy term {log_energy} is {diff} Kcal/mol')
+        print(f'The DNA type of the system is {self.dna.DNAtype}')
+        data = np.array([energies / _ef, np.array(log[log_energy]), np.array(d)])
+        results = pandas.DataFrame(data.T, columns=['Openmm energy', 'Lammps energy', 'Difference'])
+        # print(data)
+        print(pandas.DataFrame(data.T, columns=['Openmm energy', 'Lammps energy', 'Difference']))
+        # The maximum error seems to be bonds in curved BDNA (0.002)
+        assert diff < 2E-3, diff
+        assert len(results.dropna()) == len(results), results
+
+    def _test_force(self,
+                    log_energy='E_bond',
+                    log_file='Tests/adna/sim.log',
+                    traj_file='Tests/adna/traj.xyz',
+                    force='Bond', periodic_size=94.2,
+                    platform_name='Reference', dna=None, system=None):
+        self.dna = dna
+        self.system = system
+        self.system.clearForces()
+        self.system.setDefaultPeriodicBoxVectors(*np.diag([2 * periodic_size / 10] * 3))
+        log = parse_log(log_file)
+        self.system.clearForces()
+        f = forces[force]
+        tempforce = f(self.dna)
+        try:
+            tempforce.addForce(self.system)
+        except AttributeError:
+            self.system.addForce(tempforce)
+        temperature = 300 * openmm.unit.kelvin
+        integrator = openmm.LangevinIntegrator(temperature, 1 / openmm.unit.picosecond,
+                                                     2 * openmm.unit.femtoseconds)
+        platform = openmm.Platform.getPlatformByName(platform_name)
+        simulation = openmm.app.Simulation(self.system.top, self.system, integrator, platform)
+        simulation.context.setPositions(self.system.coord.getPositions())
+        energy_unit = openmm.unit.kilojoule_per_mole
+        nan_force_particles = 0
+        for i in range(10):
+            state = simulation.context.getState(getForces=True)
+            ff = state.getForces()
+            sf = (np.array(ff) ** 2).sum(axis=1) ** .5
+            for j, f in enumerate(sf):
+                if np.isnan(f.value_in_unit(unit.kilojoule_per_mole / unit.nanometer)):
+                    print(f"Particle {j + 1}/{len(sf)} has force {f} at step {i}")
+                    nan_force_particles += 1
+            simulation.step(1)
+        assert nan_force_particles == 0, "At least one particle has undefined force"
+
+    def _test_energies_slow(self):
+        test_sets = pandas.read_csv('Tests/test_cases.csv', comment='#')
+        for i, test in test_sets.iterrows():
+            dna_type = test['DNA type']
+            folder = test['Folder']
+            yield self._test_energy, test['Energy term'], f'{folder}/{test.Log}', f'{folder}/{test.Trajectory}', \
+                  test['Name'], folder, dna_type, test['periodic size'], test['Platform']
+
+if __name__=='__main__':
+    skip_platform=False
+    for _, (name, dna_type, folder, trajectory, log, energy_term, periodic_size, platform) in _test_sets.iterrows():
+        test=TestEnergies()
+        test.test_energies(name, dna_type, folder, trajectory, log, energy_term, periodic_size, platform, skip_platform)
+        test.test_forces(name, dna_type, folder, trajectory, log, energy_term, periodic_size, platform, skip_platform)
\ No newline at end of file