From: Axel Kohlmeyer (
Date: Fri May 15 2009 - 10:28:59 CDT

On Thu, 2009-05-14 at 21:01 -0400, Maxim Paliy. wrote:
> Dear VMD experts
> I want to add a few bonds to my system (join two nucleic backbones)
> and then write the corresponding psf via e.g.
> set all [atomselect top all]
> $all writepsf my.psf
> The change in bonding is reflected this way, but the angles and
> dihedrals stay intact..
> How would I force the angles and dihedrals to update as well?

as john said, there is no automatic way (yet), but it should be
possible to script it. the new topotools package in the latest
beta versions should be an option to do this in a straightforward
way. in fact, it was in part written with applications like this,
i.e. make some small modification to existing topologies in mind.

you can type: package require topotools ; topo help
to get an overview of the available commands. off the
top of my head, you could do something along the following.

you need to know the index numbers of the atoms that form
a new bond <b1> <b2>:

package require topotools 1.0
set bnd1 <b1> ; # this is the first atom of the bond
set bnd2 <b2> ; # this is the second atoms of the bond

set bond [atomselect top "index $bnd1 $bnd2"]

# get list of atoms that each bonding partner is bound to
lassign [$sel getbonds] blist1 blist2

# construct angles
foreach b $blist1 {
  topo addangle $b $bnd1 $bnd2
foreach b $blist2 {
  topo addangle $b $bnd2 $bnd1

# construct dihedrals
foreach b1 $blist1 {
    foreach b2 $blist2 {
       topo adddihedral $b1 $bnd1 $bnd2 $b2

animate write psf newstruct.psf

this assumes that there is no overlap between
the lists of atoms that the bonding atoms are
bound to, which should be true in your case.
but depending on how you hook up stuff, it may
lead to rings that require impropers instead
of dihedrals, and that you have to decide from
visual inspection.

please give it a try and let me know, if this
works out, since then i would try to include a
version of this scriptlet into as an add-on
to topotools.


> For some reason, the simplest way
> readpsf 1.psf
> coordpdb 1.pdb
> readpsf 2.psf
> coordpdb 2.pdb
> writepsf my.psf
> writepdb my.pdb
> does not work for me (probably because the ends of the chains I want to
> join are not well aligned).
> maybe there are some other ways to accomplish this task without writing
> all psf file manually?
> Thanks a lot in advance for the suggestions,
> Maxim

Axel Kohlmeyer
   Center for Molecular Modeling   --   University of Pennsylvania
Department of Chemistry, 231 S.34th Street, Philadelphia, PA 19104-6323
tel: 1-215-898-1582,  fax: 1-215-573-6233,  office-tel: 1-215-898-5425
If you make something idiot-proof, the universe creates a better idiot.