[gmx-users] some code for writing topologies
laura.leay at postgrad.manchester.ac.uk
Mon Oct 7 10:55:07 CEST 2013
Below is some code I wrote a few years ago during my PhD. I started with a list of bonds in a file called 'bonds' and t he used this list to automatically create the lists of angled and proper dihedrals. I went on to edit this code to create a topology for entire polymers based on a single monomer which you can figure out for yourself if need be.
The list of bonds was created by eye and the structure is important to ensure that the angles and dihedrals are found correctly. The list was structured so that I started with the first atom in the gro file, looked at all atoms it was bonded to and then moved on the atom number 2. I did not count bonds where the atom bonded to this 'atom of interest' had a lower number than this 'atom of interest'.
e.g. if atom 1 is bonded to atoms 2, 3, and 4 the first 3 entries in the bonds list read:
If atom 2 is bonded to atoms 1, 5 and 6 then there are only two entries for this atom:
This avoids having the same bond counted twice and it is essential that double counting of bonds is avoided. Note that if the molecule already has a sensible structure (i.e. bonds and angles are not too small or too large) then you could easily write some code to find the list of bonds for you, just specify that a bond occurs anywhere where two atoms are below a certain distance apart.
I hope this helps some of you with your topologies. I realise that there are some tools already available for this but none of them worked for the polymers I was modelling. I imagine I am not the only one to have had this problem.
Here is the code:------------------------------>
integer,parameter:: totalbonds = 603
integer,parameter:: atoms =514 !number of atoms in molecule
integer,parameter:: maxbond=4 !max no of bonds per atom
integer:: i,j,k,n !do loops
integer:: iatom(totalbonds), jatom(totalbonds) !bond list
integer:: ia(atoms,maxbond+1), ja, ka(maxbond+1) !angles
integer:: id, jd, kd, ld
integer:: count=0 !counts the number of bonds for an atom
integer:: atombonds(atoms) !list of the number of bonds for each atom
read(10,*) iatom(i), jatom(i)
! this takes the list of bonds which has been read in and finds two bonds that
! form and angle ijk. Essentially it looks through each atom in turn, finds all the
! atoms that are bonded to it and stores then in the array ia(:,:).
! This array or list is then used to construct the list of angles ijk.
write(20,*) '[ angles ]'
do ja=1,atoms !atom at centre of angle
count=0 !counter for number of bonds, this is reset for each atom
do i=1,totalbonds !loop to find all atoms in the bond list bonded to atom ja
if(iatom(i)==ja)then !look through first bond list
print*, count, ja, jatom(i) !prints out to screen so the user can check it
if(jatom(i)==ja)then !look through second bond list
print*,count, ja, iatom(i) !prints out to screen so the user can check it
end do !i loop , bonded atoms
atombonds(ja)=count !add the count of bonds to the bonds list
print*, ja, count, ' ::ia:', ia(ja,1:count) !prints out as a check if the program should stall
if(count==2) then !write the angles to the topology file, if there are only two atoms bonded to atom j then there is only one angle
write(20,*)ia(ja,1), ja, ia(ja,2)
print*,' ',ia(ja,1), ja, ia(ja,2)
else if(count.ge.3) then !if there are 3 or more atoms bonded to atom j then there are multiple angles
end do !ja loop
! this uses the list of angles ia(:,:) found above to construct the list of dihedral angles.
! It looks through the original bond list that was read in and finds how many angles (or bonds)
! are created about each atom at the end of each bond. This information is then used to construct
! a list of dihedral angles ijkl where j and k are form the original bond from the bonds list.
write(20,*) ' '
write(20,*) '[ dihedrals ]'
!ia(atoms,maxcount) records each atom that is bonded to a given atom
!atombonds(atoms) records the number of bonds for a given atom
do n=1,totalbonds !look through the bonds list
jd=iatom(n) !set jd equal to the nth entry in the bonds list
kd=jatom(n) !set kd equal to the nth entry in the bonds list
if(atombonds(jd)==1)cycle !check to see if this bond is the middle of a dihedral
if(atombonds(kd)==1)cycle !i.e. if the number of bonds on the atom is equal to 1 it must be at the end of a dihedral and so can be ignored
do i=1,atombonds(jd) !loop through the number of bonds for atom 'jd'
id=ia(jd,i) !set the atom 'id' equal to atom 'ia' which is atom i in the angle ijk from above
if(id==kd)cycle !we're not interested in 'kd' as it will be taken care of on the next line
do j=1,atombonds(kd) !loop through the number of bonds for atom 'kd'
ld=ia(kd,j) !set atom 'ld' equal to atom 'ia' which is atom i in the angle ijk from above
if(ld==jd)cycle !we've already accounted for atom 'jd' above so we can ignore it now
write(20,*)id,jd,kd,ld !writes the atom numbers for the dihedral ijkl to the topology file
end do !j loop
end do !i loop
end do !n loop, bonds
end program topol
More information about the gromacs.org_gmx-users