Hi lammps-users,

sorry for the delay, and thanks for the reply of steve. I have attached the paper of GULP, which contains how GULP coompute electrostatic potential. i am sorry i don’t really understands how ewald summation works. Actually in that simple test, i did’t exclude the electrostatic potential and L-J interaction between 1-2, 1-3, 1-4, since that is how the model is designed. i attach my input script for lammps.

regard to Andrew’s help, thanks very much, i downloaded and used the attached script, and it works exactly what i wish, instead generating totally 12 impropers for single tetrahedral, now is gives me 4 impropers for each tetrahedral. and for the reason why i am using class_2 improper for SiO4 system, that because, in the force field that i am trying to implement, there is a crossangle interaction term, which couples the deformation of two O-Si-O angle on the same side. (see the following picture), i found in lammps doc, that the definition of the second term in class_2 improper potential couples the deformation of the angles, so i use class_2 improper.

And actually, the original valence force field contains not only the coupling of the deformation of angles on the same side, but also the coupling of deformation of angles on the diagonally (angles on different sides), and the coupling of the bond length change with angle deformation, where the bond are not included in the angle, but i didn’t find how to implement these two coupling terms in lammps doc. i will be very appreciative if anyone knows hoe to do that in lammps.

Best Regards

Jiasen Guo

in.MZHB_SOD_minization (1.06 KB)

MolSimulGULP.pdf (750 KB)