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 gmx-users@gromacs.org 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 gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists