Control the temperature only for monomers

Dear LAMMPS Community,

I’m trying to model vapor cooling, more specifically nonisothermal nucleation, so I need to cool down only monomers (alone atoms).
My plan was so: do cluster analysis every Nts timesteps, get a set of monomer IDs, define a group to which only monomers belong, and set nvt fix for this group.
So, I came up with the following

variable totalRuns equal 100
variable deltaTemp equal $((v_originTemp-v_targetTemp)/v_totalRuns)
compute ClstAnlysis all cluster/atom 1.5
dump ADIOS_DUMP all custom/adios ${Nts} dumps/START0 id c_ClstAnlysis

variable runN loop ${totalRuns}


    # **magic setup of variable is_monomer for every atom**
    group monomers variable is_monomer
    variable originTempX equal $(v_originTemp - v_deltaTemp*v_runN)
    variable targetTempX equal $(v_targetTemp + v_deltaTemp*(v_totalRuns-v_runN-1))
    fix NVTFIX monomers nvt temp ${originTempX} ${targetTempX}  $(1000.0*dt)
    run ${Nts}

    next runN

    if "${runN} <= ${totalRuns}" then "jump SELF MAIN_LOOP"

unfix NVTFIX

As I understood it, the variable ‘is_monomer’ must be a boolean and True, if an atom is a monomer.
I’ve thought it will be possible to do it with python, i.e., in place of magic stuff, write python getMonomers invoke, and before loop, write

python getMonomers input 1 c_ClstAnlysis return v_is_monomer format ???? file

but recognized that in LAMMPS it is impossible to pass/receive arrays to/from a python script.

I think it is possible to set up is_monomer variable using LAMMPS loops, but how? How to predefine an empty atom-style variable and set its members like “v_is_monomer[i]=True”?

Or is there a better way to do it?

Thank you in advance for your assistance.

P.S.: I’ve even started thinking about writing a python function getMonomers that is a client that sends requests to a server, and a server in turn reads dumps/START0 dump, evaluates last c_ClstAnlysis, writes the result into a file, and then LAMMPS reads this file to set variable is_monomer. It will be really slow, I think.

I would use the compute coord/atom command — LAMMPS documentation.
This will return the number of neighboring (non-bonded) atoms within the cutoff as a per-atom vector. Thus for “monomers” the value would be 0.

So something like the following (untested) would give you what you are looking for.

compute neighs all coord/atom cutoff 1.5
variable is_monomer atom c_neighs==0
group monomers variable is_monomer

Please keep in mind that the group command is incremental, so if you are using it in a loop, you first need to delete the group, and to avoid a crash you should first create an empty group of the same name before the loop.

There is no way to do this from python since python style variables can only return a single global value.

1 Like

Thank you!

Should I unfix a fix before doing a new fix with the same name?

This is documented behavior, so please see the manual. What you should be doing depends on what you want to achieve: just update the settings of the fix or re-initialize the fix from scratch? Whether there is a difference between those two options depends on the individual fix.

1 Like