********>Bugfix 19: Author: Ross Walker Date: May 8, 2009 Programs: sander Description: The routine for identifying elements in QMMM runs can sometimes misidentify elements. This is most likely when dealing with transition metals such as Cadmium or Selenium. This patch modifies the routine to use a combination of the first letter of the atom name combined with the atomic mass. Use this patch to update amber10/src/sander/qmmm_module.f ------------------------------------------------------------------------------ --- qmmm_module.f 2009-05-07 23:29:58.000000000 -0700 +++ qmmm_module.f 2009-05-07 23:40:13.000000000 -0700 @@ -1638,7 +1638,7 @@ REQUIRE(ier == 0) do i=1,natom - call get_atomic_number(ih(m04+i-1),qmmm_struct%all_atom_numbers(i)) + call get_atomic_number(ih(m04+i-1),x(lmass+i-1),qmmm_struct%all_atom_numbers(i)) enddo endif #endif @@ -1647,10 +1647,10 @@ qmmm_nml%iqmatoms(i) = iqmatoms(i) !Get the atomic numbers (used to be done in rdparm2...) #ifdef SQM - call get_atomic_number( atnam(i), qmmm_struct%iqm_atomic_numbers(i) ) + call get_atomic_number( atnam(i), x(lmass+i-1), qmmm_struct%iqm_atomic_numbers(i) ) #else j = iqmatoms(i) - call get_atomic_number( ih(m04+j-1),qmmm_struct%iqm_atomic_numbers(i) ) + call get_atomic_number( ih(m04+j-1), x(lmass+j-1), qmmm_struct%iqm_atomic_numbers(i) ) #endif end do @@ -1838,98 +1838,288 @@ !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+Identify qm atom elements. - subroutine get_atomic_number(name,atomic_number) + subroutine get_atomic_number(atom_name,atom_mass,atomic_number) ! ! ! This subroutine assigns atomic number based upon the first -! letter of the atom symbol read in from the topology file. Because -! this symbol is set in link, this procedure is ripe for misreading -! errors, but is the best way to go about it right now. +! letter of the atom symbol read in from the topology file and the mass. +! The assumption here is that an atom name matches the first letter of the +! element. If it doesn't then this routine will need to be modified. -!NOTE: THIS ROUTINE NEEDS COMPLETELY RE-WRITING (ROSS WALKER)... -! -! name ==> character containing the atomic name +! atom_name ==> character containing the atomic name ! iqm_atomic_numbers ==> integer array of atomic numbers assigned to atoms ! implicit none !passed in integer :: atomic_number - character(len=4) :: name + _REAL_ :: atom_mass + character(len=4) :: atom_name !local integer :: i, j - character(len=2) :: elemnt(107) - data (elemnt(i),i=1,107)/'H','XX', & - 'Li','XX','B','C','N','O','F','XX', & - 'XX','Mg','Al','XX','P','S','XX','XX', & - 'XX','XX','XX','XX','XX','XX','XX','XX','XX','XX','XX', & - 'Z','Ga','XX','XX','XX','XX','XX', & - 'XX','XX','XX','XX','XX','XX','XX','XX','XX','XX','XX', & - 'XX','XX','XX','XX','XX','I','XX', & - 'XX','XX','XX','XX','XX','XX','XX','XX','XX','XX','XX','XX', & - 'XX','XX','XX','XX','XX','XX','XX','W','XX','XX','XX','XX', & - 'XX','XX','XX','XX','XX','XX','XX','XX', & - 'XX','XX','XX','XX','XX','U','XX','XX','XX','XX','XX','XX','XX',& - 'XX','XX','XX','++','+','--','-','TV'/ -! -! compare values in name to those in elemnt and assign atomic -! numbers accordingly -! - do j=1,107 - if(name(1:1) .eq. elemnt(j)(1:1)) then - if(j.eq.6) then - if(name(2:2).eq.'L' .OR. name(2:2).eq.'l') then - atomic_number = 17 - go to 25 - else - atomic_number = 6 - go to 25 - end if - elseif(j.eq.5)then - if(name(2:2).eq.'R' .OR. name(2:2).eq.'r')then - atomic_number = 35 - go to 25 - else if(name(2:2).eq.'E' .OR. name(2:2).eq.'e')then - atomic_number = 4 - go to 25 - else - atomic_number = 5 - go to 25 - endif - elseif(j.eq.16) then - if (name(2:2) .eq. 'I' .OR. name(2:2) .eq. 'i') then - atomic_number = 14 !Silicon - go to 25 - else - atomic_number = 16 !Sulphur - go to 25 - end if - elseif (j.eq.31) then - if (name(2:2) .eq. 'E' .OR. name(2:2) .eq. 'e') then - atomic_number = 32 !Germanium - go to 25 - else - atomic_number = 31 !Gallium - go to 25 - end if - elseif (j.eq.13) then - if (name(2:2) .eq. 'S' .OR. name(2:2) .eq. 's') then - atomic_number = 33 !As - go to 25 - else - atomic_number = 13 - go to 25 - end if - else - atomic_number = j - go to 25 - end if - end if - end do - write(6,*) 'Unable to correctly identify element ', name - call mexit(6,1) - 25 continue + +!Nobel Gasses are not supported. +!Lanthanides are not supported. +!Actinides are not supported. + + if(atom_name(1:1) .eq. 'a' .or. atom_name(1:1) .eq. 'A') then + if(atom_mass > 24.0d0 .and. atom_mass <= 28.0d0) then + atomic_number = 13 !Aluminium + elseif(atom_mass > 73.0d0 .and. atom_mass <= 77.0d0) then + atomic_number = 33 !Arsenic + elseif(atom_mass > 106.0d0 .and. atom_mass <= 109.0d0) then + atomic_number = 47 !Silver + elseif(atom_mass > 195.0d0 .and. atom_mass <= 199.0d0) then + atomic_number = 79 !Gold + elseif(atom_mass > 208.0d0 .and. atom_mass <= 212.0d0) then + atomic_number = 85 !Astatine + else + write(6,*) 'Unable to correctly identify element ', atom_name + call mexit(6,1) + endif + elseif(atom_name(1:1) .eq. 'b' .or. atom_name(1:1) .eq. 'B') then + if(atom_mass > 8.0d0 .and. atom_mass <= 10.0d0) then + atomic_number = 4 !Beryllium + elseif(atom_mass > 10.0d0 .and. atom_mass <= 12.0d0) then + atomic_number = 5 !Boron + elseif(atom_mass > 77.0d0 .and. atom_mass <= 81.0d0) then + atomic_number = 35 !Bromine + elseif(atom_mass > 135.0d0 .and. atom_mass <= 139.0d0) then + atomic_number = 56 !Barium + elseif(atom_mass > 207.0d0 .and. atom_mass <= 211.0d0) then + atomic_number = 83 !Bismuth + else + write(6,*) 'Unable to correctly identify element ', atom_name + call mexit(6,1) + endif + elseif(atom_name(1:1) .eq. 'c' .or. atom_name(1:1) .eq. 'C') then + if(atom_mass > 10.0d0 .and. atom_mass <= 14.0d0) then + atomic_number = 6 !Carbon + elseif(atom_mass > 33.0d0 .and. atom_mass <= 37.0d0) then + atomic_number = 17 !Chlorine + elseif(atom_mass > 38.0d0 .and. atom_mass <= 42.0d0) then + atomic_number = 20 !Calcium + elseif(atom_mass > 50.0d0 .and. atom_mass <= 54.0d0) then + atomic_number = 24 !Chromium + elseif(atom_mass > 57.0d0 .and. atom_mass <= 61.0d0) then + atomic_number = 27 !Cobalt + elseif(atom_mass > 61.0d0 .and. atom_mass <= 65.0d0) then + atomic_number = 29 !Copper + elseif(atom_mass > 110.0d0 .and. atom_mass <= 114.0d0) then + atomic_number = 48 !Cadmium + elseif(atom_mass > 131.0d0 .and. atom_mass <= 135.0d0) then + atomic_number = 55 !Cesium + else + write(6,*) 'Unable to correctly identify element ', atom_name + call mexit(6,1) + endif + elseif(atom_name(1:1) .eq. 'd' .or. atom_name(1:1) .eq. 'D') then + write(6,*) 'Unable to correctly identify element ', atom_name + call mexit(6,1) + elseif(atom_name(1:1) .eq. 'e' .or. atom_name(1:1) .eq. 'E') then + write(6,*) 'Unable to correctly identify element ', atom_name + call mexit(6,1) + elseif(atom_name(1:1) .eq. 'f' .or. atom_name(1:1) .eq. 'F') then + if(atom_mass > 17.0d0 .and. atom_mass <= 21.0d0) then + atomic_number = 9 !Fluorine + elseif(atom_mass > 54.0d0 .and. atom_mass <= 58.0d0) then + atomic_number = 26 !Iron + elseif(atom_mass > 218.0d0 .and. atom_mass <= 228.0d0) then + atomic_number = 87 !Francium + else + write(6,*) 'Unable to correctly identify element ', atom_name + call mexit(6,1) + endif + elseif(atom_name(1:1) .eq. 'g' .or. atom_name(1:1) .eq. 'G') then + if(atom_mass > 67.0d0 .and. atom_mass <= 71.0d0) then + atomic_number = 31 !Gallium + elseif(atom_mass > 71.0d0 .and. atom_mass <= 75.0d0) then + atomic_number = 32 !Germanium + else + write(6,*) 'Unable to correctly identify element ', atom_name + call mexit(6,1) + endif + elseif(atom_name(1:1) .eq. 'h' .or. atom_name(1:1) .eq. 'H') then + if(atom_mass > 0.0d0 .and. atom_mass <= 2.0d0) then + atomic_number = 1 !Hydrogen + elseif(atom_mass > 3.0d0 .and. atom_mass <= 5.0d0) then + atomic_number = 2 !Helium + elseif(atom_mass > 176.0d0 .and. atom_mass <= 180.0d0) then + atomic_number = 72 !Hafnium + elseif(atom_mass > 198.0d0 .and. atom_mass <= 202.0d0) then + atomic_number = 80 !Mercury + else + write(6,*) 'Unable to correctly identify element ', atom_name + call mexit(6,1) + endif + elseif(atom_name(1:1) .eq. 'i' .or. atom_name(1:1) .eq. 'I') then + if(atom_mass > 112.0d0 .and. atom_mass <= 116.0d0) then + atomic_number = 49 !Indium + elseif(atom_mass > 125.0d0 .and. atom_mass <= 129.0d0) then + atomic_number = 53 !Iodine + elseif(atom_mass > 190.0d0 .and. atom_mass <= 194.0d0) then + atomic_number = 77 !Iridium + else + write(6,*) 'Unable to correctly identify element ', atom_name + call mexit(6,1) + endif + elseif(atom_name(1:1) .eq. 'j' .or. atom_name(1:1) .eq. 'J') then + write(6,*) 'Unable to correctly identify element ', atom_name + call mexit(6,1) + elseif(atom_name(1:1) .eq. 'k' .or. atom_name(1:1) .eq. 'K') then + if(atom_mass > 37.0d0 .and. atom_mass <= 41.0d0) then + atomic_number = 19 !Potassium + else + write(6,*) 'Unable to correctly identify element ', atom_name + call mexit(6,1) + endif + elseif(atom_name(1:1) .eq. 'l' .or. atom_name(1:1) .eq. 'L') then + if(atom_mass > 6.0d0 .and. atom_mass <= 8.0d0) then + atomic_number = 3 !Lithium + else + write(6,*) 'Unable to correctly identify element ', atom_name + call mexit(6,1) + endif + elseif(atom_name(1:1) .eq. 'm' .or. atom_name(1:1) .eq. 'M') then + if(atom_mass > 22.0d0 .and. atom_mass <= 26.0d0) then + atomic_number = 12 !Magnesium + elseif(atom_mass > 53.0d0 .and. atom_mass <= 57.0d0) then + atomic_number = 25 !Manganese + elseif(atom_mass > 94.0d0 .and. atom_mass <= 98.0d0) then + atomic_number = 42 !Molybdenem + else + write(6,*) 'Unable to correctly identify element ', atom_name + call mexit(6,1) + endif + elseif(atom_name(1:1) .eq. 'n' .or. atom_name(1:1) .eq. 'N') then + if(atom_mass > 13.0d0 .and. atom_mass <= 15.0d0) then + atomic_number = 7 !Nitrogen + elseif(atom_mass > 21.0d0 .and. atom_mass <= 25.0d0) then + atomic_number = 11 !Sodium + elseif(atom_mass > 57.0d0 .and. atom_mass <= 61.0d0) then + atomic_number = 28 !Nickel + elseif(atom_mass > 95.0d0 .and. atom_mass <= 99.0d0) then + atomic_number = 41 !Niobium + else + write(6,*) 'Unable to correctly identify element ', atom_name + call mexit(6,1) + endif + elseif(atom_name(1:1) .eq. 'o' .or. atom_name(1:1) .eq. 'O') then + if(atom_mass > 14.0d0 .and. atom_mass <= 18.0d0) then + atomic_number = 8 !Oxygen + elseif(atom_mass > 188.0d0 .and. atom_mass <= 192.0d0) then + atomic_number = 76 !Osmium + else + write(6,*) 'Unable to correctly identify element ', atom_name + call mexit(6,1) + endif + elseif(atom_name(1:1) .eq. 'p' .or. atom_name(1:1) .eq. 'P') then + if(atom_mass > 29.0d0 .and. atom_mass <= 33.0d0) then + atomic_number = 15 !Phosphorus + elseif(atom_mass > 104.0d0 .and. atom_mass <= 108.0d0) then + atomic_number = 46 !Palladium + elseif(atom_mass > 193.0d0 .and. atom_mass <= 197.0d0) then + atomic_number = 78 !Platinum + elseif(atom_mass > 205.0d0 .and. atom_mass <= 208.0d0) then + atomic_number = 82 !Lead + elseif(atom_mass > 208.0d0 .and. atom_mass <= 212.0d0) then + atomic_number = 84 !Polonium + else + write(6,*) 'Unable to correctly identify element ', atom_name + call mexit(6,1) + endif + elseif(atom_name(1:1) .eq. 'q' .or. atom_name(1:1) .eq. 'Q') then + write(6,*) 'Unable to correctly identify element ', atom_name + call mexit(6,1) + elseif(atom_name(1:1) .eq. 'r' .or. atom_name(1:1) .eq. 'R') then + if(atom_mass > 84.0d0 .and. atom_mass <= 88.0d0) then + atomic_number = 37 !Rubidium + elseif(atom_mass > 99.0d0 .and. atom_mass <= 102.0d0) then + atomic_number = 44 !Ruthenium + elseif(atom_mass > 102.0d0 .and. atom_mass <= 105.0d0) then + atomic_number = 45 !Rhodium + elseif(atom_mass > 184.0d0 .and. atom_mass <= 188.0d0) then + atomic_number = 75 !Rhenium + elseif(atom_mass > 220.0d0 .and. atom_mass <= 230.0d0) then + atomic_number = 88 !Radium + else + write(6,*) 'Unable to correctly identify element ', atom_name + call mexit(6,1) + end if + elseif(atom_name(1:1) .eq. 's' .or. atom_name(1:1) .eq. 'S') then + if(atom_mass > 26.0d0 .and. atom_mass <= 30.0d0) then + atomic_number = 14 !Silicon + elseif(atom_mass > 30.0d0 .and. atom_mass <= 34.0d0) then + atomic_number = 16 !Sulphur + elseif(atom_mass > 43.0d0 .and. atom_mass <= 47.0d0) then + atomic_number = 21 !Scandium + elseif(atom_mass > 77.0d0 .and. atom_mass <= 81.0d0) then + atomic_number = 34 !Selenium + elseif(atom_mass > 86.0d0 .and. atom_mass <= 89.0d0) then + atomic_number = 38 !Strontium + elseif(atom_mass > 116.0d0 .and. atom_mass <= 120.0d0) then + atomic_number = 50 !Tin + elseif(atom_mass > 120.0d0 .and. atom_mass <= 124.0d0) then + atomic_number = 51 !Antimony + else + write(6,*) 'Unable to correctly identify element ', atom_name + call mexit(6,1) + end if + elseif(atom_name(1:1) .eq. 't' .or. atom_name(1:1) .eq. 'T') then + if(atom_mass > 46.0d0 .and. atom_mass <= 50.0d0) then + atomic_number = 22 !Titanium + elseif(atom_mass > 96.0d0 .and. atom_mass <= 100.0d0) then + atomic_number = 43 !Technetium + elseif(atom_mass > 125.0d0 .and. atom_mass <= 130.0d0) then + atomic_number = 52 !Tellurium + elseif(atom_mass > 179.0d0 .and. atom_mass <= 183.0d0) then + atomic_number = 73 !Tantalum + elseif(atom_mass > 201.0d0 .and. atom_mass <= 206.0d0) then + atomic_number = 81 !Thallium + else + write(6,*) 'Unable to correctly identify element ', atom_name + call mexit(6,1) + endif + elseif(atom_name(1:1) .eq. 'u' .or. atom_name(1:1) .eq. 'U') then + write(6,*) 'Unable to correctly identify element ', atom_name + call mexit(6,1) + elseif(atom_name(1:1) .eq. 'v' .or. atom_name(1:1) .eq. 'V') then + if(atom_mass > 49.0d0 .and. atom_mass <= 53.0d0) then + atomic_number = 23 !Vanadium + else + write(6,*) 'Unable to correctly identify element ', atom_name + call mexit(6,1) + endif + elseif(atom_name(1:1) .eq. 'w' .or. atom_name(1:1) .eq. 'W') then + if(atom_mass > 179.0d0 .and. atom_mass <= 183.0d0) then + atomic_number = 74 !Tungsten + else + write(6,*) 'Unable to correctly identify element ', atom_name + call mexit(6,1) + endif + elseif(atom_name(1:1) .eq. 'x' .or. atom_name(1:1) .eq. 'X') then + write(6,*) 'Unable to correctly identify element ', atom_name + call mexit(6,1) + elseif(atom_name(1:1) .eq. 'y' .or. atom_name(1:1) .eq. 'Y') then + if(atom_mass > 86.0d0 .and. atom_mass <= 90.0d0) then + atomic_number = 39 !Yttrium + else + write(6,*) 'Unable to correctly identify element ', atom_name + call mexit(6,1) + end if + elseif(atom_name(1:1) .eq. 'z' .or. atom_name(1:1) .eq. 'Z') then + if(atom_mass > 89.0d0 .and. atom_mass <= 93.0d0) then + atomic_number = 40 !Zirconium + else + write(6,*) 'Unable to correctly identify element ', atom_name + call mexit(6,1) + end if + else + write(6,*) 'Unable to correctly identify element ', atom_name + call mexit(6,1) + endif + return end subroutine get_atomic_number ------------------------------------------------------------------------------ Temporary Workarounds: None