VMD-L Mailing List

From: Evgeny Bulat (jack.bulat_at_gmail.com)
Date: Tue Jun 26 2012 - 16:42:37 CDT

Hey everyone! The issue of the autoionize plugin functioning incorrectly
for systems containing atom type names that are longer than 4 letters has
come up several times in VMD-L. In fact, it's recently come up for me as
well. I wasn't able to find code that would specifically take care of it,
and so decided to write some myself. In this message, I'm posting a Python
script that correctly adds up partial charges and provides them to the
autoionize plugin for further calculation. To make it all work, do the
following:

1) Go to the directory that contains the plugin. Wherever the VMD libraries
are stored in your system, continue on to
~/vmd/plugins/noarch/tcl/autoionize1.#. The # will be 2 or 3 depending on
the version you have (although I HIGHLY, HIGHLY suggest upgrading to vmd
1.9.1 if you haven't done so already, since there have been some crucial
updates made to the autoionize plugin, courtesy of Leonardo Trabuco). Once
you're there, paste a file called calc_charge.py that contains the
following:

///////
#calc_charge.py: Correctly calculates system charge based on charges listed
directly in the PSF
# file.
#Written for Python version 2; as far as Python 3 goes, I can't attest to
the effectiveness of this
#code in that environment.
#NOTE: Don't rely on PSF files that have been written off of the atomselect
function. In those,
#atom types with names longer than 4 letters have had their names clipped
off and, worse, their
#charges possibly ruined. Instead, use either X-Plor, converted CHARMM, or
psfgen-made files.
#
#Author: Jack Bulat, graduate student, UC Berkeley
#Date: June 26th, 2012

from sys import argv

script, psffile = argv

def start(psffile):

    f=open(psffile)
    lines = f.readlines()

    charge_count = float(0)
    start_read = False
    end_read = False

    # reads only the atom section of the PSF file, and skips everything else
    for line in lines:

        if "!NATOM" in line:
            start_read = True
            continue
        if "!NBOND" in line:
            end_read = True

        if (start_read) and (end_read):
            break

        if (start_read) and (not end_read):

            # reads partial charge and adds it on to cumulative charge
            words = line.strip().split(None, 8)
            if len(words) < 2:
                continue

            atom_charge = float(words[6])
            charge_count += atom_charge

    f.close()
    return charge_count

print start(psffile)
///////

2) Open up autoionize.tcl. Find all instances of the commands, "set
netCharge [eval "vecadd [$sel get charge]"]". There should be 2 total. In
the FIRST instance, comment off or delete the following block of code as I
have done already below:

  # Compute net charge of the system
# set sel [atomselect top all]
# set netCharge [eval "vecadd [$sel get charge]"]
# set roundNetCharge [expr round($netCharge)]
# $sel delete

Add in the following right below or in place of it:

#adjusted method for calculating charge correctly
set system $psffile
set netCharge [exec python calc_charge.py $system]
set roundNetCharge [expr round($netCharge)]

In the SECOND instance, comment off or delete the following block of code:

  # Re-compute net charge of the system
# set netCharge [eval "vecadd [$sel get charge]"]
# $sel delete

Add in the following right below or in place of it:

#adjusted method for calculating charge correctly
set system $prefix.psf
set netCharge [exec python calc_charge.py $system]

And...you're all set! Next time you run VMD, check to see that there are no
errors. This has worked on my machine, and I don't see any obvious reason
for why it can't work on others (but of course, I'm a little bit of a
novice in all this myself).

**Note: The readout for whether or not autoionize works correctly is the
net charge that's calculated by NAMD. That should be indicated in the
comments that NAMD spits out about the system. If you logged those
comments, grep "Info: TOTAL CHARGE =" and you should get it. As far as I
know, NAMD has NO issues with long atom type names, but VMD does. I could
be wrong but I don't think I am. If you know better, please let everyone
know!

I hope this helps!

Best,
Jack Bulat