akohlmey / topotools

VMD plugin for manipulating topology information
Other
28 stars 16 forks source link

Issues with Improper Atom Ordering in TopoTools 1.9 for LAMMPS #25

Open allotrust opened 11 months ago

allotrust commented 11 months ago

Hello,

I have been using TopoTools version 1.9 and encountered an issue when using the "topo guessimpropers" command. Specifically, the central atom of the improper (the one with 3 bonds) is not always listed first in the impropers list. This is problematic as LAMMPS requires the central atom to be the first one listed in the improper angle.

To investigate this behavior, I modified the 'topoimpropers.tcl' script. To simplify its logic I removed the sorting parts and atom type parts, which don't affect the inconsistencies I was observing.

Here's a simplified snippet of my code:

foreach bonds $bonddata aidx $atomindex atyp $atomtypes { set nbnd [llength $bonds] if {$nbnd == 3} { lassign $bonds b1 b2 b3 set type [join [list $atyp X X X] "-"] lappend newimproperlist [list $type $aidx $b1 $b2 $b3] } }

$aidx variable represent the first atom in the generated improper (normally) but unfortunetly it does not always appear first in the impropers list when the line lappend newimproperlist [list $type $aidx $b1 $b2 $b3] is used. However, if I alter the line to lappend newimproperlist [list $type $aidx 0 0 0], the $aidx index is consistently placed first in the generated impropers. As soon as I reintroduce $b1 $b2 $b3 into the line, the ordering inconsistency with $aidx resumes.

I have been experimenting with this a lot and have yet to understand the root of this behavior. Any insight would be greatly appreciated.

Best regards, Sam

allotrust commented 11 months ago

I've identified the issue. It doesn't originate from TopoTools.

In the TopoTools.tcl script, there's a line that utilizes the constructed list of impropers to set impropers in the structure using the following command: molinfo $mol set impropers [list $newimproperlist]

When you print newimproperlist just before this line, and the list of impropers after it, you can observe and compare how the improper list was manipulated in the process. Here's the code block:

puts "newimproperlist: $newimproperlist" molinfo $mol set impropers [list $newimproperlist] puts [molinfo $mol get impropers]

The output I receive is: newimproperlist: {CW1-X-X-X 0 56 1 2} {CS-X-X-X 2 0 5 4} {CW1-X-X-X 3 1 6 4} {CS-X-X-X 4 3 7 2} {CW1-X-X-X 6 9 10 3} {CS-X-X-X 9 15 6 14} {CS-X-X-X 14 19 16 9} {CW1-X-X-X 16 14 10 26} {CW1-X-X-X 27 29 28 97} {CS-X-X-X 29 32 31 27} {CW1-X-X-X 30 31 28 33} {CS-X-X-X 31 30 29 34} {CW1-X-X-X 33 37 36 30} {CS-X-X-X 36 33 42 41} {CS-X-X-X 41 46 36 43} {CW1-X-X-X 43 37 48 41} {CW2-X-X-X 48 53 52 43} {CS-X-X-X 53 58 48 57} {CW2-X-X-X 56 0 57 52} {CS-X-X-X 57 53 56 59} {CW1-X-X-X 68 67 70 69} {CS-X-X-X 70 68 72 73} {CW1-X-X-X 71 69 72 74} {CS-X-X-X 72 71 75 70} {CW1-X-X-X 74 77 78 71} {CS-X-X-X 77 83 82 74} {CS-X-X-X 82 87 84 77} {CW1-X-X-X 84 82 89 78} {CW2-X-X-X 89 84 94 93} {CS-X-X-X 94 98 89 99} {CW2-X-X-X 97 93 98 27} {CS-X-X-X 98 94 97 100}

{{CW1-X-X-X 2 1 56 0} {CS-X-X-X 2 0 5 4} {CW1-X-X-X 3 1 6 4} {CS-X-X-X 4 3 7 2} {CW1-X-X-X 6 9 10 3} {CS-X-X-X 14 6 15 9} {CS-X-X-X 9 16 19 14} {CW1-X-X-X 26 10 14 16} {CW1-X-X-X 97 28 29 27} {CS-X-X-X 27 31 32 29} {CW1-X-X-X 33 28 31 30} {CS-X-X-X 34 29 30 31} {CW1-X-X-X 30 36 37 33} {CS-X-X-X 36 33 42 41} {CS-X-X-X 43 36 46 41} {CW1-X-X-X 43 37 48 41} {CW2-X-X-X 43 52 53 48} {CS-X-X-X 57 48 58 53} {CW2-X-X-X 56 0 57 52} {CS-X-X-X 57 53 56 59} {CW1-X-X-X 68 67 70 69} {CS-X-X-X 70 68 72 73} {CW1-X-X-X 71 69 72 74} {CS-X-X-X 72 71 75 70} {CW1-X-X-X 74 77 78 71} {CS-X-X-X 74 82 83 77} {CS-X-X-X 77 84 87 82} {CW1-X-X-X 84 82 89 78} {CW2-X-X-X 89 84 94 93} {CS-X-X-X 99 89 98 94} {CW2-X-X-X 97 93 98 27} {CS-X-X-X 98 94 97 100}}

As you can see, for some reason, the set impropers function reverses the order of some of the impropers. For instance, look at the first improper, which I've highlighted in bold.