mdanalysis

mdanalysis copied to clipboard

import os

import MDAnalysis as mda

os.chdir("/mnt/e/research/test/poisson_potential/data")
fname = "cg_lj_coul_reduced_lammps"
universe = mda.Universe(f"{fname}.topology", f"{fname}.nc", topology_format="DATA", 
                        atom_style="id type q x y z vx vy vz")

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[35], line 7
      5 os.chdir("/mnt/e/research/test/poisson_potential/data")
      6 fname = "cg_lj_coul_reduced_lammps"
----> 7 universe = mda.Universe(f"{fname}.topology", f"{fname}.nc", topology_format="DATA", 
      8                         atom_style="id type q x y z vx vy vz")

File ~/miniconda3/envs/research/lib/python3.9/site-packages/MDAnalysis/core/universe.py:375, in Universe.__init__(self, topology, all_coordinates, format, topology_format, transformations, guess_bonds, vdwradii, fudge_factor, lower_bound, in_memory, in_memory_step, *coordinates, **kwargs)
    370 coordinates = _resolve_coordinates(self.filename, *coordinates,
    371                                    format=format,
    372                                    all_coordinates=all_coordinates)
    374 if coordinates:
--> 375     self.load_new(coordinates, format=format, in_memory=in_memory,
    376                 in_memory_step=in_memory_step, **kwargs)
    378 if transformations:
    379     if callable(transformations):

File ~/miniconda3/envs/research/lib/python3.9/site-packages/MDAnalysis/core/universe.py:580, in Universe.load_new(self, filename, format, in_memory, in_memory_step, **kwargs)
    577 # supply number of atoms for readers that cannot do it for themselves
    578 kwargs['n_atoms'] = self.atoms.n_atoms
--> 580 self.trajectory = reader(filename, format=format, **kwargs)
    581 if self.trajectory.n_atoms != len(self.atoms):
    582     raise ValueError("The topology and {form} trajectory files don't"
    583                      " have the same number of atoms!\n"
    584                      "Topology number of atoms {top_n_atoms}\n"
   (...)
    588                          fname=filename,
    589                          trj_n_atoms=self.trajectory.n_atoms))

File ~/miniconda3/envs/research/lib/python3.9/site-packages/MDAnalysis/lib/util.py:2553, in store_init_arguments.<locals>.wrapper(self, *args, **kwargs)
   2551             else:
   2552                 self._kwargs[key] = arg
-> 2553 return func(self, *args, **kwargs)

File ~/miniconda3/envs/research/lib/python3.9/site-packages/MDAnalysis/coordinates/TRJ.py:446, in NCDFReader.__init__(self, filename, n_atoms, mmap, **kwargs)
    443 super(NCDFReader, self).__init__(filename, **kwargs)
    445 # ensure maskandscale is off so we don't end up double scaling
--> 446 self.trjfile = NCDFPicklable(self.filename,
    447                              mmap=self._mmap,
    448                              maskandscale=False)
    450 # AMBER NetCDF files should always have a convention
    451 try:

File ~/miniconda3/envs/research/lib/python3.9/site-packages/scipy/io/_netcdf.py:279, in netcdf_file.__init__(self, filename, mode, mmap, version, maskandscale)
    276 self._attributes = {}
    278 if mode in 'ra':
--> 279     self._read()

File ~/miniconda3/envs/research/lib/python3.9/site-packages/scipy/io/_netcdf.py:610, in netcdf_file._read(self)
    608 # Read file headers and set data.
    609 self._read_numrecs()
--> 610 self._read_dim_array()
    611 self._read_gatt_array()
    612 self._read_var_array()

File ~/miniconda3/envs/research/lib/python3.9/site-packages/scipy/io/_netcdf.py:620, in netcdf_file._read_dim_array(self)
    618 header = self.fp.read(4)
    619 if header not in [ZERO, NC_DIMENSION]:
--> 620     raise ValueError("Unexpected header.")
    621 count = self._unpack_int()
    623 for dim in range(count):

ValueError: Unexpected header.

from netCDF4 import Dataset
from scipy.io import netcdf_file

print(Dataset(f"{fname}.nc", "r"))
print(netcdf_file(f"{fname}.nc", "r"))

<class 'netCDF4._netCDF4.Dataset'>
root group (NETCDF3_64BIT_DATA data model, file format NETCDF3):
    Conventions: AMBER
    ConventionVersion: 1.0
    program: LAMMPS
    programVersion: 3 Nov 2022
    dimensions(sizes): frame(1101), atom(1576), cell_spatial(3), cell_angular(3), label(10), spatial(3)
    variables(dimensions): |S1 spatial(spatial), |S1 cell_spatial(spatial), |S1 cell_angular(spatial, label), float32 time(frame), float32 cell_origin(frame, cell_spatial), float32 cell_lengths(frame, cell_spatial), float32 cell_angles(frame, cell_angular), int32 id(frame, atom), int32 type(frame, atom), float32 q(frame, atom), float32 coordinates(frame, atom, spatial), int32 ix(frame, atom), int32 iy(frame, atom), int32 iz(frame, atom), float32 velocities(frame, atom, spatial), float32 forces(frame, atom, spatial)
    groups: 

---------------------------------------------------------------------------
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[39], line 5
      2 from scipy.io import netcdf_file
      4 print(Dataset(f"{fname}.nc", "r"))
----> 5 print(netcdf_file(f"{fname}.nc", "r"))

File ~/miniconda3/envs/research/lib/python3.9/site-packages/scipy/io/_netcdf.py:279, in netcdf_file.__init__(self, filename, mode, mmap, version, maskandscale)
    276 self._attributes = {}
    278 if mode in 'ra':
--> 279     self._read()

File ~/miniconda3/envs/research/lib/python3.9/site-packages/scipy/io/_netcdf.py:610, in netcdf_file._read(self)
    608 # Read file headers and set data.
    609 self._read_numrecs()
--> 610 self._read_dim_array()
    611 self._read_gatt_array()
    612 self._read_var_array()

File ~/miniconda3/envs/research/lib/python3.9/site-packages/scipy/io/_netcdf.py:620, in netcdf_file._read_dim_array(self)
    618 header = self.fp.read(4)
    619 if header not in [ZERO, NC_DIMENSION]:
--> 620     raise ValueError("Unexpected header.")
    621 count = self._unpack_int()
    623 for dim in range(count):

ValueError: Unexpected header.

# ==================================================================== #
# PARAMETERS AND SETTINGS                                              #
# ==================================================================== #

shell           cd /mnt/e/research/test/poisson_potential/data
log             cg_lj_coul_reduced_lammps.log

variable        N equal 1000
variable        T_reduced equal 1
variable        mass_reduced equal 1
variable        epsilon_reduced equal 1
variable        sigma_reduced equal 1
variable        x equal 0.5
variable        N_ion equal v_x*v_N/2
variable        q_reduced equal 13.625994230719545
variable        varepsilon_r equal 78

variable        L_x_reduced equal 8
variable        L_y_reduced equal 7.794228634059947
variable        L_z_reduced equal 20.046884346862004

variable        sigma_ion_wall_reduced equal v_sigma_reduced/2
variable        q_wall_reduced equal -0.19030322569642577
variable        unit_cell_x_reduced equal v_sigma_ion_wall_reduced/2
variable        unit_cell_y_reduced equal v_sigma_ion_wall_reduced*sqrt(3)/2
variable        unit_cell_z_reduced equal v_sigma_ion_wall_reduced*sqrt(6)/3

variable        every equal 1000
variable        timesteps equal 1100000

variable        rng1 equal 12798
variable        rng2 equal 9990
variable        rng3 equal 22263
variable        rng4 equal 101763

# ==================================================================== #
# SYSTEM AND PARTICLES                                                 #
# ==================================================================== #

boundary        p p f
atom_style      charge
region          box block 0 ${L_x_reduced} 0 ${L_y_reduced} 0 ${L_z_reduced} units box
create_box      5 box

mass            * ${mass_reduced}
region          electrolyte block 0 ${L_x_reduced} 0 ${L_y_reduced} ${sigma_ion_wall_reduced} $(v_L_z_reduced-v_sigma_ion_wall_reduced) units box
create_atoms    1 random ${N_ion} ${rng1} electrolyte
create_atoms    2 random ${N_ion} ${rng2} electrolyte
create_atoms    3 random $(v_N-2*v_N_ion) ${rng3} electrolyte
set             type 1 charge ${q_reduced}
set             type 2 charge -${q_reduced}
set             type 3 charge 0
group           electrolyte type 1 2 3

lattice         hcp $((v_unit_cell_z_reduced/0.916486)^-3)
region          left_wall block 0 ${L_x_reduced} 0 ${L_y_reduced} 0 ${unit_cell_z_reduced} units box
create_atoms    4 region left_wall
lattice         hcp $((v_unit_cell_z_reduced/0.916486)^-3) origin 0 0 0.552307414567
region          right_wall block 0 ${L_x_reduced} 0 ${L_y_reduced} $(v_L_z_reduced-v_unit_cell_z_reduced) ${L_z_reduced} units box
create_atoms    5 region right_wall
set             type 4 charge 0
set             type 5 charge 0
group           electrodes type 4 5

# ==================================================================== #
# FORCE FIELD                                                          #
# ==================================================================== #

dielectric      ${varepsilon_r}
kspace_style    pppm 1e-4
kspace_modify   slab 3

pair_style      lj/cut/coul/long $(2^(1/6)*v_sigma_reduced)
pair_coeff      *3 *3 ${epsilon_reduced} ${sigma_reduced}
pair_coeff      4* 4* 0 0
pair_coeff      *3 4* $(100*v_epsilon_reduced) ${sigma_ion_wall_reduced} $(2^(1/6)*v_sigma_ion_wall_reduced)
pair_modify     shift yes

velocity        electrolyte create ${T_reduced} ${rng4}
velocity        electrodes set 0 0 0
fix             nvt electrolyte nvt temp ${T_reduced} ${T_reduced} $(100*dt)
fix             electrodes electrodes setforce 0 0 0

# ==================================================================== #
# ENERGY MINIMIZATION                                                  #
# ==================================================================== #

minimize        1e-6 1e-6 1000 1000
reset_timestep  0
set             type 4 charge ${q_wall_reduced}
set             type 5 charge $(-v_q_wall_reduced)

# ==================================================================== #
# SIMULATION                                                           #
# ==================================================================== #

write_data      cg_lj_coul_reduced_lammps.topology
dump            dump all netcdf ${every} cg_lj_coul_reduced_lammps.nc id type q x y z ix iy iz vx vy vz fx fy fz
thermo          ${every}
thermo_style    custom step spcpu cpuremain temp press ke evdwl ecoul elong
run             ${timesteps}