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
--
gmx-users mailing list [email protected]
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the
www interface or send it to [email protected].
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
