diff --git a/.gitignore b/.gitignore index b806111..447b710 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,7 @@ __pycache__ +build/ +doc/_build/ +doc/generated/ +public/ *.egg-info *.h5 diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml new file mode 100644 index 0000000..57a2171 --- /dev/null +++ b/.gitlab-ci.yml @@ -0,0 +1,22 @@ +image: python:3.11 + +test: + stage: test + script: + - pip install -U sphinx + - pip install . + - sphinx-build -b html doc public + rules: + - if: $CI_COMMIT_REF_NAME != $CI_DEFAULT_BRANCH + +pages: + stage: deploy + script: + - pip install -U sphinx + - pip install . + - sphinx-build -b html doc public + artifacts: + paths: + - public + rules: + - if: $CI_COMMIT_REF_NAME == $CI_DEFAULT_BRANCH diff --git a/.readthedocs.yaml b/.readthedocs.yaml new file mode 100644 index 0000000..aa209bb --- /dev/null +++ b/.readthedocs.yaml @@ -0,0 +1,32 @@ +# .readthedocs.yaml +# Read the Docs configuration file +# See https://docs.readthedocs.io/en/stable/config-file/v2.html for details + +# Required +version: 2 + +# Set the OS, Python version and other tools you might need +build: + os: ubuntu-22.04 + tools: + python: "3.11" + # You can also specify other tool versions: + # nodejs: "19" + # rust: "1.64" + # golang: "1.19" + +# Build documentation in the "docs/" directory with Sphinx +sphinx: + configuration: docs/conf.py + +# Optionally build your docs in additional formats such as PDF and ePub +# formats: +# - pdf +# - epub + +# Optional but recommended, declare the Python requirements required +# to build your documentation +# See https://docs.readthedocs.io/en/stable/guides/reproducible-builds.html +# python: +# install: +# - requirements: docs/requirements.txt diff --git a/LICENSE.md b/LICENSE.md index 0d468c4..ac70fea 100644 --- a/LICENSE.md +++ b/LICENSE.md @@ -1,15 +1,21 @@ # CMIclassirot licensing conditions -CMIclassirot is provided by the CFEL Controlled Molecule Imaging group as is. It is licensed under -the [GPL v3](./LICENSE-GPLv3.md) with the following additional requirements: +CMIclassirot is provided by the CFEL Controlled Molecule Imaging group as is and without any +warranty implied. It is licensed under the [GPL v3](./LICENSE-GPLv3.md) with the following +additional requirements: * Its use for scientific work is acknowledged by the appropriate reference in any resulting work, e.g., scientific publication. + For current work, please cite the preprint Muhamed Amin, Jean-Michel Hartmann, Amit K. Samanta, + Jochen Küpper, ''Laser-induced alignment of nanoparticles and macromolecules for + single-particle-imaging applications'', + [arXiv:2306.05870 \[physics\]](https://arxiv.org/abs/2306.05870) + * All corrections and enhancements must be contributed back to the development of the package within an appropriate time, i.e., no later than the publication of its first use. - Please create a pull request on GitHub (preferred) or send an appropriate patch against the latest - program version on the `develop` branch. + + Please send us an appropriate patch against the latest program version on the `develop` branch. See the documentation for a full list of contributors. diff --git a/README.devel.md b/README.devel.md index dfbb434..20eae14 100644 --- a/README.devel.md +++ b/README.devel.md @@ -1,12 +1,22 @@ # Development guidelines -Please follow some simple rules for coding on CMIstark: -- keep the code consistent across the whole package! -- use git-flow (the tools and the concept) -- consistently use two empty lines between functions/methods and three empty +Please follow some simple rules for coding on CMIclassirot: + +* keep the code consistent across the whole package! +* Document all code and adjust comments before you change the code to reflect + and discuss, even if only with yourself, the planned changes. +* use git-flow (possibly the tools and especially the concept) +* consistently use two empty lines between functions/methods and three empty lines above a class +## Developer install + +You can install the project in _editable_ mode using +``` +pip install --editable . +``` + diff --git a/Release Notes.md b/Release Notes.md new file mode 100644 index 0000000..7628dd9 --- /dev/null +++ b/Release Notes.md @@ -0,0 +1,26 @@ +# CMIclassirot release notes + +## Future Release(s) -- past 1.0 + +* Provide $`\left<\cos^2\theta\right>(t)`$ trace in HDF5 files +* Add examples for bionanoparticles + + +## Release 0.9.2 + +* Improve documentation + + +## Release 0.9.1 + +* Add installed `cmiclassirot` and `cmiclassirot-plot` user scripts to drive + operation +* Implement clean and documented examples +* Performance improvements (Field) +* Implement sphinx/online documentation +* Start documenting the code and the module + + +## Release 0.9 + +* Original functionality as described in arXiv:2306.05870 diff --git a/bin/cmiclassirot b/bin/cmiclassirot new file mode 100644 index 0000000..12ea471 --- /dev/null +++ b/bin/cmiclassirot @@ -0,0 +1,58 @@ +#!/usr/bin/env python +# -*- coding: utf-8; fill-column: 120 -*- +# +# This file is part of the CMIclassirot classical-rotation alignment simulations + + + +import click +import time + +from cmiclassirot.propagate import Propagate + + +@click.command() # help='Filename of definiton of Ensemble and Field', +@click.argument('inputfilename', required=1) +@click.option('-o', '--output', 'output', default='cmiclassirot.h5', show_default=True, + help='Write output to specified filename.') +@click.help_option('-h', '--help') +def main(inputfilename, output): + """CMIclassirot driver program: calculate the time-evolution of rigid rotors in electric fields + + This program reads an imputfile defining an Ensemble of Molecules and a Field and performs the calculation. + + Besides the options defined below, the program requires the filename of the inputfile as the only positional + argument. The input must define the following Python variables, which are used in the calculation: + + :class:`Ensemble` :param:`ensemble` Definition of the initial ensemble of molecules + + :class:`Field` :param:`field` Definition of the time-dependet electric field + + :param timerange: 2-tuple of initial and final time of integration + + :param dt_save: Tiemstep for saving the :class:`Molecule`-positions to file + + @author: Jochen Küpper + + """ + # read specification of problem + with open(inputfilename, mode='r') as inputfile: + code = inputfile.read() + exec(code, globals()) + + # perform the computation + starttime = time.time() + print('Starting propagation of molecular dynamics') + p = Propagate(ensemble, field, timerange, dt_save) + print(' Propagation took', time.time()-starttime, 's') + + # and save the results to the output file + starttime = time.time() + print('Saving data to file') + ensemble._save(output) + print(' Saving took', time.time()-starttime, 's') + + + +if __name__ == '__main__': + main() diff --git a/bin/cmiclassirot-plot b/bin/cmiclassirot-plot new file mode 100755 index 0000000..32baddf --- /dev/null +++ b/bin/cmiclassirot-plot @@ -0,0 +1,41 @@ +#!/usr/bin/env python +# -*- coding: utf-8; fill-column: 120 -*- +# +# This file is part of the CMIclassirot classical-rotation alignment simulations +# +# Plotting of degree of alignment and laser pulse + +from cmiclassirot.field import Field +from cmiclassirot.sample import Ensemble +import matplotlib.pyplot as plt +import numpy as np +import sys + +fname = sys.argv[1] + +e = Ensemble.FromFile(fname) +print("Number of Molecules:", len(e.molecules)) + +doa = np.zeros(len(e.molecules[0].pos)) +E = e.pulse + +for i in e.molecules: + t=[] + pulse = [] + for j in range(len(i.pos)): + try: + doa[j] += i.pos[j].angle.rotation_matrix[2][2]**2 + t.append(i.pos[j].time) + pulse.append(E(i.pos[j].time)) + except: + pass +doa = doa/len(e.molecules) + +fig, ax = plt.subplots(2) +ax[0].plot(np.array(t)*1e9, doa) +ax[1].plot(np.array(t)*1e9, Field.amplitude2intensity(np.array(pulse)) / (1e12 * 1e4)) +plt.xlabel('time (ns)') +ax[0].set_ylabel(r'$\left<\cos^2\theta\right>$') +ax[1].set_ylabel('$E$ ($10^{12}$~W/cm$^2$)') +plt.tight_layout() +plt.show() diff --git a/cmiclassirot/field.py b/cmiclassirot/field.py index 5e137c2..e028fb7 100644 --- a/cmiclassirot/field.py +++ b/cmiclassirot/field.py @@ -21,26 +21,56 @@ from scipy.constants import c, epsilon_0 from pyquaternion import Quaternion + + class Field(object): """Electric field class - For the beginning, this only implements a simple Gaussian pulse (in field strengths!) centered - at t=0 and along `Z`. + For the beginning, this only implements a simple Gaussian pulse (in amplitude) centered at t=0 + and along `Z`. + + Must specify one of `peak_amplitude`, `peak_intensity`, or `fieldname` to provide an + ac-electric-field envelope. In the formre two cases, one must also specify the width as `FWHM` + in amplitude or intensity, respectively, and the time of the `peak` of the laser pulse. In the + latter case, `file_content` needs to be specified. + + :param peak_amplitude: Peak amplitude of a synthetic Gaussian laser pulse (V/m) + + :param peak_intensity: Peak intensity of a synthetic Gaussian laser pulse (W/cm^2) + + :param t_peak: Time of the laser-pulse peak (s) + + :param FWHM: Full-width-at-half-maximum width of laser pulse in same units as filed + specification, i.e., :param`peak_amplitude` or :param:`peak_intensity` (s) + + :param filename: Name of file to read numerically-specified field + + :param file_content: Specify if file contains `amplitude` (V/m) or `intensity` (W/c^2) values. + This functionality is currently not implemented and the file is always expwected to contain + amplitudes. Please improve. .. todo:: Keep in mind that we want to be able to represent elliptically polaerized field by their envelope. - """ + .. todo:: document class, constructor, and methods - def __init__(self, amplitude=1.e15, sigma=1.e-9, mean=5.e-9, - intensity=True, filename=None, dc=None): - self.amplitude = amplitude - self.sigma = sigma - self.Ez = np.array([0., 0., 1.]) - self.from_file = None - self.intensity = intensity - self.mu = mean - self.dc = dc - if filename is not None: + .. todo:: Implement mixed field handling, here this is always providing (DC, AC) tuples on + __call__; maybe should be done via a derived class MixedField for performance reasons. + + """ + def __init__(self, peak_amplitude=None, peak_intensity=None, t_peak=0., FWHM=10.e-9, + filename=None, file_content='intensity'): + # store field internally as an field apmplitude + assert (peak_amplitude or peak_intensity or filename) # at least one is not None + if peak_intensity: + self.peak_amplitude = Field.intensity2amplitude(peak_intensity) + self.sigma = FWHM / 2*np.sqrt(2*np.log(2)) + self.t_peak = t_peak + elif peak_amplitude: + self.peak_amplitude = peak_amplitude + self.sigma = FWHM / 2*np.sqrt(2*np.log(2)) + self.t_peak = t_peak + else: + self.peak_amplitude = None f = open(filename) E=[] t=[] @@ -50,38 +80,34 @@ def __init__(self, amplitude=1.e15, sigma=1.e-9, mean=5.e-9, E.append(float(token[1])) t = np.array(t) E = np.array(E) - self.from_file = interpolate.interp1d(t,E) + self.amplitude = interpolate.interp1d(t,E) + # initially the field is alog :math:`Z` + self.Ez = np.array([0., 0., 1.]) + def __call__(self, t): """Field at time t - :return: Field vector at time + :return: Field amplitude vector at time :math:`t` + """ - if self.dc is not None: - if t > self.dc[0] and t < self.dc[1]: - return self.convert(self.dc[2]) - else: - return 0 - - if self.from_file is None: - E = self.amplitude * np.exp(-0.5 * ((t-self.mu)/self.sigma)**2) - if self.intensity: - return self.convert(E) - else: - return E + if self.peak_amplitude: + return self.peak_amplitude * np.exp(-0.5 * ((t - self.t_peak)/self.sigma)**2) else: - try: - E = 1.e16*self.from_file(t) - if self.intensity: - return self.convert(E) - else: - return E - except: - return 0. - - def convert(self, I): - """ Converts intensity to electric field""" - return np.sqrt(2*I/(c * epsilon_0)) + return self.amplitude(t) + + + @staticmethod + def intensity2amplitude(I): + """Converts intensity (in W/m**2) to electric field (in V/m)""" + return np.sqrt(2 * I / (c * epsilon_0)) + + + @staticmethod + def amplitude2intensity(A): + """Converts electric field amplitude (in V/m) to intensity (in W/m**2)""" + return c * epsilon_0 / 2 * A**2 + def rotate(self, quaternion=Quaternion()): """Rotate field diff --git a/cmiclassirot/propagate.py b/cmiclassirot/propagate.py index ab233db..8eeb94d 100644 --- a/cmiclassirot/propagate.py +++ b/cmiclassirot/propagate.py @@ -17,38 +17,44 @@ # see . import numpy as np -from scipy.integrate import ode +import scipy.integrate import pyquaternion as quat -from cmiclassirot.sample import Position import multiprocessing as mp +from cmiclassirot.sample import Position + + class Propagate(object): """Propagate an Ensemble in time""" - def __init__(self, ensemble, field, dt=2.e-11, timerange=(-1e-9,1e-9)): + def __init__(self, ensemble, field, timerange=(0,1e-9), dt_save=None): """Initialize propagator - :param ensemble: Ensemble with all |Molecule|s to be propagated + :param ensemble: :class:`Ensemble` with all |Molecule|s to be propagated - :param field: Field in which to propagate. This is a univariate function that returns the - field strength (V/m) for any given time (s) in the `timerange`. + :param field: :class:`Field` in which to propagate. This is a univariate function that + returns the field strength (V/m) for any given time (s) in the :param:`timerange`. - :param timerange: Timerange over which to propagate, specified as tuple (t_init, t_final) + :param timerange: Time range over which to propagate, specified as tuple (t_init, t_final) + + :param dt_save: Timesteps for saving the current phase-space positions of the + :param:`Molecule`s during propagation (default: 1 % of timerange period) """ self.ensemble = ensemble self.field = field - self.dt = dt + if dt_save: + self.dt_save = dt_save + else: + self.dt_save = timerange[1] - timerange[0] self.t_range = timerange self.run() - def run(self): - """Propagate all |Molecule|s in the current |Field| over the current time range - .. todo:: Parallelize this (which is trivial) - """ + def run(self): + """Propagate all |Molecule|s in the current |Field| over the current time range""" n_cpu = mp.cpu_count() - print("Running on:",n_cpu,"Processors") + print(f'Running on: {n_cpu} CPUs') self.ensemble.pulse = self.field pool = mp.Pool(n_cpu) results = pool.map(self._propagate, self.ensemble.molecules) @@ -58,13 +64,15 @@ def run(self): self.ensemble.molecules[i] = results[i] return self.ensemble - @classmethod - def molecule(cls, mol, field, dt, timerange): - cls.field = field - cls.dt = dt - cls.t_range = timerange - cls._propagate(cls, mol) - return mol + + # @classmethod + # def molecule(cls, mol, field, timerange, dt_save): + # cls.field = field + # cls.dt_save = dt_save + # cls.t_range = timerange + # cls._propagate(cls, mol) + # return mol + def _propagate(self, molecule): """Propagate an individual molecule in the current field and over the current time range @@ -75,7 +83,7 @@ def _propagate(self, molecule): .. todo:: Should use the modern approach, i.e., `scipy.integrate.solve_ivp` - + # create a copy store a field direction (which can the rotated) field = self.field derivative = np.empty((7,)) @@ -92,19 +100,19 @@ def _propagate(self, molecule): args=(derivative, molecule, field)) """ - integral = ode( self._derivative ) # initial ode object - integral.set_integrator('dopri5', nsteps=10000) # choose the integrator - #integral.set_integrator('lsoda') + # initialize ode object + integral = scipy.integrate.ode(self._derivative) + # choose the integrator + integral.set_integrator('dopri5', nsteps=10000) # alternatively, use integral.set_integrator('lsoda') integral.set_initial_value(np.concatenate((molecule.pos[0].angle.elements, - molecule.pos[0].velocity))).set_f_params(molecule) - while integral.successful() and integral.t <= self.t_range[1]: # start the integration - integral.integrate(integral.t + self.dt) # advance the integrator in time - molecule.pos.append( - Position(quat.Quaternion(integral.y[:4]) - , integral.y[4:7], t=integral.t)) - + molecule.pos[0].velocity)), self.t_range[0]).set_f_params(molecule) + # integrate + while integral.successful() and integral.t <= self.t_range[1]: + integral.integrate(integral.t + self.dt_save) + molecule.pos.append(Position(quat.Quaternion(integral.y[:4]), integral.y[4:7], t=integral.t)) return molecule + def _derivative(self, t, dpos, molecule): """Determine the derivative of a Molecule at a specific time. @@ -126,8 +134,8 @@ def _derivative(self, t, dpos, molecule): q = quat.Quaternion(dpos[0:4]) omega = dpos[4:7] - #E = self.field(t) * self.field.rotate() # uncomment this if you want to rotate the molecule - E = self.field(t) * self.field.rotate(q.inverse) #incase rotating the field + # E = self.field(t) * self.field.rotate() # uncomment this if you want to rotate the molecule + E = self.field(t) * self.field.rotate(q.inverse) # rotating the field position = Position(q, omega) dw = molecule.acceleration(E, position) dq = q.derivative(omega).elements diff --git a/cmiclassirot/sample.py b/cmiclassirot/sample.py index 04ec35c..b28b830 100644 --- a/cmiclassirot/sample.py +++ b/cmiclassirot/sample.py @@ -22,6 +22,8 @@ import scipy.constants import pandas as pd + + class Position(object): """Representation of a phase-space position of a molecule @@ -43,6 +45,8 @@ def __init__(self, angle=quat.Quaternion(), velocity=np.zeros((3,)), t=0.): self.velocity = velocity self.time = t + + class Molecule(object): """Object to be manipulated @@ -61,8 +65,8 @@ def __init__(self, *args, **kwargs): .. code-block:: python - __init__(self, molecule, pos=None, temp=None) - __init__(self, I, P, angle=None, velocity=None, temp=None) + __init__(self, molecule, pos=None, T=None) + __init__(self, I, P, angle=None, velocity=None, T=None, t=0.) :param molecule: Use first version of constructor: copy the specified molecule and possibly update its phase-space position @@ -70,31 +74,37 @@ def __init__(self, *args, **kwargs): :param pos: Determines position of new molecule * 'keep' or None specify to copy the phase-space position from the original Molecule * 'z' places the :class:`Molecule` along the laboratry :math:`z` axis with a angular - velocity of 0 + velocity of 0 * 'random' draws a random direction from a uniform angular distribution and an random - angular velocity according to a 3D Maxwell distribution at temperature `temp` + angular velocity according to a 3D Maxwell distribution at temperature `temp` :param I: principal moments of inertia: these are the three diagonal elements of the in the - inertial tensor in the principle axis of inertia frame (in SI units: kg * m**2) + inertial tensor in the principle axis of inertia frame (in SI units: kg * m**2) :param P: polarizabiity tensor in the inertial frame of the molecule (in SI units: ???) - :param T: Temperature of the thermal distribution from which to draw the initial momentum of - the particle; typically this should not be done here but in a (to be provided) `Source` - :param angle: initial position, as a Quaternion, if not provided a random direction id choosen :param velocity: initial angular frequency, if not provided a random vector drawn from three - one-dimensional Maxell distribution for rotations about x, y, z is drawn + one-dimensional Maxell distribution for rotations about x, y, z is drawn + + :param T: Temperature of the thermal distribution from which to draw the initial momentum of + the particle; typically this should not be done here but in a (to be provided) `Source` + + :param t: Time at which the molecule is created + """ if isinstance(args[0], Molecule): # first variant of constructor -- get args assert(len(args)==1) mol = args[0] + # get keyword args pos = kwargs.get('pos', None) - temp = kwargs.get('temp', 0) + T = kwargs.get('T', 0) + t = kwargs.get('t', 0) + # construct object self._I = mol.I self._P = mol.P @@ -103,12 +113,11 @@ def __init__(self, *args, **kwargs): elif pos == 'z': self.pos = [Position()] elif pos == 'random': - self.pos = [self._thermal_position(temp)] - elif isinstance(pos, Position): + self.pos = [self._thermal_position(T, t)] + elif isinstance(pos, Position(t=t)): self.pos = [pos] elif isinstance(pos, list): self.pos = pos - elif len(args) == 0 or isinstance(args[0], np.ndarray): # second variant of constructor if len(args) >= 1: @@ -121,17 +130,21 @@ def __init__(self, *args, **kwargs): self._P = kwargs.get('P') angle = kwargs.get('angle', None) velocity = kwargs.get('velocity', None) - temp = kwargs.get('temp', None) + T = kwargs.get('T', None) + t = kwargs.get('t', None) # set initila phase-space position self.pos = [] pos = Position() - if temp is not None: - pos = self._thermal_position(temp) + if T is not None: + pos = self._thermal_position(T) # overwrite thermal phase-space position by explicitely specified values if angle is not None: pos.angle = angle if velocity is not None: pos.velocity = velocity + if t is not None: + pass + pos.time = t self.pos.append(pos) else: # something wrong! @@ -151,20 +164,23 @@ def P(self): def acceleration(self, field, position, factor=1): """Calculate the molecule's angular acceleration at its current position in a field - :param field: Field for which to calculate the acceleration + :param field: Field (amplitude, V/m) for which to calculate the acceleration + + + .. todo:: Refactor code to get rid of the IFs. """ q = position.angle v = position.velocity - #I, P = self.rotate(q) # uncomment this line incase you want to rotate + #I, P = self.rotate(q) # uncomment this line incase you want to rotate # the molecule not the field if self._I[2][2]NUL 2>NUL +if errorlevel 9009 ( + echo. + echo.The 'sphinx-build' command was not found. Make sure you have Sphinx + echo.installed, then set the SPHINXBUILD environment variable to point + echo.to the full path of the 'sphinx-build' executable. Alternatively you + echo.may add the Sphinx directory to PATH. + echo. + echo.If you don't have Sphinx installed, grab it from + echo.https://www.sphinx-doc.org/ + exit /b 1 +) + +if "%1" == "" goto help + +%SPHINXBUILD% -M %1 %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O% +goto end + +:help +%SPHINXBUILD% -M help %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O% + +:end +popd diff --git a/doc/use.rst b/doc/use.rst new file mode 100644 index 0000000..0c0fd28 --- /dev/null +++ b/doc/use.rst @@ -0,0 +1,26 @@ +CMIclassirot user manual +======================== + + +Running CMIclassirot +-------------------- + + + +Creating input files +-------------------- + + + +Creating graphical output +------------------------- + + + + +.. comment + Local Variables: + coding: utf-8 + fill-column: 100 + truncate-lines: t + End: diff --git a/examples/CO2-impulsive-alignment.py b/examples/CO2-impulsive-alignment.py new file mode 100755 index 0000000..b75d3ca --- /dev/null +++ b/examples/CO2-impulsive-alignment.py @@ -0,0 +1,66 @@ +# -*- coding: utf-8; fill-column: 120 -*- +# +# This file is part of the CMIclassirot classical-rotation simulations + +"""Example CMIclassirot calculation for the impulsive alignment of warm CO2 + +Laser alignment of CO2 sample at room temperature 296 K with a FWHM = 100 fs Gaussian pulse with an peak intensity of +:math:`$1.5\cdot10^{14}$ W/cm$^2$` + +Compare results to Fig. 1 of [[J. M. Hartmann and C. Boulet, J. Chem. Phys. 136, 184302 +(2012)]](https://dx.doi.org/10.1063/1.4705264). + +@author: Jochen Küpper + +""" + +from math import log, sqrt +import numpy as np +import scipy +from scipy.constants import c, epsilon_0, h, pi, physical_constants +a0 = physical_constants['Bohr radius'][0] + +from cmiclassirot.field import * +from cmiclassirot.sample import * + + +# define timerange of integration and stepsize for saving the distribution +timerange = (-1e-12, 3e-12) # simulate from -1...3 ps +dt_save = 10e-15 # and save positions every 50 fs + + + +# define the alignment field +# We use a Gaussian temporal profile with a FWHM = 100 fs, centered at 0, +# and a peak intensity of 1.5 * 10^14 W/cm**2 +field = Field(peak_intensity=1.5e14 * 1e4, FWHM=100e-15, t_peak=0.) + + + +# define the polarizability tensor of CO2 in atomic units (a_0^3) [Bridge, N. J.; Buckingham, A. D. Proc. R. Soc. A, +# 295, 33 (1966)] +P = np.array([[13.018, 0., 0.], + [ 0., 13.018, 0.], + [ 0., 0., 27.250]]) +# and convert from polarizability volume (a_0^3) to polarizability (C m^2 V^-1) +P *= 4 * pi * epsilon_0 * a0**3 + +# define the inertial tensor of OCS (kg m^2) +# +# OCS rotational constant is B = 0.20286 cm-1 [Herzberg, Electronic spectra and electronic structure of polyatomic +# molecules, (van Nostrand, New York, 1966)] +B = 0.39021 * 1e-2 +I = h / (8 * pi**2 * c) * B +I = np.array([[I, 0, 0], + [0, I, 0], + [0, 0, 0]]) + +# define Molecule based on OCS values +mol = Molecule(I, P, t=timerange[0]) + +# Perform calculation for initial temperature of 298 K (room temeprature) +T = 296. + +# we calculate the alignment for an ensemble of 10,000 randomly sampled molecules +n = 10000 +ensemble = Ensemble(n, mol, T=T, t=timerange[0]) diff --git a/examples/OCS-adiabatic-alignment.py b/examples/OCS-adiabatic-alignment.py new file mode 100755 index 0000000..b85cec1 --- /dev/null +++ b/examples/OCS-adiabatic-alignment.py @@ -0,0 +1,61 @@ +# -*- coding: utf-8; fill-column: 120 -*- +# +# This file is part of the CMIclassirot classical-rotation simulations + +"""Example CMIclassirot calculation: long-pulse alignment of cold OCS with a FWHM = 1 ns Gaussian pulse + +This example uses the rotational constanf of OCS from [Herzberg, Electronic spectra and electronic structure of +polyatomic molecules, (van Nostrand, New York, 1966)] and the polarizability tensor from [Holmegaard et al., Nature +Physics 6, 428 (2010)] + +@author: Jochen Küpper + +""" + +from math import log, sqrt +import numpy as np +import scipy +from scipy.constants import c, epsilon_0, h, pi, physical_constants +a0 = physical_constants['Bohr radius'][0] + +from cmiclassirot.field import * +from cmiclassirot.sample import * + + +# define timerange of integration and stepsize for saving the distribution +timerange = (-2.5e-9, 2.5e-9) # simulate from -2.5...2.5 ns +dt_save = 50e-12 # and save positions every 50 ps + + +# define the alignment field +# We use a Gaussian temporal profile with a FWHM = 1 ns +# and a peak intensity of 10**12 W/cm**2 +field = Field(peak_intensity=1e12 * 1e4, FWHM=1e-9, t_peak=0) + + +# define the polarizability tensor of OCS in atomic units (a_0^3) [Holmegaard et al., Nature Physics 6, 428 (2010)] +P = np.array([[26.15, 0., 0.], + [0., 26.15, 0.], + [0., 0., 50.72]]) +# and convert from polarizability volume (a_0^3) to polarizability (C m^2 V^-1) +P *= 4 * pi * epsilon_0 * a0**3 + +# define the inertial tensor of OCS (kg m^2) +# +# rotational constant is B = 0.20286 cm-1 [Herzberg, Electronic spectra and electronic structure of polyatomic +# molecules, (van Nostrand, New York, 1966)] +B = 0.20286 * 1e-2 +I = h / (8 * pi**2 * c) * B +I = np.array([[I, 0, 0], + [0, I, 0], + [0, 0, 0]]) + +# define Molecule based on OCS values +mol = Molecule(I, P, t=timerange[0]) + +# Perform calculation for initial temperature of 0 K +T = 0. + +# we calculate the alignment for an ensemble of 1000 randomly sampled molecules +n = 1000 +ensemble = Ensemble(n, mol, T=T, t=timerange[0]) diff --git a/examples/OCS-impulsive-alignment.py b/examples/OCS-impulsive-alignment.py new file mode 100755 index 0000000..7c371c7 --- /dev/null +++ b/examples/OCS-impulsive-alignment.py @@ -0,0 +1,61 @@ +# -*- coding: utf-8; fill-column: 120 -*- +# +# This file is part of the CMIclassirot classical-rotation simulations + +"""Example CMIclassirot calculation for the impulsive alignment of cold OCS + +Laser alignment of OCS sample at 2 K with a FWHM = 500 fs Gaussian pulse with an peak intensity of $10^{13}$ W/cm$^2$ + +@author: Jochen Küpper + +""" + +from math import log, sqrt +import numpy as np +import scipy +from scipy.constants import c, epsilon_0, h, pi, physical_constants +a0 = physical_constants['Bohr radius'][0] + +from cmiclassirot.field import * +from cmiclassirot.sample import * + + +# define timerange of integration and stepsize for saving the distribution +timerange = (-10e-12, 40e-12) # simulate from -10...40 ps +dt_save = 250e-15 # and save positions every 250 fs + + + +# define the alignment field +# We use a Gaussian temporal profile with a FWHM = 500 fs, centered at 0, +# and a peak intensity of 10^13 W/cm**2 +field = Field(peak_intensity=1e13 * 1e4, FWHM=500e-15, t_peak=0.) + + + +# define the polarizability tensor of OCS in atomic units (a_0^3) [Holmegaard et al., Nature Physics 6, 428 (2010)] +P = np.array([[26.15, 0., 0.], + [0., 26.15, 0.], + [0., 0., 50.72]]) +# and convert from polarizability volume (a_0^3) to polarizability (C m^2 V^-1) +P *= 4 * pi * epsilon_0 * a0**3 + +# define the inertial tensor of OCS (kg m^2) +# +# OCS rotational constant is B = 0.20286 cm-1 [Herzberg, Electronic spectra and electronic structure of polyatomic +# molecules, (van Nostrand, New York, 1966)] +B = 0.20286 * 1e-2 +I = h / (8 * pi**2 * c) * B +I = np.array([[I, 0, 0], + [0, I, 0], + [0, 0, 0]]) + +# define Molecule based on OCS values +mol = Molecule(I, P, t=timerange[0]) + +# Perform calculation for initial temperature of 2 K +T = 2. + +# we calculate the alignment for an ensemble of 10,000 randomly sampled molecules +n = 10000 +ensemble = Ensemble(n, mol, T=T, t=timerange[0]) diff --git a/examples/co2.py b/examples/co2.py deleted file mode 100644 index a00ccfd..0000000 --- a/examples/co2.py +++ /dev/null @@ -1,29 +0,0 @@ -from cmiclassirot.field import * -from cmiclassirot.sample import * -from scipy.constants import epsilon_0 -from cmiclassirot.propagate import * -import numpy as np -P = 4 * np.pi * 1.e-30 * epsilon_0 * np.array([[ 1.93, 0.E2, 0.E2], - [ 0.E2, 1.93, 0.E2], - [ 0.E2, 0.E2, 4.03]]) -#Polarizabilities are taken from -#https://faculty.missouri.edu/~glaserr/vitpub/jp002927r.pdf - -Iz = 0. -Ix = Iy = 7.14987936e-46 -I = np.array([[Ix ,0.,0.], - [0., Iy,0.], - [0., 0.,Iz]]) -T = 298. -sigma = 0.1e-12/ 2.355 -# define a field -E = Field(amplitude=1.5e18, sigma=sigma, mean=0.25e-12) -n = 1500 -# define a molecule -mol = Molecule(I, P) -# define an ensemble -ensemble = Ensemble(n, mol, temperature=T) -# start simulation -p = Propagate(ensemble, E, timerange=(0.,1.5e-12), dt=1.e-14) -# save results -ensemble._save('CO2_150TWcm2_298K.h5') diff --git a/examples/ocs.py b/examples/ocs.py deleted file mode 100644 index 7d92aa0..0000000 --- a/examples/ocs.py +++ /dev/null @@ -1,21 +0,0 @@ -from cmiclassirot.field import * -from cmiclassirot.sample import * -from scipy.constants import epsilon_0 -from cmiclassirot.propagate import * -import numpy as np - -P = 1.e-30 * epsilon_0 * np.array([[ 3.06, 0.E2, 0.E2], - [ 0.E2, 3.06, 0.E2], - [ 0.E2, 0.E2, 7.46]]) -Ix = Iy = 1.38493163e-45 -I = np.array([[Ix ,0.,0.], - [0., Iy,0.], - [0., 0.,0.]]) -T = 0. -sigma = 1.e-10 / 2.355 -E = Field(amplitude=1.e14/2., sigma=sigma, mean=2.5e-10) -n = 50 -mol = Molecule(I, P) -ensemble = Ensemble(n, mol, temperature=T) -p = Propagate(ensemble, E, timerange=(0.,1.e-9), dt=1.e-11) -ensemble._save('OCS_0.01TWcm2_0K.h5') diff --git a/scripts/cmiclassirot b/scripts/cmiclassirot deleted file mode 100644 index 19f8290..0000000 --- a/scripts/cmiclassirot +++ /dev/null @@ -1,52 +0,0 @@ -#!/usr/bin/env python -# -*- coding: utf-8; fill-column: 120 -*- -# -# This file is part of CMIrotMD -# -# This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public -# License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later -# version. -# -# If you use this programm for scientific work, you must correctly reference it; see LICENSE file for details. -# -# This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied -# warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License along with this program. If not, see -# . -# -# This file shall eventually serve as a driver for the calculation and postprocessing -# for rotMD simulations - - -import h5py -import numpy as np -import matplotlib.pyplot as plt -import sys -fname = sys.argv[1] -rho = 19300. -M = 3.79347234e-18 -L = 50.e-9 -R = 5.e-9 -V = np.pi*L*R**2 -M = V*rho -KB = 1.38064852e-23 -I = np.array([ (0.083*M*L**2 + 0.25*M*R**2), (0.083*M*L**2 + 0.25*M*R**2), 0.5*M*R**2]) -f = h5py.File(fname, 'r') - -time = f['time'] -cos = f['cos'] -E = f['pulse'] -Tt = f['pulseT'] -t = max(np.max(Tt), np.max(time)) -fig, (ax1, ax2) = plt.subplots(2) -ax1.plot(time, cos) -ax2.plot(Tt, E) -ax1.set_xlabel('Time (s)') -ax2.set_xlabel('Time (s)') -ax1.set_ylabel('cos2') -ax2.set_ylabel('I TW/cm2') -ax1.set_xlim(0., t) -ax2.set_xlim(0., t) -plt.tight_layout() -plt.show() diff --git a/setup.py b/setup.py index 411fd79..0871cd1 100644 --- a/setup.py +++ b/setup.py @@ -8,27 +8,32 @@ extra_compile_args = [] library_dirs = [] -long_description = """CMI classical-physics rotational molecular-dynamics simulations +long_description = """CMI classical-physics rotational molecular-dynamics simulations (CMIclassirot) This package provides facilities for the simulation of the rotation dynamics of molecules and particles in external electric fields. """ +name = 'CMIclassirot' +release = '0.9.1' +version = '0.9.1' -setup(name="cmiclassirot", - author = "CFEL Controlled Molecule Imaging group", - maintainer = "CFEL Controlled Molecule Imaging group", - maintainer_email = "jochen.kuepper@cfel.de", - url = "https://www.controlled-molecule-imaging.org/research/further_projects/software", - description = "CMI classical-physics rotational molecular-dynamics simulations", - version = "0.1.dev0", - long_description = long_description, - license = "GPL", - packages = ['cmiclassirot'], - scripts = ['scripts/cmiclassirot'], - python_requires = '>=3.5', - install_requires = ['numpy>=1.16.0', - 'pyquaternion', - 'scipy>=1.3.0'], - ) +setup(name=name, + author = "CFEL Controlled Molecule Imaging group", + maintainer = "CFEL Controlled Molecule Imaging group", + maintainer_email = "jochen.kuepper@cfel.de", + url = "https://www.controlled-molecule-imaging.org/research/further_projects/software", + license = "GPLv3", + description = "CMI classical-physics rotational molecular-dynamics simulations", + long_description = long_description, + version = version, + packages = ['cmiclassirot'], + scripts = ['bin/cmiclassirot', + 'bin/cmiclassirot-plot'], + python_requires = '>=3.9', + install_requires = ['numpy>=1.16.0', + 'pyquaternion', + 'scipy>=1.3.0', + 'tables'], + ) diff --git a/utilities/plot.py b/utilities/plot.py deleted file mode 100644 index 01b037b..0000000 --- a/utilities/plot.py +++ /dev/null @@ -1,32 +0,0 @@ -from cmiclassirot.sample import Ensemble -import matplotlib.pyplot as plt -import numpy as np -import sys - -fname = sys.argv[1] -t=[] -cos2=[] -e = Ensemble.FromFile(fname) -r = np.zeros(len(e.molecules[0].pos)) -E = e.pulse -print("Number of Molecules:", len(e.molecules)) -for i in e.molecules: - t=[] - pulse = [] - for j in range(len(i.pos)): - try: - r[j] = r[j] + i.pos[j].angle.rotation_matrix[2][2]**2 - t.append(i.pos[j].time) - pulse.append(E(i.pos[j].time)) - except: - pass - -r=r/len(e.molecules) -fig, ax = plt.subplots(2) -ax[0].plot(t, r) -ax[1].plot(t, pulse) -plt.xlabel("Time (s)") -ax[0].set_ylabel("Cos^2") -ax[1].set_ylabel("E (V/m)") -plt.tight_layout() -plt.show() diff --git a/utilities/time_of_alignment.py b/utilities/time_of_alignment.py index 74d4242..d8af039 100644 --- a/utilities/time_of_alignment.py +++ b/utilities/time_of_alignment.py @@ -2,31 +2,31 @@ # -*- coding: utf-8; fill-column: 120 -*- import numpy as np from cmiext import const - + # physical constants ε0 = const.vacuum_permittivity c = const.speed_of_light # needed for unit conversion - - + + # definition of an Au cylinder in its principal axis system rho = 19300. # density of gold (kg/m**3) L = 2.e-9 # length of the cylinder (m) R = 0.2e-9 # radius of the cylinder (m) V = np.pi*L*R**2 # volume (m**3) M = V*rho # mass (kg) - + # inertia tensor for a cylinder with the figure axis along z # see ... (link to Wikipedia) I = np.array([(0.083*M*L**2 + 0.25*M*R**2), (0.083*M*L**2 + 0.25*M*R**2), 0.5*M*R**2]) - + ### *** potential issue *** ### ### The following might not be the correct polarizability tensor -- which ### should be diagonal in the principal axis frame, not? ### Or is there something "special" about a *rod*? - + # polarizability tensor for a cylinder with the figure axis along z (m**3) # obtained as polarizability volume from a ZENO calucation (angstrom**3), # see ... (link to blogpost with input and documentation) @@ -36,8 +36,8 @@ # conversion to polariability in SI units (C m**2 / V) # see https://en.wikipedia.org/wiki/Polarizability P *= ε0 / 1e30 - - + + # definition of control field # we assume a constant ("dc") field corresponding to a laser of intensity 10**10 W/cm**2 I0 = 1e14/2. # intensity (W/m**2) @@ -45,16 +45,16 @@ A = np.sqrt((2*I0) / (c * ε0 * 1)) # applied electric field at 45º from Z in XZ plane E = A * np.array((1,0,1))/np.sqrt(2) - + # induced dipole: dot product between of polarizability and electric field dipole = np.dot(P, E) # torque: cross product between the dipole moment and E torque = np.cross(dipole, E) # angular acceleration (with center of mass at origin) a = torque/I - + # simple estimate of the time needed for the rod to rotate Δθ = 1/4 pi in space (Δθ = 1/2 a t**2) t = np.sqrt(np.pi/4 / (np.abs(a)/2)) - + # print the result for the print('Estimated time for 1/4 pi rotation about y (from θ=0..45º in the XZ plane): %.2g s' % t[1])