[gmx-users] some code for writing topologies

[Laura Leay]

All,

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:

1         2

1         3
1      4
If atom 2 is bonded to atoms 1, 5 and 6 then there are only two entries for this atom:

2         5

2         6
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:------------------------------>



program topol

implicit none

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

open(10,file='bonds')

do i=1,totalbonds
  read(10,*) iatom(i), jatom(i)
  !print*,iatom(i), jatom(i)
end do

ia(:,:)=0
open(20,file='topology')


!----------------------------------------------------angles
!  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,*)' '
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
  ! ka(:)=0
   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
         count=count+1
         print*, count, ja, jatom(i)    !prints out to screen so the user can check it
         ia(ja,count)=jatom(i)
       end if
       if(jatom(i)==ja)then     !look through second bond list
         count=count+1
         print*,count, ja, iatom(i)             !prints out to screen so the user can check it
         ia(ja,count)=iatom(i)
       end if
   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
     do i=1,count
       do k=1,count
         if(k.le.i)cycle
         write(20,*)ia(ja,i),ja,ia(ja,k)
         print*,'  ',ia(ja,i),ja,ia(ja,k)
       end do
     end do
   end if
print*, '-----------------------'
end do !ja loop


!------------------------------------------------------------dihedrals
! 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
     ! print*,id,jd,kd,ld
    end do      !j loop
  end do !i loop
  !pause
end do !n loop, bonds

end program topol