Skip to content
This repository was archived by the owner on Dec 8, 2024. It is now read-only.

Commit d570637

Browse files
larsbratholmandersx
authored andcommitted
option to calculate only subset of atoms in local representations (#26)
* option to calculate only subset of atoms in local representations * added reference calculation of representations * finalized reference testing * fixed error in test
1 parent 370ea56 commit d570637

13 files changed

+953
-391
lines changed

qml/compound.py

+12-5
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
# MIT License
22
#
3-
# Copyright (c) 2016 Anders Steen Christensen, Felix Faber
3+
# Copyright (c) 2016-2017 Anders Steen Christensen, Felix Faber, Lars Andersen Bratholm
44
#
55
# Permission is hereby granted, free of charge, to any person obtaining a copy
66
# of this software and associated documentation files (the "Software"), to deal
@@ -70,7 +70,7 @@ def __init__(self, xyz = None):
7070
if xyz is not None:
7171
self.read_xyz(xyz)
7272

73-
def generate_coulomb_matrix(self, size = 23, sorting = "row-norm"):
73+
def generate_coulomb_matrix(self, size = 23, sorting = "row-norm", indices = None):
7474
""" Creates a Coulomb Matrix representation of a molecule.
7575
A matrix :math:`M` is constructed with elements
7676
@@ -132,7 +132,8 @@ def generate_eigenvalue_coulomb_matrix(self, size = 23):
132132
self.nuclear_charges, self.coordinates, size = size)
133133

134134
def generate_atomic_coulomb_matrix(self, size = 23, sorting = "row-norm",
135-
central_cutoff = 1e6, central_decay = -1, interaction_cutoff = 1e6, interaction_decay = -1):
135+
central_cutoff = 1e6, central_decay = -1, interaction_cutoff = 1e6, interaction_decay = -1,
136+
indices = None):
136137
""" Creates a Coulomb Matrix representation of the local environment of a central atom.
137138
For each central atom :math:`k`, a matrix :math:`M` is constructed with elements
138139
@@ -178,6 +179,11 @@ def generate_atomic_coulomb_matrix(self, size = 23, sorting = "row-norm",
178179
179180
The upper triangular of M, including the diagonal, is concatenated to a 1D
180181
vector representation.
182+
183+
The representation can be calculated for a subset by either specifying
184+
``indices = [0,1,...]``, where :math:`[0,1,...]` are the requested atom indices,
185+
or by specifying ``indices = 'C'`` to only calculate central carbon atoms.
186+
181187
The representation is calculated using an OpenMP parallel Fortran routine.
182188
183189
:param size: The size of the largest molecule supported by the representation
@@ -194,6 +200,8 @@ def generate_atomic_coulomb_matrix(self, size = 23, sorting = "row-norm",
194200
:type interaction_cutoff: float
195201
:param interaction_decay: The distance over which the the coulomb interaction decays from full to none
196202
:type interaction_decay: float
203+
:param indices: Subset indices or atomtype
204+
:type indices: Nonetype/array/string
197205
198206
199207
:return: nD representation - shape (:math:`N_{atoms}`, size(size+1)/2)
@@ -202,7 +210,7 @@ def generate_atomic_coulomb_matrix(self, size = 23, sorting = "row-norm",
202210

203211

204212
self.representation = generate_atomic_coulomb_matrix(
205-
self.nuclear_charges, self.coordinates, size = size,
213+
self.nuclear_charges, self.coordinates, size = size, indices = indices,
206214
sorting = sorting, central_cutoff = central_cutoff, central_decay = central_decay,
207215
interaction_cutoff = interaction_cutoff, interaction_decay = interaction_decay)
208216

@@ -284,7 +292,6 @@ def generate_slatm(self, mbtypes,
284292
if local: slatm = np.asarray(slatm)
285293
self.representation = slatm
286294

287-
288295
def read_xyz(self, filename):
289296
"""(Re-)initializes the Compound-object with data from an xyz-file.
290297

qml/frepresentations.f90

+91-74
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
! MIT License
22
!
3-
! Copyright (c) 2016 Anders Steen Christensen
3+
! Copyright (c) 2016-2017 Anders Steen Christensen, Lars Andersen Bratholm
44
!
55
! Permission is hereby granted, free of charge, to any person obtaining a copy
66
! of this software and associated documentation files (the "Software"), to deal
@@ -204,18 +204,21 @@ subroutine fgenerate_unsorted_coulomb_matrix(atomic_charges, coordinates, nmax,
204204

205205
end subroutine fgenerate_unsorted_coulomb_matrix
206206

207-
subroutine fgenerate_local_coulomb_matrix(atomic_charges, coordinates, natoms, nmax, &
208-
& cent_cutoff, cent_decay, int_cutoff, int_decay, cm)
207+
subroutine fgenerate_local_coulomb_matrix(central_atom_indices, central_natoms, &
208+
& atomic_charges, coordinates, natoms, nmax, cent_cutoff, cent_decay, &
209+
& int_cutoff, int_decay, cm)
209210

210211
implicit none
211212

213+
integer, intent(in) :: central_natoms
214+
integer, dimension(:), intent(in) :: central_atom_indices
212215
double precision, dimension(:), intent(in) :: atomic_charges
213216
double precision, dimension(:,:), intent(in) :: coordinates
214217
integer,intent(in) :: natoms
215218
integer, intent(in) :: nmax
216219
double precision, intent(inout) :: cent_cutoff, cent_decay, int_cutoff, int_decay
217220

218-
double precision, dimension(natoms, ((nmax + 1) * nmax) / 2), intent(out):: cm
221+
double precision, dimension(central_natoms, ((nmax + 1) * nmax) / 2), intent(out):: cm
219222

220223
integer :: idx
221224

@@ -232,10 +235,11 @@ subroutine fgenerate_local_coulomb_matrix(atomic_charges, coordinates, natoms, n
232235
double precision, allocatable, dimension(:, :) :: distance_matrix
233236
double precision, allocatable, dimension(:, :) :: distance_matrix_tmp
234237

235-
integer i, j, m, n, k
238+
integer i, j, m, n, k, l
236239

237240
double precision, parameter :: pi = 4.0d0 * atan(1.0d0)
238241

242+
239243
if (size(coordinates, dim=1) /= size(atomic_charges, dim=1)) then
240244
write(*,*) "ERROR: Coulomb matrix generation"
241245
write(*,*) size(coordinates, dim=1), "coordinates, but", &
@@ -287,25 +291,28 @@ subroutine fgenerate_local_coulomb_matrix(atomic_charges, coordinates, natoms, n
287291
enddo
288292
!$OMP END PARALLEL DO
289293

290-
do i = 1, natoms
291-
if (cutoff_count(i) > nmax) then
294+
do i = 1, central_natoms
295+
j = central_atom_indices(i)
296+
if (cutoff_count(j) > nmax) then
292297
write(*,*) "ERROR: Coulomb matrix generation"
293298
write(*,*) nmax, "size set, but", &
294-
& cutoff_count(i), "size needed!"
299+
& cutoff_count(j), "size needed!"
295300
stop
296301
endif
297302
enddo
298303

299304
! Allocate temporary
300-
allocate(pair_distance_matrix(natoms, natoms, natoms))
301-
allocate(row_norms(natoms, natoms))
305+
allocate(pair_distance_matrix(natoms, natoms, central_natoms))
306+
allocate(row_norms(natoms, central_natoms))
302307

303308
pair_distance_matrix = 0.0d0
304309
row_norms = 0.0d0
305310

306-
!$OMP PARALLEL DO PRIVATE(pair_norm, prefactor) REDUCTION(+:row_norms) COLLAPSE(2)
311+
312+
!$OMP PARALLEL DO PRIVATE(pair_norm, prefactor, k) REDUCTION(+:row_norms) COLLAPSE(2)
307313
do i = 1, natoms
308-
do k = 1, natoms
314+
do l = 1, central_natoms
315+
k = central_atom_indices(l)
309316
! self interaction
310317
if (distance_matrix(i,k) > cent_cutoff) then
311318
cycle
@@ -318,8 +325,8 @@ subroutine fgenerate_local_coulomb_matrix(atomic_charges, coordinates, natoms, n
318325
endif
319326

320327
pair_norm = prefactor * prefactor * 0.5d0 * atomic_charges(i) ** 2.4d0
321-
pair_distance_matrix(i,i,k) = pair_norm
322-
row_norms(i,k) = row_norms(i,k) + pair_norm * pair_norm
328+
pair_distance_matrix(i,i,l) = pair_norm
329+
row_norms(i,l) = row_norms(i,l) + pair_norm * pair_norm
323330

324331
do j = i+1, natoms
325332
if (distance_matrix(j,k) > cent_cutoff) then
@@ -344,32 +351,34 @@ subroutine fgenerate_local_coulomb_matrix(atomic_charges, coordinates, natoms, n
344351
& * (distance_matrix(j,k) - cent_cutoff + cent_decay) / cent_decay) + 1)
345352
endif
346353

347-
pair_distance_matrix(i, j, k) = pair_norm
348-
pair_distance_matrix(j, i, k) = pair_norm
354+
pair_distance_matrix(i, j, l) = pair_norm
355+
pair_distance_matrix(j, i, l) = pair_norm
349356
pair_norm = pair_norm * pair_norm
350-
row_norms(i,k) = row_norms(i,k) + pair_norm
351-
row_norms(j,k) = row_norms(j,k) + pair_norm
357+
row_norms(i,l) = row_norms(i,l) + pair_norm
358+
row_norms(j,l) = row_norms(j,l) + pair_norm
352359
enddo
353360
enddo
354361
enddo
355362
!$OMP END PARALLEL DO
356363

357364
! Allocate temporary
358-
allocate(sorted_atoms_all(natoms, natoms))
365+
allocate(sorted_atoms_all(natoms, central_natoms))
359366

360-
!$OMP PARALLEL DO
361-
do k = 1, natoms
362-
row_norms(k,k) = huge_double
367+
!$OMP PARALLEL DO PRIVATE(k)
368+
do l = 1, central_natoms
369+
k = central_atom_indices(l)
370+
row_norms(k,l) = huge_double
363371
enddo
364372
!$OMP END PARALLEL DO
365373

366-
!$OMP PARALLEL DO PRIVATE(j)
367-
do k = 1, natoms
374+
!$OMP PARALLEL DO PRIVATE(j,k)
375+
do l = 1, central_natoms
376+
k = central_atom_indices(l)
368377
!$OMP CRITICAL
369378
do i = 1, cutoff_count(k)
370-
j = maxloc(row_norms(:,k), dim=1)
371-
sorted_atoms_all(i, k) = j
372-
row_norms(j,k) = 0.0d0
379+
j = maxloc(row_norms(:,l), dim=1)
380+
sorted_atoms_all(i, l) = j
381+
row_norms(j,l) = 0.0d0
373382
enddo
374383
!$OMP END CRITICAL
375384
enddo
@@ -383,14 +392,15 @@ subroutine fgenerate_local_coulomb_matrix(atomic_charges, coordinates, natoms, n
383392
! Fill coulomb matrix according to sorted row-2-norms
384393
cm = 0.0d0
385394

386-
!$OMP PARALLEL DO PRIVATE(i, j, idx)
387-
do k = 1, natoms
395+
!$OMP PARALLEL DO PRIVATE(i, j, k, idx)
396+
do l = 1, central_natoms
397+
k = central_atom_indices(l)
388398
do m = 1, cutoff_count(k)
389-
i = sorted_atoms_all(m, k)
399+
i = sorted_atoms_all(m, l)
390400
idx = (m*m+m)/2 - m
391401
do n = 1, m
392-
j = sorted_atoms_all(n, k)
393-
cm(k, idx+n) = pair_distance_matrix(i,j,k)
402+
j = sorted_atoms_all(n, l)
403+
cm(l, idx+n) = pair_distance_matrix(i,j,l)
394404
enddo
395405
enddo
396406
enddo
@@ -403,18 +413,20 @@ subroutine fgenerate_local_coulomb_matrix(atomic_charges, coordinates, natoms, n
403413

404414
end subroutine fgenerate_local_coulomb_matrix
405415

406-
subroutine fgenerate_atomic_coulomb_matrix(atomic_charges, coordinates, natoms, nmax, &
407-
& cent_cutoff, cent_decay, int_cutoff, int_decay, cm)
416+
subroutine fgenerate_atomic_coulomb_matrix(central_atom_indices, central_natoms, atomic_charges, &
417+
& coordinates, natoms, nmax, cent_cutoff, cent_decay, int_cutoff, int_decay, cm)
408418

409419
implicit none
410420

421+
integer, dimension(:), intent(in) :: central_atom_indices
422+
integer, intent(in) :: central_natoms
411423
double precision, dimension(:), intent(in) :: atomic_charges
412424
double precision, dimension(:,:), intent(in) :: coordinates
413425
integer,intent(in) :: natoms
414426
integer, intent(in) :: nmax
415427
double precision, intent(inout) :: cent_cutoff, cent_decay, int_cutoff, int_decay
416428

417-
double precision, dimension(natoms, ((nmax + 1) * nmax) / 2), intent(out):: cm
429+
double precision, dimension(central_natoms, ((nmax + 1) * nmax) / 2), intent(out):: cm
418430

419431
integer :: idx
420432

@@ -430,7 +442,7 @@ subroutine fgenerate_atomic_coulomb_matrix(atomic_charges, coordinates, natoms,
430442
double precision, allocatable, dimension(:, :) :: distance_matrix
431443
double precision, allocatable, dimension(:, :) :: distance_matrix_tmp
432444

433-
integer i, j, m, n, k
445+
integer i, j, m, n, k, l
434446

435447
double precision, parameter :: pi = 4.0d0 * atan(1.0d0)
436448

@@ -485,11 +497,12 @@ subroutine fgenerate_atomic_coulomb_matrix(atomic_charges, coordinates, natoms,
485497
enddo
486498
!$OMP END PARALLEL DO
487499

488-
do i = 1, natoms
489-
if (cutoff_count(i) > nmax) then
500+
do i = 1, central_natoms
501+
k = central_atom_indices(i)
502+
if (cutoff_count(k) > nmax) then
490503
write(*,*) "ERROR: Coulomb matrix generation"
491504
write(*,*) nmax, "size set, but", &
492-
& cutoff_count(i), "size needed!"
505+
& cutoff_count(k), "size needed!"
493506
stop
494507
endif
495508
enddo
@@ -520,20 +533,22 @@ subroutine fgenerate_atomic_coulomb_matrix(atomic_charges, coordinates, natoms,
520533
!$OMP END PARALLEL DO
521534

522535
! Allocate temporary
523-
allocate(sorted_atoms_all(natoms, natoms))
536+
allocate(distance_matrix_tmp(natoms, natoms))
537+
allocate(sorted_atoms_all(natoms, central_natoms))
524538

525539
distance_matrix_tmp = distance_matrix
526540
!Generate sorted list of atom ids by distance matrix
527-
!$OMP PARALLEL DO PRIVATE(j)
528-
do k = 1, natoms
529-
!$OMP CRITICAL
530-
do i = 1, cutoff_count(k)
531-
j = minloc(distance_matrix_tmp(:,k), dim=1)
532-
sorted_atoms_all(i, k) = j
533-
distance_matrix_tmp(j, k) = huge_double
534-
enddo
535-
!$OMP END CRITICAL
536-
enddo
541+
!$OMP PARALLEL DO PRIVATE(j, k)
542+
do l = 1, central_natoms
543+
k = central_atom_indices(l)
544+
!$OMP CRITICAL
545+
do i = 1, cutoff_count(k)
546+
j = minloc(distance_matrix_tmp(:,k), dim=1)
547+
sorted_atoms_all(i, l) = j
548+
distance_matrix_tmp(j, k) = huge_double
549+
enddo
550+
!$OMP END CRITICAL
551+
enddo
537552
!$OMP END PARALLEL DO
538553

539554
! Clean up
@@ -544,34 +559,36 @@ subroutine fgenerate_atomic_coulomb_matrix(atomic_charges, coordinates, natoms,
544559

545560
pair_norm = 0.0d0
546561

547-
do k = 1, natoms
548-
do m = 1, cutoff_count(k)
549-
i = sorted_atoms_all(m, k)
562+
!$OMP PARALLEL DO PRIVATE(i, prefactor, idx, j, pair_norm, k)
563+
do l = 1, central_natoms
564+
k = central_atom_indices(l)
565+
do m = 1, cutoff_count(k)
566+
i = sorted_atoms_all(m, l)
550567

551-
if (distance_matrix(i,k) > cent_cutoff) then
552-
cycle
553-
endif
554-
prefactor = 1.0d0
555-
if (distance_matrix(i,k) > cent_cutoff - cent_decay) then
556-
prefactor = 0.5d0 * (cos(pi &
557-
& * (distance_matrix(i,k) - cent_cutoff + cent_decay) &
558-
& / cent_decay) + 1.0d0)
559-
endif
568+
if (distance_matrix(i,k) > cent_cutoff) then
569+
cycle
570+
endif
571+
prefactor = 1.0d0
572+
if (distance_matrix(i,k) > cent_cutoff - cent_decay) then
573+
prefactor = 0.5d0 * (cos(pi &
574+
& * (distance_matrix(i,k) - cent_cutoff + cent_decay) &
575+
& / cent_decay) + 1.0d0)
576+
endif
560577

561-
idx = (m*m+m)/2 - m
562-
do n = 1, m
563-
j = sorted_atoms_all(n, k)
564-
565-
pair_norm = prefactor * pair_distance_matrix(i, j)
566-
if (distance_matrix(j,k) > cent_cutoff - cent_decay) then
567-
pair_norm = pair_norm * 0.5d0 * (cos(pi &
568-
& * (distance_matrix(j,k) - cent_cutoff + cent_decay) &
569-
& / cent_decay) + 1)
570-
endif
571-
cm(k, idx+n) = pair_norm
572-
enddo
578+
idx = (m*m+m)/2 - m
579+
do n = 1, m
580+
j = sorted_atoms_all(n, l)
581+
582+
pair_norm = prefactor * pair_distance_matrix(i, j)
583+
if (distance_matrix(j,k) > cent_cutoff - cent_decay) then
584+
pair_norm = pair_norm * 0.5d0 * (cos(pi &
585+
& * (distance_matrix(j,k) - cent_cutoff + cent_decay) &
586+
& / cent_decay) + 1)
587+
endif
588+
cm(l, idx+n) = pair_norm
573589
enddo
574590
enddo
591+
enddo
575592

576593
! Clean up
577594
deallocate(distance_matrix)

0 commit comments

Comments
 (0)