Hello everyone

I’m still having problem with Buckingham potential for my system.

I used the LJ potential and I see energy transfer from my excited molecule to the bath gas but when I switch to the Buckingham I don’t see the same results. I have to mention that I converted the LJ parameters to the Buckingham potential and then optimized the Buckingham potential so it gives me the same LJ potential up to ~2 kcal/mol at the repulsive wall ( I attached one plot for Ar-H and I did the same for Ar-O) but even with this optimized potential I don’t see the same results.

I used pair-write to write the energy and force for a range of distances for both Ar-H and Ar-O and what LAMMPS calculates and What I calculated are the same. So it seems there is a problem when the trajectories are running.

I looked at each trajectory from Exp6 and LJ individually and I realized that for all LJ trajectories I see relaxation but for Exp-6 it’s not the case. for example if I compare 10 Exp-6 trajectories only two of them show relaxation and the rest show small changes. I also attached two figures for LJ and optimized Exp-6 trajectories (one trajectory) for this section. The initial conditions for both LJ and EXP-6 are the same but it seems when the trajectories are running something wired happens for Exp-6 simulations!

NOW I want to calculate the net forces on each atom in the molecule due to pair interactions and compare these forces between LJ and optimized Exp-6 but I have problem calculating the force on each atom.

1- I have a diatomic molecule(OH) in a bath gas of rare gas. I want to calculate the net force due to pair interactions (Ar-H force and Ar-O force). I used the following:

Compute ID all pair/local force

Compute ID all property/local ptype1 ptype2

.

.

dump 6 all local 1000 dump.pairforce index c_4[1] c_4[2] c_3

but it gives me weird output:

index(integer number) 26682720 26682720 force (real number)

The total index is 2691 and in each line I have the same number for c_4[1] and c_4[2]. I changed the ptype to patom and the same results, only the c_4[] would be different number but again the same large number for both c_4[1] and c_4[2]!

2- If I use “fix ID hydroxyl store/force”, I can calculate to total force on each atom (H and O) and print them via dump and if I use “compute bond/local” I can calculate the force on each atom due to bond interactions and if I subtract this two the results would be the net force on H and O due to pair interactions. right ?

I do appreciate your help.

Best,

Leila

Ar-H.pdf (64 KB)