Move cvm::atom member data to synchronized vectors within cvm::atom_group

[giacomofiorin]

This issue can be broken up in multiple steps:

• Transform the atom_group class into a container for vectors of shape xxxx..., yyyy... and zzzz... for positions, velocities and total forces.
• Change the atom class into an iterator of the atom_group vectors, which implements for each property: (a) a property() accessor returning a copy of the corresponding property as a 3-vector, (b) a set_property() function that takes a 3-vector and writes it to corresponding location in the atom_group vector, and (c) an incr_property() function.
• Add pointers inside atom_group to vectors containing the gradients and fit gradients of each CVC and of each of its components (e.g. x through z for vector-valued CVCs). For CVCs that are affected by the fit, gradients and fit gradients should be requested at the same time to the atom group. For those that don't.
• [ ] Add laboratory-frame versions of the above if useful. (removed because no clear use case)
• Implement exceptions to the above for those CVCs like distancePairs where the number of components is prohibitive.
• Collect pointers to all of the above (positions, total forces, gradients, fit gradients) into a single class (nested in the atom_group namespace), together with a unique integer allowing a CVC to search the corresponding object in the atom group. This should allow quick access and garbage collection when a CVC is deleted but not the atom group.

[jhenin]

Another route would be moving gradients to the cvc objects - at least that seems more intuitive.

[giacomofiorin]

Yeah, but then the cvc objects have to deal with changes in frame of reference, fit gradients, and such. With no solution giving a perfect division of labor, I prefer having a lightweight CVC class and a heavy atom group class, instead of two classes that are both heavy.

[jhenin]

This means that at some point, when a CVC applies a force to an atom group, the atom group will have to look up among its stored gradients which ones are those of that particular CVC. My impression is that it would be simpler for the CVC to pass a reference to its own gradients, together with the force.

One way or another, a reference to the gradients will be needed, but I suppose it doesn't make a difference which object they are actually stored in... Except, say, when the CVC is deleted, in which case its gradients should be deleted as well. If the gradients were stored in the CVC, the atom group wouldn't need to know that a CVC has come and gone. It would just receive forces, or not. So to me the current setup seems like the atom group is doing someone else's work for them. Of course that is natural in the current situation where the atom group belongs to the CVC anyway, but not if that changes.

[giacomofiorin]

The atom group can easily store pointers to the arrays of gradients, and set a pointers to NULL when the corresponding CVC are deleted. To store atomic data (positions vs gradients) in two different classes sounds to me more complex, actually.

…

On Tue, Mar 14, 2017 at 4:25 PM, Jérôme Hénin ***@***.***> wrote: This means that at some point, when a CVC applies a force to an atom group, the atom group will have to look up among its stored gradients which ones are those of that particular CVC. My impression is that it would be simpler for the CVC to pass a reference to its own gradients, together with the force. One way or another, a reference to the gradients will be needed, but I suppose it doesn't make a difference which object they are actually stored in... Except, say, when the CVC is deleted, in which case its gradients should be deleted as well. If the gradients were stored in the CVC, the atom group wouldn't need to know that a CVC has come and gone. It would just receive forces, or not. So to me the current setup seems like the atom group is doing someone else's work for them. Of course that is natural in the current situation where the atom group belongs to the CVC anyway, but not if that changes. — You are receiving this because you were assigned. Reply to this email directly, view it on GitHub <#61 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AETZWWg0QClN_VRZ_ROCO_2RZbHmI68Zks5rlvefgaJpZM4IxYGN> .

-- Giacomo Fiorin Associate Professor of Research, Temple University, Philadelphia, PA Contractor, National Institutes of Health, Bethesda, MD http://goo.gl/Q3TBQU https://github.com/giacomofiorin

[HanatoK]

Yeah, but then the cvc objects have to deal with changes in frame of reference, fit gradients, and such. With no solution giving a perfect division of labor, I prefer having a lightweight CVC class and a heavy atom group class, instead of two classes that are both heavy.

Here's my thoughts. I think if the nesting CVs are allowed, then the gradients can be put into CVC class while still being lightweight. Let's say if there is a class derived from CVC named AlignedAtoms, then you can put all fitting algorithm into it. Apart from the values and gradients of atomic coordinates, AlignedAtoms also provides interfaces to query atomic indexes and names. It just like a moving frame of virtual atoms. Since CVCs can be nested, so in the end we can have some configuration like:

AlignedAtoms {
    name        aligned_reference
    indexGroup  reference
    atoms {
        centerReference    on
        rotateReference    on
        fittingGroup {
            indexGroup  protein
        }
        refpositionsfile  ./complex_largeBox.xyz
    }
}

AlignedAtoms {
    name        aligned_ligand
    indexGroup  ligand
    atoms {
        centerReference    on
        rotateReference    on
        fittingGroup {
            indexGroup  protein
        }
        refpositionsfile  ./complex_largeBox.xyz
    }
}

colvar {
    name polarPhi
    customFunction atan2(i3, i1) * 180 / 3.1415926
    period  360
    wrapAround 0.0
    distanceDir {
        name  i
        group1 {
            virtualAtoms aligned_reference
        }
        group2 {
            virtualAtoms aligned_ligand
        }
    }
}

In this way, you can move the gradients into CVCs (such as polarPhi above), and don't need to handle the fitting gradients in it. The fitting gradients can be handled by AlignedAtoms. When applying forces the chain rule should be OK.

[jhenin]

A related point: I've noted that in practice, fitting is nearly often done per-CVC and not per atom group. I can't think of an example where different atom groups are fitted in a different way within a CVC, and intuitively it doesn't make much sense to me.
So maybe the task could be offloaded to the CVC, avoiding redundant calculations.

[giacomofiorin]

@HanatoK Whichever way you look at it, gradients and fit gradients are specific to the combination of CVC and atom group, so storage for those will be shared by both, requiring changes to the existing layout of both classes. It would be good not to reinvent the wheel, so at the moment I'm trying to see whether we can finally use shared_ptr or something similar, and the CVC gradient computation is not something that one could #ifdef out of VMD.

@jhenin Again along those lines, a fit is determined by the CVC-group combination, so storage in either class would make sense to me. But the most immediate case of redundant calculations that you could think of are the path CVs, where you fit to a different reference structure each time and so the fit is different for each: I don't see how can you save work by moving the fitting around, as most of the terms are not redundant after all. Am I missing something?

[jhenin]

I agree that gradients and fit gradients are specific to the combination of CVC and atom group. But while an atom group could potentially be shared and re-used between different CVCs, a single CVC cannot be "shared/re-used between different atom groups". So to me the only unbreakable link is between gradient and CVC.

About redundant fits, I meant fitting several atom groups independently within the same CVC. The most common case of this happens with distanceVec and distanceZ.

[giacomofiorin]

I think it's more than that. With distanceVec or distanceZ it's not just the fitting part (which isn't even there all the time) but the entire function that is computed redundantly. Hence #232 no?

[jhenin]

There are cases for distanceZ when several projections are needed, yes - that would be good to address as well. For distanceVec, I can't think of a case. Here I was only addressing the case of fitting both groups the same way.

PS: if this is going to be an architecture change, it might as well address all these issues at once if possible.

[giacomofiorin]

I can't think of a way that an architecture change, no matter how limited, wouldn't be involved at some point. We're discussing redundancies between different groups of the same CVC, different CVCs of the same colvar and different colvars altogether. That's three levels of hierarchy. Issues like this one don't stay open for 5+ years based on procrastination alone.

So any discussion is of course useful/welcome to see all sides of it. More automated (unit) tests would also help, even if those cover code that is very stable by now, but would need to be changed for performance. Admittedly that's less visible work, but no less important.

[jhenin]

Re-reading this conversation in light of #784 , I am even more convinced by @HanatoK 's point above that the fitting operation can be expressed as a CVC. In a sense it would generalize the cartesian CVC that just forwards unmodified atom coordinates.
The only thing that is missing from this formalism is factorizing the alignment process for the fitting group, but if atom groups could be made reusable separately, we'd have all the pieces.

[HanatoK]

This issue can be broken up in multiple steps:

• Transform the atom_group class into a container for vectors of shape xxxx..., yyyy... and zzzz... for positions, velocities and total forces.

This should be completed in #788.

• Change the atom class into an iterator of the atom_group vectors, which implements for each property: (a) a property() accessor returning a copy of the corresponding property as a 3-vector, (b) a set_property() function that takes a 3-vector and writes it to corresponding location in the atom_group vector, and (c) an incr_property() function.

#788 does implement getter and setter for xyz components separately, but it does not provide incr_property(). I am not sure what incr_property() should do.

• Add pointers inside atom_group to vectors containing the gradients and fit gradients of each CVC and of each of its components (e.g. x through z for vector-valued CVCs). For CVCs that are affected by the fit, gradients and fit gradients should be requested at the same time to the atom group. For those that don't.

#788 uses grad_x(i), grad_y(i) and grad_z(i) (the same for fit gradients). Maybe this is not the best design.

• Add laboratory-frame versions of the above if useful.

What are the "laboratory-frame versions"?

• Implement exceptions to the above for those CVCs like distancePairs where the number of components is prohibitive.

I don't understand what exceptions are here. Should I throw exceptions in distancePairs? Why is the number of components prohibitive?

• Collect pointers to all of the above (positions, total forces, gradients, fit gradients) into a single class (nested in the atom_group namespace), together with a unique integer allowing a CVC to search the corresponding object in the atom group. This should allow quick access and garbage collection when a CVC is deleted but not the atom group.

Unfortunately this is not implemented in #788. When would the CVC search the atom group?

These steps were proposed before I started my PhD, and I knew nothing about Colvars before June 2016. There might be something that I overlooked.

[giacomofiorin]

This should be completed in #788.

:-)

#788 does implement getter and setter for xyz components separately, but it does not provide incr_property(). I am not sure what incr_property() should do.

That's for incrementing gradients of forces with nested loops, but a different design can work as well.

#788 uses grad_x(i), grad_y(i) and grad_z(i) (the same for fit gradients). Maybe this is not the best design.

However, it exists, so it is de facto the best design :-)

What are the "laboratory-frame versions"?

Those that are not in the rotated frame.

I don't understand what exceptions are here. Should I throw exceptions in distancePairs? Why is the number of components prohibitive?

The language implied that we could have explicit gradients for vector variables as well, provided that these are of low enough dimensions (i.e. 3 or 4). In #713 you worked around that issue by not requiring the explicit gradients at all.

Unfortunately this is not implemented in #788. When would the CVC search the atom group?

That would be useful for sharing atom groups between multiple CVCs, which is a separate improvement from the SOA layout. As long as this can be implemented separately, it's all good. We do not need to check all the boxes at the same time.

These steps were proposed before I started my PhD, and I knew nothing about Colvars before June 2016. There might be something that I overlooked.

I don't think so, consider also that this issue/discussion hasn't been kept in sync with code changes in the intervening years.

My more pressing concern is that in #788 we have two coexisting code paths. If we agree that the SOA code path that you added works, do you see any reason for keeping the old AOS one?

[HanatoK]

The new SOA code path should cover all the functionalities that the AOS code path has, except that the insertion, delete and moving of atoms are not that straightforward, since I used a single buffer ("straight array") for all three xyz components of atomic coordinates.

I kept the AOS code path and just thought that keeping it might be more acceptable to you, and make the PR easier to approve. Do you think that I need to remove the AOS code path entirely?

[giacomofiorin]

I'm super-happy with removing old code, as long as we are not (1) losing features that already worked or (2) making it more difficult to finish others that were in the making.

It is with that intent that I opened #789, which was in good shape before NAMD 3, when it became moot, and now maybe useful again. On a quick look, #788 and #789 do not seem to conflict too much, and #788 is more universally useful anyway so I'd choose to have that, if forced to. But I would like to hear your thoughts on how these pieces fit together with what you had in mind.

[HanatoK]

#788 has been merged, so the first step should be completed. There are alternative for other steps. Shall we close this issue and discuss the GPU work in #652?

[giacomofiorin]

I mean, it's not required to be necessarily up to date like an Issue would

[jhenin]

We're beginning to accumulate many open discussions, so I'm in favor of closing them once they become stale.

[giacomofiorin]

We're beginning to accumulate many open discussions, so I'm in favor of closing them once they become stale.

Well, in that case better close them, then. @HanatoK can you start a new one?

[HanatoK]

Well, I think the other steps may not be necessary, or perhaps I still don't quite understand them and your original goals. Now cvm::atom has been changed to cvm::atom_group::simple_atom and there is no cvm::atom anymore. cvm::atom_group::simple_atom (

colvars/src/colvaratoms.h

Lines 106 to 115 in d06a030

	struct simple_atom {
	int proxy_index;
	int id;
	cvm::real mass;
	cvm::real charge;
	cvm::atom_pos pos;
	cvm::rvector vel;
	cvm::rvector total_force;
	cvm::rvector grad;
	};

) is intended to be POD and may be used in the GPU code path in the future. The atom_modifier is responsible for synchronizing the changes of "AOS atom groups" to the SOA data structure, so I think the goal indicated in the discussion title has been completed.

What should be included if I start a new one? Could you clarify and elaborate your longterm roadmap for improving the atom group data structure?

[giacomofiorin]

I was referring to the discussion of GPU porting work. You had suggested to continue that in #652, which is also an old discussion that is not as easy to navigate. I think it would be clearer if you started a new one.

As far as this discussion goes, I agree that several of the points are outdated, and more importantly the main goal is taken care of by #788

[HanatoK]

I continue the discussion of GPU related tasks in #816.