add Xygauss surface

[lmiq]

@sschott

I merged the changes here, and now you can clone the repo and checkout this branch.

The package compiles. Can you give some simple example of the usage of this constraint, so I can add to the tests?

I have tested this:

filetype pdb
output xzgauss.pdb

tolerance 2.0

structure CLA.pdb
    number 50
    inside box -50 -50 -50 50 50 50
    over xygauss 21.0 5.0 10.0 20.0 -23.0 15.0
end structure

where CLA.pdb is:

HEADER
HETATM    1 CLA  CLA    2        0.000   0.000   0.000

but I got this error:

Program received signal SIGFPE: Floating-point exception - erroneous arithmetic operation.

Backtrace for this error:
#0  0x7f2d132c1ad0 in ???
#1  0x7f2d132c0c35 in ???
#2  0x7f2d12fd151f in ???
        at ./signal/../sysdeps/unix/sysv/linux/x86_64/libc_sigaction.c:0
#3  0x7f2d13227a00 in xflow
        at ../sysdeps/ieee754/dbl-64/math_err.c:40
#4  0x7f2d131eef02 in __GI___exp
        at ./w_exp_template.c:32
#5  0x559bef155cbd in comprest_
        at src/comprest.f90:170
#6  0x559bef17703d in computef_
        at src/computef.f90:78
#7  0x559bef0df8c4 in initial_
        at src/initial.f90:211
#8  0x559bef168e61 in packmol
        at app/packmol.f90:686
#9  0x559bef16bd6a in main
        at app/packmol.f90:36
Floating point exception (core dumped)

Is there some parameter missing in the example of the xygauss surface? Or do you see any other obvious problem?

[SSchott]

I don't see any particular problem, but noticed that after the merge, I get Note: The following floating-point exceptions are signalling: IEEE_UNDERFLOW_FLAG, which seems to be a fairly common numerical warning (at least in the Amber world :grimacing: ). Anything changed in particular in the initial molecule placement? Your example works if you compile without the error traps (so just make, not make devel)

I will see if I can trace where this originates.

[lmiq]

It is probably some error in the constraint parameters. I may have missed some change. It is better to try with the minimal example possible.

[SSchott]

ok, I backtraced it to the exponent function; when the argument is about <= -500, which essentially makes it 0, the error gets triggered as division by zero, so it is not really a problem, but I get that it is not elegant as it is.

As I am not really a Fortran coder, maybe you can give me a hint here. Is there any way of dealing with this? I tested placing a "if argument <= -500" check, and that works without errors, but I guess that when the error gets triggered might depend on the precision of the machine where this is compiled. If such an if statement would do it, let me know and I can pass a commit.

[lmiq]

The best is to compute the denominator first and, if it is smaller than, say, 1e-10, just set the result to a large number. The problem is that the function can go to infinity, which will be always bad in terms of modeling with this function.

It may be reasonable to throw an informative error the first time such values are encountered, to indicate to the user that the extent of the constraint is not reasonable anyways.

[lmiq]

I think that this is overall fine. I won´t be able to merge it and test it carefully for the following two weeks. It would be to have some examples of the surfaces being produced with the correct geometry. With the parameters I've tested I could not get anything that looked differently from a plane. @SSchott

[SSchott]

If I use the very test you did up here, but increasing the number of CLA to 10000, then, at least on my side, the curve is noticeable.

filetype pdb
output xzgauss.pdb

tolerance 2.0

structure CLA.pdb
    number 10000
    inside box -50 -50 -50 50 50 50
    below xygauss 21.0 5.0 10.0 20.0 -23.0 15.0
end structure

Here a pic of the result. Can you confirm you get something similar?

[lmiq]

Sorry for the long delay in moving this forward.

Actually I'm still having a hard time getting that result. This is the input file:

filetype pdb
output xyzgauss.pdb

tolerance 2.0

structure CLA.pdb
    number 10000
    inside box -50 -50 -50 50 50 50
    below xygauss 21.0 5.0 10.0 20.0 -23.0 15.0
end structure

And I get:

There is some distortion, but I cannot see that there is a gaussian surface there. It would be nice to have a more explicit example, just to be sure.

[SSchott]

Ok, this is a head scratcher for me, but I somehow managed to get your behavior. In my case, the CLA input files give the expected result (xygauss), so I played around with a "low" and "high" constrained water packing with a "tougher" gaussian shape. When there is enough spacing in the initial approximation to make it flat, then you don't get the curved surface. On the other hand, if you force some clashes by increasing the tolerance, you geat a curved surface as expected. Interestingly, if you open the intermediate pdb files, the "Gaussian peak" gets "populated" "bottom up". This doesn't happen if you have more complicated packing scenarios that disrupt a perfect cuboid, or if you use the inverted shape, so, "above" instead of "below", despite having comparable densities.

You can find all of the examples in the attached tar file, with the output pdbs in the saved subfolder. I am not sure how the initial guess is constructed, but could that have something to do, or do you have another idea where the initial bias could be coming from? to_send.tar.gz

[lmiq]

The input file is a disperse distribution of the particles, more or less homogeneous in the space defined by maximum and minimum coordinates that are determined from the constraints. Thus, in initial packing steps, many particles will be moved rapidly to the boundaries of the constraints, in such a way that the constraints may become more apparent. Packmol only requires that the constraints are fulfilled, not that they are fulfilled in any particular way, so if a flat surface is enough, it can be a solution.

Since you mentioned that packing density may be an issue, then everything seems fine. With 10k particles the packing was too easy, and there was room enough for the surface not being visible. With 30k particles it is clear:

I think everything is correct. I will update the package and the user guide accordingly. Thank you for your input and patience. I am sorry for letting this PR be stalled for so long.

[SSchott]

Great, nice to see that it worked :+1:

[lmiq]

Just to mention that it is now released in v20.13.0, and documented in the new page (https://m3g.github.io/packmol).