Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Aspirin example #3

Open
qzhu2017 opened this issue Oct 3, 2023 · 6 comments
Open

Aspirin example #3

qzhu2017 opened this issue Oct 3, 2023 · 6 comments

Comments

@qzhu2017
Copy link
Contributor

qzhu2017 commented Oct 3, 2023

from ost.build import Builder
from pyxtal.db import database
import numpy as np
import os

# Set the crystal model
data = [
        ('ACSALA', [[2,0,0], [0,1,0], [0,0,1]]), # Plastic, 93.99, ac
        ('ACSALA13', [[2,0,1], [0,1,0], [0,0,1]]), # Plastic, 93.99, ac
       ]

style = 'gaff'
db = database('dataset/ht.db')
#dim = [10, 10, 10]

for d in data:
    (code, matrix) = d
    matrix = np.array(matrix)
    print(code)
    xtal = db.get_pyxtal(code)
    smiles = [mol.smile for mol in xtal.molecules]
    bu = Builder(smiles=smiles, style=style)
    bu.set_xtal(xtal, para_min=5.0, dimer=True)
    # Apply the orientation

    # Directory
    folder = code+'-'+style
    if not os.path.exists(folder): os.makedirs(folder)
    cwd = os.getcwd()
    os.chdir(folder)

    # Example 1: tensile, assuming z direction
    task1 = {'type': 'tensile',
             'temperature': 200,
             'pressure': 1.0,
             'max_strain': 0.2,
             'rate': 2e+8, 
           }

    bu.set_slab(bu.xtal, bu.xtal_mol_list, matrix=matrix, dimer=True)
    print('Supercell:  ', bu.ase_slab.get_cell_lengths_and_angles())
    bu.set_task(task1)
    
    bu.lammps_slab.write_lammps()
    bu.ase_slab.write('test.xyz', format='extxyz')
    os.chdir(cwd)
@qzhu2017
Copy link
Contributor Author

qzhu2017 commented Oct 3, 2023

  • compute VACF for dimer centers

@qzhu2017
Copy link
Contributor Author

qzhu2017 commented Oct 3, 2023

@pedroantoniosantosf

Below is the script to compute averaged position and velocities for each molecular center at every 1000 steps.
How can I modify the following script so that I can print out the instantaneous position and velocities for each molecular center at every 5 steps.

compute         cmol all chunk/atom molecule # Define the chunk as a molecule
compute         com all com/chunk cmol       # Calculate the center-of-mass per molecule
compute         vcm all vcm/chunk cmol       # Calculate the velocity of center-of-mass per molecule
fix             3 all ave/time 2 100 1000 c_com[*] c_vcm[*] file center.dat mode vector

@pedroantoniosantosf
Copy link
Collaborator

pedroantoniosantosf commented Oct 4, 2023

@qzhu2017

You can define each valeu as a variable and print the variables in a file:

variable       step         equal step
variable       com_x        equal c_com[1]
variable       com_y        equal c_com[2]
variable       com_z        equal c_com[3]
variable       vcm_x        equal c_vcm[1]
variable       vcm_y        equal c_vcm[2]
variable       vcm_z        equal c_vcm[3]

fix            4 all print 5 "${step} ${com_x} ${com_y} ${com_z} ${vcm_x} ${vcm_y} ${vcm_z}" file center_2.dat screen no title "# step com_x com_y com_z vcm_x vcm_y vcm_z"

@qzhu2017
Copy link
Contributor Author

qzhu2017 commented Oct 4, 2023

@pedroantoniosantosf

ERROR: Variable com_x: Compute global vector in equal-style variable formula (../variable.cpp:1533)

com is vector with two dimensions, [number of molecules, 3]

@qzhu2017
Copy link
Contributor Author

qzhu2017 commented Oct 4, 2023

@pedroantoniosantosf
I think I will just use

fix             3 all ave/time 1 1 1 c_com[*] c_vcm[*] file center.dat mode vector

@pedroantoniosantosf
Copy link
Collaborator

@qzhu2017
Oh, it is true. If you want to print each 5 steps, you can use:

fix             3 all ave/time 1 1 5 c_com[*] c_vcm[*] file center.dat mode vector

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants