Reading smallbox.nc file for system mp0.
smallbox.nc does not contain coordinates data.
smallbox.nc does not contain forces data.
smallbox.nc does not contain reference energies data.
Data of system mp0 was read from file smallbox.nc
smallbox.zip

Hi,
I am trying to make use of this wonderful tool, but I keep getting errors for the provided examples, mostly because of the dependence of ParaMol on many external helper codes, like ASE, DFTB+, etc. I wonder if you have recently tested the examples and if you can list the exact versions of the auxiliary codes you have used and if you can provide some interactive Jupiter notebooks with these examples,

Thanks,
Marawan

Hi,
I am trying to optimize an existing ligand amber parameter topology file using the adaptive parameterization scheme (DFTB+). At the updated parameters writing step, I am observing that the saved updated amber prmtop file stays the same as the original prmtop file, even though the 'lig_adaptive_param.ff' got updated, i.e. it is slightly different from the file saved before the adaptive parameterization. I am wondering if I am missing something. Here is my script:

#########################
import time
from time import process_time
t1_start = process_time()

print(t1_start)

from ParaMol.System.system import *
from ParaMol.MM_engines.openmm import *
from ParaMol.QM_engines.qm_engine import *
from ParaMol.Tasks.adaptive_parametrization import *
from ParaMol.Utils.settings import *
from ParaMol.Utils.amber_symmetrizer import *

---------------------------------------------------------

Preparation

---------------------------------------------------------

Create the OpenMM engine for lig

openmm_system = OpenMMEngine(init_openmm=True, topology_format='AMBER', top_file='lig.prmtop', crd_file='lig.inpcrd')

Create ParaMol System

lig = ParaMolSystem(name="lig", engine=openmm_system, n_atoms=67, n_cpus=4)

Create ParaMol's force field representation and ask to parametrize bonds, angles and torsions

lig.force_field.create_force_field(opt_bonds=True, opt_angles=True, opt_torsions=True)

Create ParaMol settings instance

paramol_settings = Settings()

The objective function will contain a energy, force and regularization term

paramol_settings.properties["include_energies"] = True
paramol_settings.properties["include_forces"] = True
paramol_settings.properties["include_regularization"] = True

Create the ASE calculator

from ase.calculators.dftb import *

calc = Dftb(Hamiltonian_='DFTB', # line is included by default
Hamiltonian_MaxAngularMomentum_='',
Hamiltonian_MaxAngularMomentum_H='s',
Hamiltonian_MaxAngularMomentum_O='p',
Hamiltonian_MaxAngularMomentum_C='p',
Hamiltonian_MaxAngularMomentum_N="p",
Hamiltonian_MaxAngularMomentum_S="p",
Hamiltonian_MaxAngularMomentum_Br='p',
Hamiltonian_MaxAngularMomentum_F='p',
Hamiltonian_MaxAngularMomentum_Cl='p',
Hamiltonian_MaxAngularMomentum_I='p',
Hamiltonian_MaxAngularMomentum_P='p',
Hamiltonian_Dispersion="DftD3 { \n s6=1.000 \n s8=0.5883 \n Damping = BeckeJohnson { \n a1=0.5719 \n a2=3.6017 \n } \n }",
Hamiltonian_SCC='Yes',
Hamiltonian_SCCTolerance=1e-8,slako_dir='~/dftbplus-20.2.1/slako/3ob/3ob-3-1/')

Alternative, we could set the calculator in the settings

paramol_settings.qm_engine["qm_engine"] = "ase"
paramol_settings.qm_engine["ase"]["calculator"] = calc

amber_symmetrizer = AmberSymmetrizer("lig.prmtop")
amber_symmetrizer.get_amber_symmetries(lig.force_field)

Write symmetrized Force-Field to file

lig.force_field.write_ff_file("lig_sym.ff")

---------------------------------------------------------

Adaptive Parametrization

---------------------------------------------------------

adaptive_parametrization = AdaptiveParametrization()
systems, parameter_space, objective_function, optimizer = adaptive_parametrization.run_task(paramol_settings, [lig], rmsd_tol=0.01, max_iter=50, structures_per_iter=100, )

Write ParaMol Force Field file with final parameters

lig.force_field.write_ff_file("lig_adaptive_param.ff")

Update AMBER symmetrizer with new parameters

amber_symmetrizer.update_term_types_parameters(parameter_space.optimizable_parameters)

Write AMBER topology file (.prmtop)

amber_symmetrizer.save_prmtop("lig_opt.prmtop")
amber_symmetrizer.save_frcmod("lig_opt.frcmod")

t1_stop = process_time()
print("Elapsed time:", t1_stop, t1_start)
print("Elapsed time during the whole program in hours:",(t1_stop-t1_start)/3600)
#########################

Any clue,

Thanks,
Marawan

Hi I'm having some trouble running ParaMol in parallel mode:
objective_function uses OpenCL as default, even if CUDA is given
platform properties (like OpenCLPlatformIndex) cannot be passed via objective_function presumably leading to the error
--------------------------------------------------------------------------- RemoteTraceback Traceback (most recent call last) RemoteTraceback: """ Traceback (most recent call last): File "/home/ac127777/anaconda3/envs/mdanalysis-dev/lib/python3.8/multiprocessing/pool.py", line 125, in worker result = (True, func(*args, **kwds)) File "/home/ac127777/anaconda3/envs/mdanalysis-dev/lib/python3.8/multiprocessing/pool.py", line 51, in starmapstar return list(itertools.starmap(args[0], args[1])) File "/home/ac127777/anaconda3/envs/mdanalysis-dev/lib/python3.8/site-packages/ParaMol/Objective_function/gpu_objective_function.py", line 80, in f_callback_function context = Context(system, copy.deepcopy(integrator), platform) File "/home/ac127777/anaconda3/envs/mdanalysis-dev/lib/python3.8/site-packages/simtk/openmm/openmm.py", line 13232, in __init__ _openmm.Context_swiginit(self, _openmm.new_Context(*args)) simtk.openmm.OpenMMException: Specified DeviceIndex but not OpenCLPlatformIndex. When multiple platforms are available, a platform index is needed to specify a device. """
`The above exception was the direct cause of the following exception:

OpenMMException Traceback (most recent call last)
in
7 adaptive_parametrization = AdaptiveParametrization()
8 os.system('rm -rf ASEWorkDir_mp0_test')
----> 9 adaptive_parametrization.run_task(paramol_settings, [mp0_test2], rmsd_tol=0.0001, max_iter=2, structures_per_iter=1)

~/anaconda3/envs/mdanalysis-dev/lib/python3.8/site-packages/ParaMol/Tasks/adaptive_parametrization.py in run_task(self, settings, systems, parameter_space, objective_function, optimizer, interface, rmsd_tol, max_iter, steps_integrator, structures_per_iter, wham_reweighing, restart)
192
193 # Perform parametrization
--> 194 systems, parameter_space, objective_function, optimizer = parametrization.run_task(settings=settings,
195 systems=systems,
196 parameter_space=parameter_space,

~/anaconda3/envs/mdanalysis-dev/lib/python3.8/site-packages/ParaMol/Tasks/parametrization.py in run_task(self, settings, systems, parameter_space, objective_function, optimizer, interface, adaptive_parametrization, apply_charge_correction, restart)
107 '''
108 # Print Initial Info of Objective Function
--> 109 objective_function.f(parameter_space.optimizable_parameters_values_scaled, opt_mode=False)
110
111 # Create optimizer

~/anaconda3/envs/mdanalysis-dev/lib/python3.8/site-packages/ParaMol/Objective_function/objective_function.py in f(self, parameters_values, opt_mode)
318
319 if self._parallel:
--> 320 fmm, emm = run_parallel()
321 esp = None
322 else:

~/anaconda3/envs/mdanalysis-dev/lib/python3.8/site-packages/ParaMol/Objective_function/objective_function.py in run_parallel()
294 with mp.get_context("fork").Pool(processes=self._total_n_batches,
295 initializer=self._parallel_function.init_worker) as pool:
--> 296 data = pool.starmap(self._parallel_function.f_callback_function, args)
297 pool.terminate()
298

~/anaconda3/envs/mdanalysis-dev/lib/python3.8/multiprocessing/pool.py in starmap(self, func, iterable, chunksize)
370 func and (a, b) becomes func(a, b).
371 '''
--> 372 return self._map_async(func, iterable, starmapstar, chunksize).get()
373
374 def starmap_async(self, func, iterable, chunksize=None, callback=None,

~/anaconda3/envs/mdanalysis-dev/lib/python3.8/multiprocessing/pool.py in get(self, timeout)
769 return self._value
770 else:
--> 771 raise self._value
772
773 def _set(self, i, obj):

OpenMMException: Specified DeviceIndex but not OpenCLPlatformIndex. When multiple platforms are available, a platform index is needed to specify a device.`

This is my input:

`import matplotlib.pyplot as plt
import numpy as np

ParaMol imports

from ParaMol.Utils.settings import *
from ParaMol.System.system import *
from ParaMol.MM_engines.openmm import *
from ParaMol.QM_engines.qm_engine import *
from ParaMol.Tasks.adaptive_parametrization import *
from ParaMol.Utils.conformational_sampling import conformational_sampling

openmm_system = OpenMMEngine(init_openmm=True, topology_format='CHARMM', top_file='sys_psf_fixed_edit.psf', crd_format='PDB', charmm_param_file='mp0.str', crd_file='sys_fixed_box3.pdb', platform_name='OpenCL', platform_properties = {'Precision': 'single','DeviceIndex': '0', 'OpenCLPlatformIndex': '0'}, periodic_sys=True, x=1.65, y=1.65, z=1.65, angle_x=90.0, angle_y=90.0, angle_z=90.0, create_system_params={'nonbondedMethod': app.PME, 'nonbondedCutoff': 0.8 * unit.nanometer, 'switchDistance': 0.7*unit.nanometer, 'constraints': app.HBonds, 'rigidWater': True} , integrator_params={'temperature': 310.0 * unit.kelvin, 'stepSize' : 0.001 * unit.picoseconds, 'frictionCoeff' : 5.0 / unit.picoseconds})

Create ParaMol System

mp0_test2 = ParaMolSystem(name="mp0_test", engine=openmm_system, n_atoms=332, n_cpus=4)

Create ParaMol's force field representation and ask to parametrize bonds, angles and torsions

mp0_test2.force_field.create_force_field(opt_bonds=True, opt_angles=True, opt_torsions=True, opt_lj=True, opt_sc=True)

Create ParaMol settings instance

paramol_settings = Settings()

The objective function will contain a energy, force, electrostat. potential, and regularization term

paramol_settings.properties["include_energies"] = True
paramol_settings.properties["include_forces"] = True
paramol_settings.properties["include_regularization"] = True

from ase.calculators.cp2k import CP2K

CP2K.command = 'env OMP_NUM_THREADS=4 cp2k_shell.ssmp'

cp2k config input

inp = """
&GLOBAL
RUN_TYPE ENERGY
&END GLOBAL

&FORCE_EVAL
METHOD Quickstep ! Electronic structure method (DFT,...)
&DFT
BASIS_SET_FILE_NAME ../../cp2k_geom_opts/mp0/PBE_TZV2P_M_D3/BASIS_MOLOPT
POTENTIAL_FILE_NAME ../../cp2k_geom_opts/mp0/PBE_TZV2P_M_D3/POTENTIAL
CHARGE 0
&MGRID
NGRIDS 5
CUTOFF 750
REL_CUTOFF 80
&END MGRID
&QS
METHOD GPW
EPS_DEFAULT 1.0E-12
&END QS
&POISSON ! Solver requested for non periodic calculations
PERIODIC XYZ
PSOLVER PERIODIC ! Type of solver
&END POISSON
&SCF ! Parameters controlling the convergence of the scf. This section should not be changed.
SCF_GUESS ATOMIC
EPS_SCF 1.0E-5
MAX_SCF 300
&MIXING
ALPHA 0.4
&END MIXING
&OT
MINIMIZER CG
&END OT
&END SCF
&XC ! Parameters needed to compute the electronic exchange potential
&VDW_POTENTIAL
DISPERSION_FUNCTIONAL PAIR_POTENTIAL
&PAIR_POTENTIAL
TYPE DFTD3
PARAMETER_FILE_NAME ../../cp2k_geom_opts/mp0/PBE_TZV2P_M_D3/dftd3.dat
REFERENCE_FUNCTIONAL PBE
CALCULATE_C9_TERM TRUE
REFERENCE_C9_TERM TRUE
&END PAIR_POTENTIAL
&END VDW_POTENTIAL
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&END XC
&END DFT

&SUBSYS
&CELL
SYMMETRY CUBIC
ABC 16.5 16.5 16.5
ALPHA_BETA_GAMMA 90.0 90.0 90.0
&END CELL
&KIND H
ELEMENT H
BASIS_SET DZVP-MOLOPT-GTH-q1
POTENTIAL GTH-PBE-q1
&END KIND
&KIND C
ELEMENT C
BASIS_SET DZVP-MOLOPT-GTH-q4
POTENTIAL GTH-PBE-q4
&END KIND
&KIND P
ELEMENT P
BASIS_SET DZVP-MOLOPT-GTH-q5
POTENTIAL GTH-PBE-q5
&END KIND
&KIND O
ELEMENT O
BASIS_SET DZVP-MOLOPT-GTH-q6
POTENTIAL GTH-PBE-q6
&END KIND
&PRINT
&ATOMIC_COORDINATES LOW
&END ATOMIC_COORDINATES
&INTERATOMIC_DISTANCES LOW
&END INTERATOMIC_DISTANCES
&KINDS
POTENTIAL
&END KINDS
&TOPOLOGY_INFO
PSF_INFO
&END TOPOLOGY_INFO
&END PRINT
&END SUBSYS

&PRINT
&FORCES ON
&END FORCES
&END PRINT

&END FORCE_EVAL

&MOTION
&GEO_OPT
OPTIMIZER LBFGS
MAX_DR 3.0E-03
MAX_FORCE 4.5E-04
RMS_DR 1.5E-03
RMS_FORCE 3.0E-04
&END
&CONSTRAINT
&FIXED_ATOMS
COMPONENTS_TO_FIX XYZ
LIST 12..332
&END
&END CONSTRAINT
&PRINT
&FORCES
UNIT N
FORMAT XYZ
FILENAME geom_opt_forces
&END FORCES
&END PRINT
&END MOTION
"""

ase calc config

calc = CP2K(basis_set=None,
basis_set_file=None,
max_scf=None,
charge=None,
cutoff=None,
force_eval_method=None,
potential_file=None,
poisson_solver=None,
pseudo_potential=None,
stress_tensor=False,
uks=False,
xc=None,
inp=inp)

from scipy.optimize import BFGS as BFGS_scipy

paramol_settings.optimizer['scipy']['method'] = BFGS_scipy

paramol_settings.qm_engine["ase"]['calculator'] = calc

paramol_settings.objective_function["parallel"] = True # should be true if n_cpus > 1

paramol_settings.objective_function["platform_name"] = 'OpenCL' # = 'CUDA' funnily gives the exact same error, also specifying platform properties for init_openmm is apparently independent, also doesn't change the error.

Make a copy of the original force_field dict.

It will be useful later to assess the quality of the re-parametrized force field.

force_field_original = copy.deepcopy(mp0_test2.force_field.force_field)

---------------------------------------------------------

Adaptive Parametrization

---------------------------------------------------------

adaptive_parametrization = AdaptiveParametrization()
adaptive_parametrization.run_task(paramol_settings, [mp0_test2], rmsd_tol=0.0001, max_iter=2, structures_per_iter=1)`

printout before error:

`# --------------------------------------------------------- #

Adaptive Parametrization

---------------------------------------------------------

adaptive_parametrization = AdaptiveParametrization()

#os.system('rm -rf ASEWorkDir_mp0_test')

adaptive_parametrization.run_task(paramol_settings, [mp0_test2], rmsd_tol=0.0001, max_iter=2, structures_per_iter=1)

!=================================================================================!
! ADAPTIVE PARAMETRIZATION !
!=================================================================================!
mp0_test
Number of parallel threads inherited from n_cpus: 4
Number of parallel threads set to: 1.
Please set cp2k threads using CP2K.command
!=================================================================================!
! Scaling Constants !
! Term type Value !
!---------------------------------------------------------------------------------!
! bond_eq 0.10009122 !
! bond_k 364230.66156522 !
! angle_eq 1.83622943 !
! angle_k 462.52405246 !
! torsion_phase 0.39269908 !
! torsion_k 2.40580000 !
! lj_sigma 0.13566853 !
! lj_eps 0.34719513 !
! scee 1.00000000 !
! scnb 0.83951436 !
!=================================================================================!
!=================================================================================!
! Prior Widths !
! Term type Value !
!---------------------------------------------------------------------------------!
! bond_eq 0.05000000 !
! bond_k 100000.00000000 !
! angle_eq 0.19634954 !
! angle_k 100.00000000 !
! torsion_phase 3.14159265 !
! torsion_k 16.73600000 !
! lj_sigma 0.30000000 !
! lj_eps 0.20000000 !
! scee 1.00000000 !
! scnb 1.00000000 !
!=================================================================================!
Iteration no. 0 of the adaptive parametrization loop.
Generating new configuration 1.
Atoms(symbols='CH219O111P', pbc=True, cell=[16.5, 16.5, 16.5], calculator=CP2K(...))
Generated new 1 MM structures for system mp0_test.

Writing restart_paramol/mp0_test_data_restart.nc file for system mp0_test.
Data of system mp0_test was written to file restart_paramol/mp0_test_data_restart.nc
!=================================================================================!
! PARAMETRIZATION !
!=================================================================================!
`
(I have modified the MM_engines.openmm to support my periodic box system and the QM_engines.ase_warpper
to avoid the deepcopy of my calculator object which has a thread lock in it somewhere)
modified_engines.zip

Hi, I did a torsion scan with my own molecule with paramol but the following error occured:

Traceback (most recent call last):
  File "/root/repo/prepare/paramol.py", line 85, in <module>
    systems, parameter_space, objective_function, optimizer = parametrization.run_task(paramol_settings, [norfloxacin])
  File "/data/miniconda3/envs/uii/lib/python3.7/site-packages/ParaMol/Tasks/parametrization.py", line 109, in run_task
    objective_function.f(parameter_space.optimizable_parameters_values_scaled, opt_mode=False)
  File "/data/miniconda3/envs/uii/lib/python3.7/site-packages/ParaMol/Objective_function/objective_function.py", line 325, in f
    objective_function = self.calculate_objective_function(parameters_values, fmm, emm, esp)
  File "/data/miniconda3/envs/uii/lib/python3.7/site-packages/ParaMol/Objective_function/objective_function.py", line 196, in calculate_objective_function
    energy_decomposed = property.calculate_property(emm)
  File "/data/miniconda3/envs/uii/lib/python3.7/site-packages/ParaMol/Objective_function/Properties/energy_property.py", line 109, in calculate_property
    obj_fun_energies.append(obj_fun / var)
FloatingPointError: invalid value encountered in double_scalars

I've print the obj_fun and var, they were 0.

Here is the code and files I've used, which was basically copied from the example.
RB2345_AC.zip

    openmm_system = OpenMMEngine(init_openmm=True, topology_format='AMBER',
                                    
 top_file='RB2345_AC.prmtop', 
                                    crd_format='AMBER',
                                    crd_file='RB2345_AC.inpcrd')
    # Create ParaMol System
    norfloxacin = ParaMolSystem(name="norfloxacin", engine=openmm_system, n_atoms=54)

    # Create ParaMol's force field representation
    norfloxacin.force_field.create_force_field()

    # Create ParaMol settings instance
    paramol_settings = Settings()

    # The objective function will contain an energy and a regularization term
    paramol_settings.properties["include_energies"] = True
    paramol_settings.properties["include_forces"] = False # One should not include forces when a torsional scan
    paramol_settings.properties["include_regularization"] = True

    # --------------------------------------------------------- #
    #                Symmetrize ParaMol ForceField              #
    # --------------------------------------------------------- #
    # Symmetry ParaMol ForceField so that it respects atom-type symmetries
    amber_symmetrizer = AmberSymmetrizer(top_file="RB2345_AC.prmtop")
    amber_symmetrizer.get_symmetries()
    amber_symmetrizer.symmetrize_force_field(norfloxacin.force_field)

    # --------------------------------------------------------- #
    #                     Set the QM Engine                     #
    # --------------------------------------------------------- #
    # Create the ASE calculator
    from ase.calculators.dftb import *

    calc = Dftb(Hamiltonian_='DFTB',  # line is included by default
                Hamiltonian_MaxAngularMomentum_='',
                Hamiltonian_MaxAngularMomentum_H='s',
                Hamiltonian_MaxAngularMomentum_O='p',
                Hamiltonian_MaxAngularMomentum_C='p',
                Hamiltonian_MaxAngularMomentum_N="p",
                Hamiltonian_MaxAngularMomentum_F="p",
                Hamiltonian_Dispersion="DftD3 { \n s6=1.000 \n s8=0.5883 \n Damping = BeckeJohnson { \n a1=0.5719 \n a2=3.6017 \n } \n }",
                Hamiltonian_SCC='Yes',
                Hamiltonian_SCCTolerance=1e-8, 
                slako_dir='/data/halorg-0-1/')

    # Alternative, we could set the calculator in the settings
    paramol_settings.qm_engine["ase"]["calculator"] = calc

    # --------------------------------------------------------- #
    #                  Perform the Torsional Scan               #
    # --------------------------------------------------------- #
    torsion_to_scan = [[27, 28, 29, 30]]
    scan_settings = [[-180.0, 180.0, 90.0]]
    torsion_scan = TorsionScan()
    torsion_scan.run_task(paramol_settings, [norfloxacin], torsion_to_scan, scan_settings, optimize_qm_before_scan=True)

    # Save coordinates and energies into .nc file
    norfloxacin.write_data("norfloxacin_scan.nc")
    # --------------------------------------------------------- #
    #                   Parametrize the Torsion                 #
    # --------------------------------------------------------- #
    norfloxacin.force_field.optimize_torsions_by_symmetry(torsion_to_scan)
    norfloxacin.force_field.write_ff_file("norfloxacin_symm.ff")

    parametrization = Parametrization()
    systems, parameter_space, objective_function, optimizer = parametrization.run_task(paramol_settings, [norfloxacin])

    # Update AMBER symmetrizer with new parameters
    amber_symmetrizer.update_term_types_parameters(parameter_space.optimizable_parameters)

    # Write AMBER topology file (.prmtop) and .frcmod file
    amber_symmetrizer.save("norfloxacin_opt.prmtop")
    amber_symmetrizer.save_frcmod("norfloxacin_opt.frcmod")

Hi, I'm currently doing a torsion scan to improve the dihedrals of some molecules but I'm getting an error regarding the units of the position of the coordinates. I've tried a couple of approaches to fix the problem because the code was working with an earlier release in April or May. I'm currently trying to fix the issue using a diphenylmethane but I haven't been successful. Any input as to what the issue maybe would be useful. Below I pasted the error that I'm getting. Thanks!!

BFGS: 60 09:36:31 -13673.084137 0.0124
BFGS: 61 09:36:31 -13673.084155 0.0096

TypeError Traceback (most recent call last)
/ipykernel_22881/1791874995.py in
14 mol, paramol_settings, amber_symmetrizer = generate_fresh_system(mol_name)
15 torsion_scan = TorsionScan()
---> 16 torsion_scan.run_task(paramol_settings,
17 [mol],
18 [torsion_to_scan],

python3.8/site-packages/ParaMol/Tasks/torsions_scan.py in run_task(self, settings, systems, torsions_to_scan, scan_settings, interface, torsions_to_freeze, ase_constraints, optimize_qm, optimize_qm_before_scan, optimize_mm, optimize_mm_before_scan, optimize_mm_type, sampling, rotate_from_initial, n_structures_to_sample, dihedral_conservation_threshold, mm_opt_force_constant, mm_opt_tolerance, mm_opt_max_iter, rdkit_conf, restart)
170 if torsional_scan_dim == 1:
171 # Perform 1D Scan
--> 172 scan_angles, qm_energies_list, qm_forces_list, qm_conformations_list, mm_energies_list, mm_conformations_list = self.scan_1d(
173 interface, settings.restart, system, conf, torsions_to_scan[0], torsions_to_freeze, scan_settings[0], optimize_qm, optimize_qm_before_scan,
174 optimize_mm, optimize_mm_before_scan, optimize_mm_type, ase_constraints, sampling, rotate_from_initial, n_structures_to_sample, dihedral_conservation_threshold,

python3.8/site-packages/ParaMol/Tasks/torsions_scan.py in scan_1d(self, interface, restart_settings, system, rdkit_conf, torsion_to_scan, torsions_to_freeze, scan_settings, optimize_qm, optimize_qm_before_scan, optimize_mm, optimize_mm_before_scan, optimize_mm_type, ase_constraints, sampling, rotate_from_initial, n_structures_to_sample, threshold, mm_opt_force_constant, mm_opt_tolerance, mm_opt_max_iter, restart)
468
469 # Set positions of context to last position
--> 470 dummy_context.setPositions(positions * unit.nanometers)
471
472 # Set RDKit geometry to the current in the OpenMM context

python3.8/site-packages/openmm/openmm.py in setPositions(self, positions)
3676
3677 if unit.is_quantity(positions):
-> 3678 positions = positions.value_in_unit(unit.nanometer)
3679
3680

python3.8/site-packages/openmm/unit/quantity.py in value_in_unit(self, unit)
621 Returns underlying value, in the specified units.
622 """
--> 623 val = self.in_units_of(unit)
624 if is_quantity(val):
625 return val._value

python3.8/site-packages/openmm/unit/quantity.py in in_units_of(self, other_unit)
657 """
658 if not self.unit.is_compatible(other_unit):
--> 659 raise TypeError('Unit "%s" is not compatible with Unit "%s".' % (self.unit, other_unit))
660 f = self.unit.conversion_factor_to(other_unit)
661 return self._change_units_with_factor(other_unit, f)

TypeError: Unit "nanometer**2" is not compatible with Unit "nanometer".

I'm interested in using ParaMol for small molecule parameterization. For molecule such as tert-butyl alcohol, the three methyl group should have the same partial charge. I wonder how should one implement this in the RESP routine? Thank you.