Automatic setting of LJ pair coefficients is not evaluated by lj.pair_coeff.set

[Nando V]

I am trying to include in the python code a loop to automatically set the lennard-jones coefficients for a system with many types of atoms.

The mixing rule is by Lorentz-Berthelot.

I have a list with all the epsilon, another with all the sigmas, and one more with the r_cut (to set as 0 if there is no force).

Usually, I set as:

  lj.pair_coeff.set('O', 'Ow', epsilon=sqrt(O_LJ[0]*Ow_LJ[0])/28.0, sigma=(O_LJ[1]+Ow_LJ[1])/2, r_cut=14.0)

which needs to be repeated in many lines for a few atoms.

Since this is very prone to error, I tried to make a loop.

And the loop I made is:

for i in range(0,len(snapshot.particles.types)-2):

  for j in range(0,len(snapshot.particles.types)-1-i):

    print(

      snapshot.particles.types[i],

      snapshot.particles.types[j],

      sqrt(epsilon_particles[i]*epsilon_particles[j]),

      (sigma_particles[i]+sigma_particles[j])/2,

      (rcut_particles[i] and rcut_particles[j])

      )

    lj.pair_coeff.set(      

      snapshot.particles.types[i],       

      snapshot.particles.types[j],      

      epsilon=sqrt(epsilon_particles[i]*epsilon_particles[j]),

      sigma=(sigma_particles[i]+sigma_particles[j])/2,      

      r_cut=(rcut_particles[i] and rcut_particles[j])

      )

with the vectors defined as:

epsilon_particles=[28.0, 85.55, 15.1, 0.10, 89.5736, 0, 0]

sigma_particles=[3.36, 3.07, 2.42, 2.0, 3.097, 1 , 1]

rcut_particles=[14.0, 14.0, 14.0, 14.0, 14.0, 0.0 , 0.0]

after running, the print part (blue) prints the correct values in the screen

However, in the pair coefficient setting (red), none of the values are evaluated (hoomd just takes it as text).

Then, before running, hoomd complains that the type pairs are not found.

Do anyone know which part I'm missing??

[wsxu...@gmail.com]

Did you set the coefficients for all possible pairs?  I simulated a system with many particle types, and the following scheme worked for me.

PType = system.particles.types

NType = len(PType)

for i in range(NType):

for j in range(i, NType):

DIJ = 0.5*0.01*(float(PType[i]) + float(PType[j]))

lj.pair_coeff.set(PType[i], PType[j], epsilon=1.0, sigma=DIJ, r_cut=2.0**(1.0/6.0)*DIJ)

In this case, the diameters of particle types are used for their names. One can certainly define a separate list for the diameter or epsilon values.