Skip to content

Commit

Permalink
Fixed blasted VDW bonding test; added TD_SCF stubs
Browse files Browse the repository at this point in the history
  • Loading branch information
@dylan-jayatilaka
dylan-jayatilaka committed Nov 24, 2024
1 parent 88c89e5 commit 0d738af
Show file tree
Hide file tree
Showing 9 changed files with 463 additions and 313 deletions.
4 changes: 3 additions & 1 deletion CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -251,8 +251,9 @@ set(FOO_SRC
${FOO_DIR}/crystal.foo
${FOO_DIR}/crystal23.foo
${FOO_DIR}/spherical_harmonic.foo
${FOO_DIR}/table_column.foo
${FOO_DIR}/t_tensor.foo
${FOO_DIR}/table_column.foo
${FOO_DIR}/td_scf.foo
${FOO_DIR}/atom_group.foo
${FOO_DIR}/vec{atom_group}.foo
${FOO_DIR}/vec{atom}.foo
Expand Down Expand Up @@ -408,6 +409,7 @@ set(FORTRAN_SRC
${CMAKE_CURRENT_BINARY_DIR}/time.F90
${CMAKE_CURRENT_BINARY_DIR}/types.F90
${CMAKE_CURRENT_BINARY_DIR}/t_tensor.F90
${CMAKE_CURRENT_BINARY_DIR}/td_scf.F90
${CMAKE_CURRENT_BINARY_DIR}/unit_cell.F90
${CMAKE_CURRENT_BINARY_DIR}/vec_atom_group.F90
${CMAKE_CURRENT_BINARY_DIR}/vec_atom.F90
Expand Down
236 changes: 137 additions & 99 deletions foofiles/atom.foo
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6040,7 +6040,7 @@ contains

! Bond distances and bonding

bond_distance_to(that) result (res) ::: PURE
bond_distance_to(that) result (res) ::: pure
! Return the distance between "self" and "that" other atom.
self :: IN
that :: ATOM, IN
Expand All @@ -6053,125 +6053,163 @@ contains

end

is_bonded_to(b,range_factor) result (res) ::: pure
! Return true if "self" is bonded to "b". If present, "range_factor" is used
! to determine the distance range in which the atoms are regarded as bonded.
! This uses the CCDC method, as documented on their web page.
is_nearby_to(that,dist) result (res) ::: pure
! Return TRUE if atom "self" and "that" atom are nearby
! i.e. within length "dist".
self :: IN
b :: ATOM, IN
range_factor :: REAL, IN, optional
res :: BIN
that :: ATOM, IN
dist :: REAL, IN
res :: BIN

tmp :: VEC{REAL}(3)
r2 :: REAL

tmp = .position - that.position
tmp = abs(tmp)

res = FALSE

if (tmp(1)>dist) return
if (tmp(2)>dist) return
if (tmp(3)>dist) return

r2 = dot_product(tmp,tmp)
res = r2 < dist*dist

end

t,bond,bond_min,bond_max,dx,dy,dz,r2 :: REAL
is_bonded_to(that,range_factor) result (res) ::: pure
! Return true if "self" is bonded to atom "that".
! "range_factor" is the distance range around the
! covalent radii sum where the atoms are still
! regarded as covalently bonded. This is the CCDC
! method. Radii from Brun et al (2011) Struct.
! Sci 67 p. 333-49.
self :: IN
that :: ATOM, IN
range_factor :: REAL, optional, IN
res :: BIN

dx,dy,dz,r2,t :: REAL
bond,bmin,bmax :: REAL
Za,Zb,as,bs :: INT

Za = .atomic_number
Zb = b.atomic_number
as = .site_disorder_group
bs = b.site_disorder_group
Za = self.atomic_number
Zb = that.atomic_number

if (Za<1) then; res = FALSE
else if (Zb<1) then; res = FALSE
else if (Za>covalent_radii_ccdc.dim) then; res = FALSE
else if (Zb>covalent_radii_ccdc.dim) then; res = FALSE
else if (as*bs>0 AND as/=bs) then; res = FALSE
else
res = FALSE

t = atom_bonded_range_factor
if (present(range_factor)) t = range_factor
if (Za<1) return
if (Zb<1) return

if (Za>covalent_radii_ccdc.dim) return
if (Zb>covalent_radii_ccdc.dim) return

as = self.site_disorder_group
bs = that.site_disorder_group
if (as*bs>0 AND as/=bs) return

bond = .covalent_radius_ccdc + b.covalent_radius_ccdc
bond_min = max(bond - t,ZERO)
bond_max = bond + t
t = atom_bonded_range_factor
if (present(range_factor)) t = range_factor

! For HH, nothing closer than 0.9A
if (self.atomic_number==1 AND b.atomic_number==1) bond_min = 0.7d0*BOHR_PER_ANGSTROM
bond = self.covalent_radius_ccdc &
+ that.covalent_radius_ccdc
bmin = max(bond - t,ZERO)
bmax = bond + t

dx = abs(.position(1) - b.position(1))
dy = abs(.position(2) - b.position(2))
dz = abs(.position(3) - b.position(3))
if (dx>bond_max) then; res = FALSE
else if (dy>bond_max) then; res = FALSE
else if (dz>bond_max) then; res = FALSE
else
r2 = dx*dx + dy*dy + dz*dz
res = (r2 < bond_max*bond_max) AND (r2 > bond_min*bond_min)
end
! For HH, nothing closer than 0.7A
if (Za==1 AND Zb==1) bmin = 0.7d0*BOHR_PER_ANGSTROM

end
dx = abs(.position(1) - that.position(1))
dy = abs(.position(2) - that.position(2))
dz = abs(.position(3) - that.position(3))

end
if (dx>bmax) return
if (dy>bmax) return
if (dz>bmax) return

r2 = dx*dx + dy*dy + dz*dz
res = r2<bmax*bmax AND r2>bmin*bmin

is_nearby_to(b,dist) result (res)
! Return TRUE if atom "self" and atom "b" are nearby, i.e. within length
! "dist".
b :: ATOM
dist :: REAL
res :: BIN
tmp :: VEC{REAL}(3)
r2 :: REAL
tmp = .position - b.position
tmp = abs(tmp)
if (tmp(1)>dist) then; res = FALSE
else if (tmp(2)>dist) then; res = FALSE
else if (tmp(3)>dist) then; res = FALSE
else
r2 = dot_product(tmp,tmp)
res = (r2 < dist*dist)
end
end

is_vdw_bonded_to(b,range_factor,vdw_range_pc) result (res) ::: pure
! Return true if "self" is vdw bonded to "b". Atoms which are
! covalently bonded are not vand-der-waals bonded. If present,
! "range_factor" is used to determine the distance range in which
! the atoms are regarded as covalent bonded. If present
! "vdw_range_pc" is the percentage increase in the vdw radius to be
! tolerated while still regarding the atoms vdw bonded (default
! 10%). This uses the CCDC method, as documented in
! Brun et al (2011) Struct. Sci 67 p. 333-49.
is_vdw_bonded_to(that,range_factor,vdw_range_pc) result (res) ::: pure
! Return true if "self" is VDW bonded to "that" atom.
! Atoms covalently bonded are not regarded VDW bonded.
! "range_factor" is the distance range around the
! covalent radii sum where the atoms are still
! regarded as covalently bonded. "vdw_range_pc" is the
! *percentage* increase in the vdw radius sum to be
! tolerated so that atoms are vdw bonded (default 10%).
self :: IN
b :: ATOM, IN
range_factor,vdw_range_pc :: REAL, IN, optional
that :: ATOM, IN
range_factor :: REAL, optional, IN
vdw_range_pc :: REAL, optional, IN
res :: BIN

t,bond,bond_max,dx,dy,dz,r2 :: REAL
Za,Zb,as,bs,ag,bg :: INT

Za = .atomic_number
Zb = b.atomic_number
as = .site_disorder_group
bs = b.site_disorder_group
ag = .group
bg = b.group

if (Za<1) then; res = FALSE
else if (Zb<1) then; res = FALSE
else if (Za>covalent_radii_ccdc.dim) then; res = FALSE
else if (Zb>covalent_radii_ccdc.dim) then; res = FALSE
else if (as*bs>0 AND as==bs) then; res = FALSE
else if (ag*bg>0 AND ag==bg) then; res = FALSE
else if (.is_bonded_to(b,range_factor)) then; res = FALSE
else
dx,dy,dz,r2, t :: REAL
bond,bmin,bmax :: REAL
Za,Zb, as,bs :: INT
! ag,bg :: INT
cov :: BIN

t = ONE+0.01*atom_vdw_bonded_range_pc
if (present(vdw_range_pc)) t = ONE + 0.01*vdw_range_pc
Za = self.atomic_number
Zb = that.atomic_number

bond = .vdw_radius_ccdc + b.vdw_radius_ccdc
bond_max = bond*t
res = FALSE

dx = abs(.position(1) - b.position(1))
dy = abs(.position(2) - b.position(2))
dz = abs(.position(3) - b.position(3))
if (dx>bond_max) then; res = FALSE
else if (dy>bond_max) then; res = FALSE
else if (dz>bond_max) then; res = FALSE
else
r2 = dx*dx + dy*dy + dz*dz
res = (r2 < bond_max*bond_max)
end
if (Za<1) return
if (Zb<1) return

end
if (Za>vdw_radii_ccdc.dim) return
if (Zb>vdw_radii_ccdc.dim) return

as = self.site_disorder_group
bs = that.site_disorder_group
if (as*bs>0 AND as==bs) return

! This does not work for 1 group!
! ag = self.group
! bg = that.group
! if (ag*bg>0 AND ag==bg) return

! Is it covalently bonded?
! WARNING: must be the same as .is_bonded_to
t = atom_bonded_range_factor
if (present(range_factor)) t = range_factor

bond = self.covalent_radius_ccdc &
+ that.covalent_radius_ccdc
bmin = max(bond - t,ZERO)
bmax = bond + t

! For HH, nothing closer than 0.7A
if (Za==1 AND Zb==1) bmin = 0.7d0*BOHR_PER_ANGSTROM

dx = abs(.position(1) - that.position(1))
dy = abs(.position(2) - that.position(2))
dz = abs(.position(3) - that.position(3))

! Can't be vdw if covalent
r2 = dx*dx + dy*dy + dz*dz
cov = r2<bmax*bmax AND r2>bmin*bmin
if (cov) return

! Still here? Now check if VDW bonded
t = ONE + TOL(2)*atom_vdw_bonded_range_pc
if (present(vdw_range_pc)) t = ONE + TOL(2)*vdw_range_pc

bond = self.vdw_radius_ccdc &
+ that.vdw_radius_ccdc
bmax = bond*t

if (dx>bmax) return
if (dy>bmax) return
if (dz>bmax) return

! Finally ...
r2 = dx*dx + dy*dy + dz*dz
res = r2 < bmax*bmax

end

Expand Down
42 changes: 21 additions & 21 deletions foofiles/crystal.foo
View file
Original file line number Diff line number Diff line change
Expand Up @@ -743,27 +743,27 @@ contains

end

read_asymmetric_unit_geometry ::: private
! Read in the asymmetric unit geometry in crystal coordinates.
! The coordinates are read in as a single vector ordered as x,y,z
! incrementing fastest, for the first to the last atom.
self :: INOUT
geometry :: VEC{REAL}@

! WARN_IF(.asymmetric_unit_geometry.allocated,"asymmetric_unit_geometry exists!")

stdin.read_all(geometry)
ENSURE(modulo(geometry.dim,3)==0,"number of coordinates must be divisible by 3")

.n_asymmetric_unit_atoms = geometry.dim/3
.asymmetric_unit_geometry.destroy
.asymmetric_unit_geometry.create(3,.n_asymmetric_unit_atoms)
.asymmetric_unit_geometry = reshape(geometry,[3,.n_asymmetric_unit_atoms])
geometry.destroy

.asymmetric_unit_source = "manual-input"

end
! read_asymmetric_unit_geometry ::: private
! ! Read in the asymmetric unit geometry in crystal coordinates.
! ! The coordinates are read in as a single vector ordered as x,y,z
! ! incrementing fastest, for the first to the last atom.
! self :: INOUT
! geometry :: VEC{REAL}@
!
! ! WARN_IF(.asymmetric_unit_geometry.allocated,"asymmetric_unit_geometry exists!")
!
! stdin.read_all(geometry)
! ENSURE(modulo(geometry.dim,3)==0,"number of coordinates must be divisible by 3")
!
! .n_asymmetric_unit_atoms = geometry.dim/3
! .asymmetric_unit_geometry.destroy
! .asymmetric_unit_geometry.create(3,.n_asymmetric_unit_atoms)
! .asymmetric_unit_geometry = reshape(geometry,[3,.n_asymmetric_unit_atoms])
! geometry.destroy
!
! .asymmetric_unit_source = "manual-input"
!
! end

! Cluster keywords

Expand Down
1 change: 1 addition & 0 deletions foofiles/molecule.main.foo
View file
Original file line number Diff line number Diff line change
Expand Up @@ -545,6 +545,7 @@ contains
! case ("spherically_averaged_sf "); .MISC:spherically_averaged_sf
! case ("test_dump_file "); MOLECULE.MAIN:test_dump_file
! case ("test_plot_info "); .:test_plot_info
case ("test_vdw_bonding "); .MISC:test_vdw_bonding

case default ; UNKNOWN(keyword)

Expand Down
Loading

0 comments on commit 0d738af

Please sign in to comment.