Script for calculating static structure factors of molecular liquids

dynasor
1

Hi Everyone,

I am hoping to use dynasor to calculate the static structure factors of molecular liquid using a CHARMM/NAMD DCD file. I want to make sure I’m calling the options correctly. I build a dictionary of the atomic indices by element from a PSF file, but otherwise it reads the coordinates, units, and cell vectors from the DCD. Can someone comment as to whether I’m calling dynasor correctly (see script below)? The profiles are reasonable but I’m not sure the discrepancies are because I’m calling dynasor incorrectly or the underlying potential.

Thanks!

import os
from ase import io
import dynasor
import numpy as np
import json
import MDAnalysis as mda
import re
import sys
from dynasor.post_processing.neutron_scattering_lengths import NeutronScatteringLengths

from dynasor.post_processing import get_weighted_sample
from dynasor.post_processing import get_spherically_averaged_sample_smearing

def extract_atomic_symbols(psf_file):
# Load the PSF file
u = mda.Universe(psf_file)

element_dict = {}

for atom in u.atoms:
    # Extract element symbol using regex (assuming format {element}{number})
    match = re.match(r'([A-Za-z]+)\d*', atom.name)
    if match:
        element = match.group(1)
        if element == 'CL':
            element = 'Cl'
        if element not in element_dict:
            element_dict[element] = []
        element_dict[element].append(atom.index)

return element_dict

molname=sys.argv[1]

Get all .traj files in the current directory

all_atoms =

indices = extract_atomic_symbols(molname + ‘.psf’)

trajectory_format=‘dcd’
traj=dynasor.Trajectory(molname + ‘.’ + trajectory_format,
trajectory_format,
atomic_indices=indices,
frame_start=0,
frame_stop=None,
frame_step=1)

n_atoms=len(indices)
frame=traj.next()
q_max = 10.0
q_points = dynasor.get_spherical_qpoints(traj.cell, q_max=q_max)
sample = dynasor.compute_static_structure_factors(traj, q_points)

q_linspace = np.linspace(0, q_max, 101)
q_width = 0.01

compute Sq

sample_averaged = get_spherically_averaged_sample_smearing(sample, q_norms=q_linspace, q_width=q_width)
print(sample_averaged.atom_types)
nsl = NeutronScatteringLengths(sample_averaged.atom_types)

sq=sample.Sq / n_atoms
sample_weighted = get_weighted_sample(sample_averaged, nsl)

fhout=open(molname + ‘_neutron_averaged.csv’, ‘w’)

fhout.write(‘q,sq\n’)

for i in range(0, len(sample_weighted.q_norms)):
fhout.write( str(sample_weighted.q_norms[i]) + ‘,’ + str(sample_weighted.Sq[i][0]) + ‘\n’)

2

Hi Christopher!

I think your script looks reasonable. Note that there are a few parameters that affect the convergence of your result. These include the length of your trajectory (#number of frames), the step between read frames (frame_step), the number of q-points, and the smearing factor q_width. In particular, if you have MD system of “regular” size, only using 101 q-points up to 10.0 1/Angstrom would yield a very sparse sampling in q.

By the way, why do you call frame = traj.next()? I don’t see that being used anywhere.

Here is a link to a tutorial which nicely demonstrates the effect of these different convergence parameters when computing the static structucture factor: Static structure factor in halide perovskite (CsPbI3) — dynasor documentation