SCM-NV / pyZacros

Python Library for Automating Zacros Simulations
Other
7 stars 2 forks source link

Labels for Cluster and ElementaryReaction are now automatically set. #14

Closed nfaguirrec closed 4 years ago

lopeztarifa commented 4 years ago

Hi @nfaguirrec

Thanks for this enhancement, it is really useful that the labels are generated automatically. It will actually support a common nomenclature for reactions .. ours ;)

I have two comments:

nfaguirrec commented 4 years ago
  1. This case should be like H+H->H2. The H2 without *. However, I just realize that in order to get an unique name, we have to include the neighboring information as well. That is already included in this branch. Let me know your thoughts about it.
  2. Yes. I'll change it
lopeztarifa commented 4 years ago
  1. This case should be like H+H->H2. The H2 without *. However, I just realize that in order to get a unique name, we have to include the neighboring information as well. That is already included in this branch. Let me know your thoughts about it.

Yes, I got it, gas species should not come with an asterisk.... but my problem is more general. I was wondering how to define a ReactionMechanism where one of the species is in the gas phase. Given the example, H + H -> H2, the object will be something like this:

        myReaction = ElementaryReaction( site_types=( "f", "f" ),
                                         neighboring=[ (1,2) ],
                                         initial=myCluster1,
                                         final=myCluster2,
                                         reversible=False,
                                         pre_expon=1e+13,
                                         pe_ratio=0.676,
                                         activation_energy = 0.2 )

I think there "should not" be any problem with the neighboring [ (1, 2) ] as it is applied/understood only to/for the first cluster, right? But what about the final cluster (myCluster2)? How can it be defined? The underlying problem is that we are mixing clusters with something that isn't (H2).

I guess a workaround would be to declare initial and final clusters as optional arguments, leaving the option that they could be None. Then reading the string of the gas species .... somehow. It is a great idea to construct the mechanism starting from clusters, so it worths the effort to come up with a solution :)

nfaguirrec commented 4 years ago

Oh, Thanks! Now I get it! You are right; for gas phases, we must do something special. But, I think it is not related to a None value for the final cluster. Because the final Cluster is, in fact, a Cluster with all sites empty, but with a particular neighbors structure, however, it may have one or more gas species. I would like to keep the mass balance between the initial and final Clusters. I mean:

  H*   H*                 *    *   H2
  |    |                  |    |
  f -- f    g      -->    f -- f    g

|----initial----|        |-----final-----|

In that way, I feel tempted to include another attribute like "gas_species" in the Cluster class. The following would be an example of the final interface.

s0 = Species( "*" )                    # Empty adsorption site
s1 = Species( "H*", 1 )                # H adsorbed with dentation 1
s2 = Species( "H2", gas_energy=0.0 )   # H2 gas phase

clus1 = Cluster( site_types=( "f", "f" ),
                 neighboring=[ (1,2) ],
                 species=( s1, s1 ),
                 multiplicity=2,
                 cluster_energy=0.1 )

clus2 = Cluster( site_types=( "f", "f" ),
                 neighboring=[ (1,2) ],
                 species=( s0, s0 ),
                 gas_species=( s2 ),
                 multiplicity=2,
                 cluster_energy=0.1 )

reaction = ElementaryReaction( site_types=( "f", "f" ),
                               neighboring=[ (1,2) ],
                               initial=clus1,
                               final=clus2,
                               reversible=True,
                               pre_expon=1e+13,
                               pe_ratio=0.676,
                               activation_energy=0.2 )

However, after this modification, I'm no longer 100% comfortable calling this class "Cluster." What do you think?

The previous code should generate something like this (I hope it is correct):

step H2_desorp
  gas_reacs_prods H2 1
  sites 2
  neighboring 1-2
  initial
    1 H* 1
    2 H* 1
  final
    1 * 1
    2 * 1
  site_types f f
  pre_expon 1.000e+13
  pe_ratio 0.676
  activ_eng 1.000
end_step
lopeztarifa commented 4 years ago

Hi @nfaguirrec,

I agree with the two suggestions you made:

  1. Include this empty-cluster possibility, I think calling it cluster is not fully right, but it's ok. It is an exception we must add if we want to construct the ReactionMechanism from clusters.

  2. About the gas_species attribute, go for it! I think there is no way we can infer it from the initial and final clusters. The user must provide it.

nfaguirrec commented 4 years ago

Done!

Take a look at the following example (available in the ElementaryReaction.test):

        s0 = Species( "*" )      # Empty adsorption site
        s1 = Species( "H*", 1 )  # H adsorbed with dentation 1
        s2 = Species( "H2*", 1 ) # H2 adsorbed with dentation 1
        s3 = Species( "H2", gas_energy=0.0 ) # H2(gas)

        myCluster1 = Cluster( site_types=( "f", "f" ),
                              neighboring=[ (1,2) ],
                              species=( s1, s1 ),
                              multiplicity=2,
                              cluster_energy=0.1 )

        myCluster3 = Cluster( site_types=( "f", "f" ),
                              neighboring=[ (1,2) ],
                              species=[ s0, s0 ],
                              gas_species=[ s3 ],
                              multiplicity=2,
                              cluster_energy=0.1 )

        myReaction2 = ElementaryReaction( site_types=( "f", "f" ),
                                          neighboring=[ (1,2) ],
                                          initial=myCluster1,
                                          final=myCluster3,
                                          reversible=False,
                                          pre_expon=1e+13,
                                          pe_ratio=0.676,
                                          activation_energy = 0.2 )

        print( myReaction2 )

The produced output is:

step H*-f,H*-f:(1,2)-->*-f,*-f:H2,H2:(1,2)
  gas_species H2 2
  sites 2
  neighboring 1-2
  initial
    1 H* 1
    2 H* 1
  final
    1 * 1
    2 * 1
  site_types f f
  pre_expon 1.000000e+13
  pe_ratio 0.676
  activ_eng 0.2
end_step

The code takes into account the stoichiometry of gas-species and if they are either reactives or products following the definition of the Zacros' input files. I mean, if you use gas_species=[ s4,s4 ] in the reactives, it will produce gas_species H2 -2. But if you use it in the products, it will produce gas_species H2 2.

By the way, the label of the cluster follows this format: <adsorbed_species>:<gas_species>:<neighboring>. For example: *-f,*-f:H2:(1,2). Inside each section, items are separated by commas.

On the other hand, I also changed the attribute energy of the Cluster to cluster_energy.

lopeztarifa commented 4 years ago

Great! @nfaguirrec

Thanks, I merge it.