forked from libAtoms/GAP
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathgap_fit_module.f95
2336 lines (1893 loc) · 106 KB
/
gap_fit_module.f95
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
! HND XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
! HND X
! HND X GAP (Gaussian Approximation Potental)
! HND X
! HND X
! HND X Portions of GAP were written by Albert Bartok-Partay, Gabor Csanyi,
! HND X and Sascha Klawohn. Copyright 2006-2021.
! HND X
! HND X Portions of GAP were written by Noam Bernstein as part of
! HND X his employment for the U.S. Government, and are not subject
! HND X to copyright in the USA.
! HND X
! HND X GAP is published and distributed under the
! HND X Academic Software License v1.0 (ASL)
! HND X
! HND X GAP is distributed in the hope that it will be useful for non-commercial
! HND X academic research, but WITHOUT ANY WARRANTY; without even the implied
! HND X warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! HND X ASL for more details.
! HND X
! HND X You should have received a copy of the ASL along with this program
! HND X (e.g. in a LICENSE.md file); if not, you can write to the original licensors,
! HND X Gabor Csanyi or Albert Bartok-Partay. The ASL is also published at
! HND X http://github.com/gabor1/ASL
! HND X
! HND X When using this software, please cite the following reference:
! HND X
! HND X A. P. Bartok et al Physical Review Letters vol 104 p136403 (2010)
! HND X
! HND X When using the SOAP kernel or its variants, please additionally cite:
! HND X
! HND X A. P. Bartok et al Physical Review B vol 87 p184115 (2013)
! HND X
! HND XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
#include "error.inc"
module gap_fit_module
use error_module
use libatoms_module
use descriptors_module
use gp_predict_module
use gp_fit_module
use fox_wxml
use potential_module
use ScaLAPACK_module
use task_manager_module
use MPI_context_module, only : is_root
implicit none
integer, parameter :: SPARSE_LENGTH = 10000
integer, parameter :: THETA_LENGTH = 10000
integer, parameter :: E0_ISOLATED = 1
integer, parameter :: E0_AVERAGE = 2
integer, parameter :: EXCLUDE_LOC = -1
#ifdef GAP_VERSION
integer, parameter, private :: gap_version = GAP_VERSION
#else
integer, parameter, private :: gap_version = huge(1)
#endif
type gap_fit
!% everything from the command line
type(Atoms), dimension(:), allocatable :: at
character(len=STRING_LENGTH) :: at_file='', core_ip_args = '', e0_str, local_property0_str, &
energy_parameter_name, local_property_parameter_name, force_parameter_name, virial_parameter_name, &
stress_parameter_name, hessian_parameter_name, config_type_parameter_name, sigma_parameter_name, &
config_type_sigma_string, core_param_file, gp_file, template_file, force_mask_parameter_name, &
condition_number_norm, linear_system_dump_file, config_file
character(len=10240) :: command_line = ''
real(dp), dimension(total_elements) :: e0, local_property0
real(dp) :: max_cutoff
real(dp), dimension(4) :: default_sigma
real(dp) :: default_local_property_sigma
real(dp) :: sparse_jitter, e0_offset, hessian_delta
integer :: e0_method = E0_ISOLATED
logical :: do_core = .false., do_copy_at_file, has_config_type_sigma, sigma_per_atom = .true.
logical :: sparsify_only_no_fit = .false.
logical :: dryrun = .false.
integer :: n_frame = 0
integer :: n_coordinate = 0
integer :: n_ener = 0
integer :: n_force = 0
integer :: n_virial = 0
integer :: n_hessian = 0
integer :: n_local_property = 0
integer :: n_species = 0
integer :: min_save
integer :: mpi_blocksize_rows = 0
integer :: mpi_blocksize_cols = 0
type(extendable_str) :: quip_string, config_string
type(Potential) :: core_pot
type(gpFull) :: my_gp
type(gpSparse) :: gp_sp
type(MPI_Context) :: mpi_obj
type(ScaLAPACK) :: ScaLAPACK_obj
type(task_manager_type) :: task_manager
type(descriptor), dimension(:), allocatable :: my_descriptor
character(len=STRING_LENGTH), dimension(200) :: gap_str
real(dp), dimension(:), allocatable :: delta, f0, theta_uniform, zeta, unique_hash_tolerance, unique_descriptor_tolerance !, theta
real(dp), dimension(:,:), allocatable :: sigma
integer, dimension(:), allocatable :: n_sparseX, sparse_method, target_type, n_cross, n_descriptors, species_Z, covariance_type
integer, dimension(:,:), allocatable :: config_type_n_sparseX
character(len=STRING_LENGTH), dimension(:), allocatable :: theta_file, sparse_file, theta_fac_string, config_type, config_type_n_sparseX_string, print_sparse_index
logical, dimension(:), allocatable :: mark_sparse_atoms, add_species, has_theta_fac, has_theta_uniform, has_theta_file, has_zeta
logical :: sparseX_separate_file
logical :: sparse_use_actual_gpcov
logical :: has_template_file, has_e0, has_local_property0, has_e0_offset, has_linear_system_dump_file, has_config_file
endtype gap_fit
private
public :: fit_n_from_xyz
public :: fit_data_from_xyz
public :: e0_from_xyz
public :: w_Z_from_xyz
public :: gap_fit
public :: gap_fit_print_xml
public :: file_print_xml
! public :: print_sparse
public :: set_baselines
public :: parse_config_type_sigma
public :: parse_config_type_n_sparseX
public :: read_fit_xyz
public :: read_descriptors
public :: get_species_xyz
public :: add_multispecies_gaps
public :: add_template_string
public :: gap_fit_parse_command_line
public :: gap_fit_parse_gap_str
public :: gap_fit_read_core_param_file
public :: gap_fit_init_mpi_scalapack
public :: gap_fit_init_task_manager
public :: gap_fit_distribute_tasks
public :: gap_fit_set_mpi_blocksizes
public :: gap_fit_is_root
public :: gap_fit_print_linear_system_dump_file
public :: gap_fit_estimate_memory
contains
subroutine gap_fit_parse_command_line(this)
!% This subroutine parses the main command line options.
type(gap_fit), intent(inout), target :: this
type(Dictionary) :: params
character(len=STRING_LENGTH), pointer :: at_file, e0_str, local_property0_str, &
core_param_file, core_ip_args, &
energy_parameter_name, local_property_parameter_name, force_parameter_name, &
virial_parameter_name, stress_parameter_name, hessian_parameter_name, &
config_type_parameter_name, sigma_parameter_name, config_type_sigma_string, &
gp_file, template_file, force_mask_parameter_name, condition_number_norm, &
linear_system_dump_file, config_file
character(len=STRING_LENGTH) :: gap_str, verbosity, sparse_method_str, covariance_type_str, e0_method, &
parameter_name_prefix
logical, pointer :: sigma_per_atom, do_copy_at_file, sparseX_separate_file, sparse_use_actual_gpcov, sparsify_only_no_fit
logical, pointer :: dryrun
logical :: do_ip_timing, has_sparse_file, has_theta_uniform, has_at_file, has_gap, has_config_file, has_default_sigma
logical :: mpi_print_all, file_exists
real(dp), pointer :: e0_offset, sparse_jitter, hessian_delta
real(dp), dimension(:), pointer :: default_sigma
real(dp), pointer :: default_local_property_sigma
integer :: rnd_seed
integer, pointer :: mpi_blocksize_rows, mpi_blocksize_cols
config_file => this%config_file
at_file => this%at_file
e0_str => this%e0_str
local_property0_str => this%local_property0_str
e0_offset => this%e0_offset
default_sigma => this%default_sigma
default_local_property_sigma => this%default_local_property_sigma
sparse_jitter => this%sparse_jitter
hessian_delta => this%hessian_delta
core_param_file => this%core_param_file
core_ip_args => this%core_ip_args
energy_parameter_name => this%energy_parameter_name
local_property_parameter_name => this%local_property_parameter_name
force_parameter_name => this%force_parameter_name
virial_parameter_name => this%virial_parameter_name
stress_parameter_name => this%stress_parameter_name
hessian_parameter_name => this%hessian_parameter_name
config_type_parameter_name => this%config_type_parameter_name
sigma_parameter_name => this%sigma_parameter_name
force_mask_parameter_name => this%force_mask_parameter_name
config_type_sigma_string => this%config_type_sigma_string
sigma_per_atom => this%sigma_per_atom
do_copy_at_file => this%do_copy_at_file
sparseX_separate_file => this%sparseX_separate_file
sparse_use_actual_gpcov => this%sparse_use_actual_gpcov
gp_file => this%gp_file
template_file => this%template_file
sparsify_only_no_fit => this%sparsify_only_no_fit
dryrun => this%dryrun
condition_number_norm => this%condition_number_norm
linear_system_dump_file => this%linear_system_dump_file
mpi_blocksize_rows => this%mpi_blocksize_rows
mpi_blocksize_cols => this%mpi_blocksize_cols
call initialise(params)
call param_register(params, 'config_file', '', config_file, has_value_target=has_config_file, &
help_string="File as alternative input (newlines converted to spaces)")
! check if config file is given, ignore everything else
! prepare parsing of config file or command line string later
if (param_read_args(params, ignore_unknown=.true., command_line=this%command_line)) then
if (has_config_file) then
inquire(file=config_file, exist=file_exists)
if (.not. file_exists) call system_abort("Config file does not exist: "//config_file)
call read(this%config_string, config_file, keep_lf=.true., mpi_comm=this%mpi_obj%communicator, mpi_id=this%mpi_obj%my_proc)
end if
end if
if (.not. has_config_file) this%config_string = this%command_line
this%has_config_file = has_config_file
call param_register(params, 'atoms_filename', '//MANDATORY//', at_file, has_value_target = has_at_file, help_string="XYZ file with fitting configurations", altkey="at_file")
call param_register(params, 'gap', '//MANDATORY//', gap_str, has_value_target = has_gap, help_string="Initialisation string for GAPs")
call param_register(params, 'e0', '0.0', e0_str, has_value_target = this%has_e0, &
help_string="Atomic energy value to be subtracted from energies before fitting (and added back on after prediction). &
& Specifiy a single number (used for all species) or by species: {Ti:-150.0:O:-320...}. energy = baseline + GAP + e0")
call param_register(params, 'local_property0', '0.0', local_property0_str, has_value_target = this%has_local_property0, &
help_string="Local property value to be subtracted from the local property before fitting (and added back on after prediction). &
& Specifiy a single number (used for all species) or by species: {H:20.0:Cl:35.0...}.")
call param_register(params, 'e0_offset', '0.0', e0_offset, has_value_target = this%has_e0_offset, &
help_string="Offset of baseline. If zero, the offset is the average atomic energy of the input data or the e0 specified manually.")
call param_register(params, 'e0_method','isolated',e0_method, &
help_string="Method to determine e0, if not explicitly specified. Possible options: isolated (default, each atom &
present in the XYZ needs to have an isolated representative, with a valid energy), average (e0 is the average of &
all total energies across the XYZ)")
call param_register(params, 'default_kernel_regularisation', '//MANDATORY//', default_sigma, has_value_target = has_default_sigma, &
help_string="error in [energies forces virials hessians]", altkey="default_sigma")
call param_register(params, 'default_kernel_regularisation_local_property', '0.001', default_local_property_sigma, &
help_string="error in local_property", altkey="default_local_property_sigma")
call param_register(params, 'sparse_jitter', "1.0e-10", sparse_jitter, &
help_string="Extra regulariser used to regularise the sparse covariance matrix before it is passed to the linear solver. Use something small, it really shouldn't affect your results, if it does, your sparse basis is still very ill-conditioned.")
call param_register(params, 'hessian_displacement', "1.0e-2", hessian_delta, &
help_string="Finite displacement to use in numerical differentiation when obtaining second derivative for the Hessian covariance", altkey="hessian_delta")
call param_register(params, 'baseline_param_filename', 'quip_params.xml', core_param_file, &
help_string="QUIP XML file which contains a potential to subtract from data (and added back after prediction)", altkey="core_param_file")
call param_register(params, 'baseline_ip_args', '', core_ip_args, has_value_target = this%do_core, &
help_string=" QUIP init string for a potential to subtract from data (and added back after prediction)", altkey="core_ip_args")
call param_register(params, 'energy_parameter_name', 'energy', energy_parameter_name, &
help_string="Name of energy property in the input XYZ file that describes the data")
call param_register(params, 'local_property_parameter_name', 'local_property', local_property_parameter_name, &
help_string="Name of local_property (column) in the input XYZ file that describes the data")
call param_register(params, 'force_parameter_name', 'force', force_parameter_name, &
help_string="Name of force property (columns) in the input XYZ file that describes the data")
call param_register(params, 'virial_parameter_name', 'virial', virial_parameter_name, &
help_string="Name of virial property in the input XYZ file that describes the data")
call param_register(params, 'stress_parameter_name', 'stress', stress_parameter_name, &
help_string="Name of stress property (6-vector or 9-vector) in the input XYZ file that describes the data - stress values only used if virials are not available (opposite sign, standard Voigt order)")
call param_register(params, 'hessian_parameter_name', 'hessian', hessian_parameter_name, &
help_string="Name of hessian property (column) in the input XYZ file that describes the data")
call param_register(params, 'config_type_parameter_name', 'config_type', config_type_parameter_name, &
help_string="Allows grouping on configurations into. This option is the name of the key that indicates the configuration type in the input XYZ file. With the default, the key-value pair config_type=blah would place that configuration into the group blah.")
call param_register(params, 'kernel_regularisation_parameter_name', 'sigma', sigma_parameter_name, &
help_string="kernel regularisation parameters for a given configuration in the database. &
Overrides the command line values (both defaults and config-type-specific values). In the input XYZ file, it must be prepended by energy_, force_, virial_ or hessian_", altkey="sigma_parameter_name")
call param_register(params, 'force_mask_parameter_name', 'force_mask', force_mask_parameter_name, &
help_string="To exclude forces on specific atoms from the fit. In the XYZ, it must be a logical column.")
call param_register(params, 'parameter_name_prefix', '', parameter_name_prefix, &
help_string="Prefix that gets uniformly appended in front of {energy,local_property,force,virial,...}_parameter_name")
call param_register(params, 'config_type_kernel_regularisation', '', config_type_sigma_string, has_value_target = this%has_config_type_sigma, &
help_string="What kernel regularisation values to choose for each type of data, when the configurations are grouped into config_types. Format: {configtype1:energy:force:virial:hessian:config_type2:energy:force:virial:hessian...}", altkey="config_type_sigma")
call param_register(params, 'kernel_regularisation_is_per_atom', 'T', sigma_per_atom, &
help_string="Interpretation of the energy and virial sigmas specified in >>default_kernel_regularisation<< and >>config_type_kernel_regularisation<<. &
If >>T<<, they are interpreted as per-atom errors, and the variance will be scaled according to the number of atoms in the configuration. &
If >>F<< they are treated as absolute errors and no scaling is performed. &
NOTE: values specified on a per-configuration basis (see >>kernel_regularisation_parameter_name<<) are always absolute, not per-atom.", altkey="sigma_per_atom")
call param_register(params, 'do_copy_atoms_file', 'T', do_copy_at_file, &
help_string="Copy the input XYZ file into the GAP XML file (should be set to False for NetCDF input).", altkey="do_copy_at_file")
call param_register(params, 'sparse_separate_file', 'T', sparseX_separate_file, &
help_string="Save sparse point data in separate file in binary (use it for large datasets)")
call param_register(params, 'sparse_use_actual_gpcov', 'F', sparse_use_actual_gpcov, &
help_string="Use actual GP covariance for sparsification methods")
call param_register(params, 'gap_file', 'gap_new.xml', gp_file, &
help_string="Name of output XML file that will contain the fitted potential", altkey="gp_file")
call param_register(params, 'verbosity', 'NORMAL', verbosity, &
help_string="Verbosity control. Options: NORMAL, VERBOSE, NERD, ANALYSIS.") ! changed name to ANALYSIS now that we are grown up
call param_register(params, "rnd_seed", "-1", rnd_seed, &
help_string="Random seed.")
call param_register(params, "openmp_chunk_size", "0", openmp_chunk_size, &
help_string="Chunk size in OpenMP scheduling; 0: each thread gets a single block of similar size (default)")
call param_register(params, 'do_ip_timing', 'F', do_ip_timing, &
help_string="To enable or not timing of the interatomic potential.")
call param_register(params, 'template_file', 'template.xyz', template_file, has_value_target=this%has_template_file, &
help_string="Template XYZ file for initialising object")
call param_register(params, 'sparsify_only_no_fit', 'F', sparsify_only_no_fit, &
help_string="If true, sparsification is done, but no fitting. print the sparse index by adding print_sparse_index=file.dat to the descriptor string.")
call param_register(params, 'dryrun', 'F', dryrun, &
help_string="If true, exits after memory estimate, before major allocations.")
call param_register(params, 'condition_number_norm', ' ', condition_number_norm, &
help_string="Norm for condition number of matrix A; O: 1-norm, I: inf-norm, <space>: skip calculation (default)")
call param_register(params, 'linear_system_dump_file', '', linear_system_dump_file, has_value_target=this%has_linear_system_dump_file, &
help_string="Basename prefix of linear system dump files. Skipped if empty (default).")
call param_register(params, 'mpi_blocksize_rows', '0', mpi_blocksize_rows, &
help_string="Blocksize of MPI distributed matrix rows. Affects efficiency and memory usage slightly. Max if 0 (default).")
call param_register(params, 'mpi_blocksize_cols', '100', mpi_blocksize_cols, &
help_string="Blocksize of MPI distributed matrix cols. Affects efficiency and memory usage considerably. Max if 0. Default: 100")
call param_register(params, 'mpi_print_all', 'F', mpi_print_all, &
help_string="If true, each MPI processes will print its output. Otherwise, only the first process does (default).")
if (.not. param_read_line(params, replace(string(this%config_string), quip_new_line, ' '))) then
call system_abort('Exit: Mandatory argument(s) missing...')
endif
call print_title("Input parameters")
call param_print(params)
call print_title("")
call finalise(params)
if (mpi_print_all) then
call mpi_all_inoutput(mainlog, .true.)
call activate(mainlog)
call mpi_all_inoutput(errorlog, .true.)
call activate(errorlog)
end if
if (len_trim(parameter_name_prefix) > 0) then
energy_parameter_name = trim(parameter_name_prefix) // trim(energy_parameter_name)
local_property_parameter_name = trim(parameter_name_prefix) // trim(local_property_parameter_name)
force_parameter_name = trim(parameter_name_prefix) // trim(force_parameter_name)
virial_parameter_name = trim(parameter_name_prefix) // trim(virial_parameter_name)
hessian_parameter_name = trim(parameter_name_prefix) // trim(hessian_parameter_name)
stress_parameter_name = trim(parameter_name_prefix) // trim(stress_parameter_name)
config_type_parameter_name = trim(parameter_name_prefix) // trim(config_type_parameter_name)
sigma_parameter_name = trim(parameter_name_prefix) // trim(sigma_parameter_name)
force_mask_parameter_name = trim(parameter_name_prefix) // trim(force_mask_parameter_name)
endif
if (sparsify_only_no_fit) then
force_parameter_name = '//IGNORE//'
virial_parameter_name = '//IGNORE//'
hessian_parameter_name = '//IGNORE//'
stress_parameter_name = '//IGNORE//'
call print_warning("sparsify_only_no_fit == T: force, virial, hessian, stress parameters are ignored.")
end if
if( len_trim(this%gp_file) > 216 ) then ! The filename's length is limited to 255 char.s in some filesystem.
! Without this check, the fit would run but produce a core file and only a temporary xml file.
! The limit is set to 216 as the sparse file can be 39 characters longer.
call system_abort("gap_file's name "//this%gp_file//" is too long. Please start the fit again with a shorter name.")
endif
if(do_ip_timing) call enable_timing()
select case(verbosity)
case ("NORMAL")
call verbosity_push(PRINT_NORMAL)
case ("VERBOSE")
call verbosity_push(PRINT_VERBOSE)
case ("NERD")
call verbosity_push(PRINT_NERD)
case ("ANALYSIS") ! changed name now that we are grown up
call verbosity_push(PRINT_ANALYSIS) ! changed name now that we are grown up
case default
call system_abort("confused by verbosity " // trim(verbosity))
end select
select case(lower_case(e0_method))
case ("isolated")
this%e0_method = E0_ISOLATED
case ("average")
this%e0_method = E0_AVERAGE
case default
call system_abort("confused by e0_method " // trim(e0_method))
end select
if (rnd_seed >= 0) call system_set_random_seeds(rnd_seed)
call print_title('Gaussian Approximation Potentials - Database fitting')
call print('')
call print('Initial parsing of command line arguments finished.')
call split_string(gap_str,':;','{}',this%gap_str(:),this%n_coordinate,matching=.true.)
call print('Found '//this%n_coordinate//' GAPs.')
endsubroutine gap_fit_parse_command_line
subroutine set_baselines(this)
type(gap_fit), intent(inout) :: this
integer :: i
this%e0 = 0.0_dp
if( count( (/this%has_e0, this%has_e0_offset/) ) > 1 ) then
call print_warning('Both e0 and e0_offset has been specified. That means your atomic energy is e0 + e0_offset')
endif
if( this%has_e0 ) then
call parse_atomtype_value_str(this%e0_str,this%e0)
else
call e0_from_xyz(this) ! calculates the average atomic energy so it can be subtracted later.
endif
if( this%has_e0_offset ) this%e0 = this%e0 + this%e0_offset
if( .not. this%has_e0 ) then
do i = 1, size(this%e0)
if( all(i/=this%species_Z) ) this%e0(i) = 0.0_dp
enddo
call print('E0/atom = '//this%e0)
endif
if( this%has_local_property0 ) then
call parse_atomtype_value_str(this%local_property0_str,this%local_property0)
this%e0 = 0.0_dp
else
this%local_property0 = 0.0_dp
endif
endsubroutine set_baselines
subroutine parse_atomtype_value_str(this,values,error)
character(len=STRING_LENGTH), intent(in) :: this
real(dp), dimension(total_elements), intent(out) :: values
integer, intent(out), optional :: error
integer :: n_string_array, i, z
character(len=STRING_LENGTH), dimension(2*total_elements) :: string_array
INIT_ERROR(error)
call split_string(this,':','{}',string_array(:),n_string_array,matching=.true.)
if(n_string_array == 1) then
values = string_to_real(trim(string_array(1)))
elseif(mod(n_string_array,2) == 0) then
values = 0.0_dp
do i = 1, n_string_array / 2
z = atomic_number(trim( string_array((i-1)*2+1) ))
if( z==0 ) then
RAISE_ERROR("parse_atomtype_value_str: invalid atomic symbol "//trim(string_array((i-1)*2+1)),error)
endif
values(z) = string_to_real(trim( string_array(2*i) ))
enddo
else
RAISE_ERROR("parse_atomtype_value_str: number of fields is an odd number. It must be a list of pairs of values, such as {Ti:-150.4:O:-345.1}",error)
endif
endsubroutine parse_atomtype_value_str
subroutine gap_fit_parse_gap_str(this)
!% This subroutine parses the options given in the gap string, for each GAP.
type(gap_fit), intent(inout), target :: this
type(Dictionary) :: params
integer :: i_coordinate
real(dp) :: delta, f0, theta_uniform, zeta, unique_hash_tolerance, unique_descriptor_tolerance
integer :: n_sparseX, sparse_method, covariance_type
character(len=STRING_LENGTH) :: config_type_n_sparseX_string, theta_fac_string, theta_file, sparse_file, print_sparse_index, &
covariance_type_str, sparse_method_str
logical :: mark_sparse_atoms, add_species, has_sparse_file
allocate(this%delta(this%n_coordinate))
allocate(this%f0(this%n_coordinate))
allocate(this%n_sparseX(this%n_coordinate))
allocate(this%config_type_n_sparseX_string(this%n_coordinate))
allocate(this%theta_fac_string(this%n_coordinate))
allocate(this%theta_uniform(this%n_coordinate))
allocate(this%theta_file(this%n_coordinate))
allocate(this%has_theta_fac(this%n_coordinate))
allocate(this%has_theta_uniform(this%n_coordinate))
allocate(this%has_theta_file(this%n_coordinate))
allocate(this%sparse_file(this%n_coordinate))
allocate(this%mark_sparse_atoms(this%n_coordinate))
allocate(this%sparse_method(this%n_coordinate))
allocate(this%add_species(this%n_coordinate))
allocate(this%covariance_type(this%n_coordinate))
allocate(this%zeta(this%n_coordinate))
allocate(this%has_zeta(this%n_coordinate))
allocate(this%print_sparse_index(this%n_coordinate))
allocate(this%unique_hash_tolerance(this%n_coordinate))
allocate(this%unique_descriptor_tolerance(this%n_coordinate))
do i_coordinate = 1, this%n_coordinate
call initialise(params)
call param_register(params, 'energy_scale', "//MANDATORY//", delta, &
help_string="Set the typical scale of the function you are fitting (or the specific term if you use multiple descriptors). It is equivalent to the standard deviation of the Gaussian process in the probabilistic view, and typically this would be &
set to the standard deviation (i.e. root mean square) of the function &
that is approximated with the Gaussian process. ", altkey="delta")
call param_register(params, 'f0', '0.0', f0, &
help_string="Set the mean of the Gaussian process. Defaults to 0.")
call param_register(params, 'n_sparse', "0", n_sparseX, &
help_string="Number of sparse points to use in the sparsification of the Gaussian process")
call param_register(params, 'config_type_n_sparse', '', config_type_n_sparseX_string, &
help_string="Number of sparse points in each config type. Format: {type1:50:type2:100}")
call param_register(params, 'sparse_method', 'RANDOM', sparse_method_str, &
help_string="Sparsification method. RANDOM(default), PIVOT, CLUSTER, UNIFORM, KMEANS, COVARIANCE, NONE, FUZZY, FILE, &
INDEX_FILE, CUR_COVARIANCE, CUR_POINTS")
call param_register(params, 'lengthscale_factor', '1.0', theta_fac_string, has_value_target = this%has_theta_fac(i_coordinate), &
help_string="Set the width of Gaussians for the Gaussian and PP kernel by multiplying the range of each descriptor by lengthscale_factor. &
Can be a single number or different for each dimension. For multiple theta_fac separate each value by whitespaces.", altkey="theta_fac")
call param_register(params, 'lengthscale_uniform', '0.0', theta_uniform, has_value_target = this%has_theta_uniform(i_coordinate), &
help_string="Set the width of Gaussians for the Gaussian and PP kernel, same in each dimension.", altkey="theta_uniform")
call param_register(params, 'lengthscale_file', '', theta_file, has_value_target = this%has_theta_file(i_coordinate), &
help_string="Set the width of Gaussians for the Gaussian kernel from a file. &
There should be as many real numbers as the number of dimensions, in a single line", altkey="theta_file")
call param_register(params, 'sparse_file', '', sparse_file, has_value_target = has_sparse_file, &
help_string="Sparse points from a file. If sparse_method=FILE, descriptor values (real) listed in a text file, one &
& >>element<< per line. If sparse_method=INDEX_FILE, 1-based index of sparse points, one per line.")
call param_register(params, 'mark_sparse_atoms', 'F', mark_sparse_atoms, &
help_string="Reprints the original xyz file after sparsification process. &
sparse propery added, true for atoms associated with a sparse point.")
call param_register(params, 'add_species', 'T', add_species, &
help_string="Create species-specific descriptor, using the descriptor string as a template.")
call param_register(params, 'covariance_type', "//MANDATORY//", covariance_type_str, &
help_string="Type of covariance function to use. Available: Gaussian, DOT_PRODUCT, BOND_REAL_SPACE, PP (piecewise polynomial)")
!call param_register(params, 'theta', '1.0', main_gap_fit%theta(i_coordinate), &
!help_string="Width of Gaussians for use with bond real space covariance.")
call param_register(params, 'soap_exponent', '1.0', zeta, has_value_target = this%has_zeta(i_coordinate), &
help_string="Exponent of soap type dot product covariance kernel", altkey="zeta")
call param_register(params, 'print_sparse_index', '', print_sparse_index, &
help_string="If given, after determinining the sparse points, their 1-based indices are appended to this file")
call param_register(params, 'unique_hash_tolerance', '1.0e-10', unique_hash_tolerance, &
help_string="Hash tolerance when filtering out duplicate data points")
call param_register(params, 'unique_descriptor_tolerance', '1.0e-10', unique_descriptor_tolerance, &
help_string="Descriptor tolerance when filtering out duplicate data points")
if (.not. param_read_line(params, this%gap_str(i_coordinate), ignore_unknown=.true., task='main program gap_str('//i_coordinate//')')) then
call system_abort("main program failed to parse gap string ("//i_coordinate//")='"//trim(this%gap_str(i_coordinate))//"'")
endif
call finalise(params)
this%delta(i_coordinate) = delta
this%f0(i_coordinate) = f0
this%n_sparseX(i_coordinate) = n_sparseX
this%config_type_n_sparseX_string(i_coordinate) = config_type_n_sparseX_string
this%theta_fac_string(i_coordinate) = theta_fac_string
this%theta_uniform(i_coordinate) = theta_uniform
this%theta_file(i_coordinate) = theta_file
this%sparse_file(i_coordinate) = sparse_file
this%mark_sparse_atoms(i_coordinate) = mark_sparse_atoms
this%add_species(i_coordinate) = add_species
this%zeta(i_coordinate) = zeta
this%print_sparse_index(i_coordinate) = print_sparse_index
this%unique_hash_tolerance(i_coordinate) = unique_hash_tolerance
this%unique_descriptor_tolerance(i_coordinate) = unique_descriptor_tolerance
select case(lower_case(trim(sparse_method_str)))
case('random')
this%sparse_method(i_coordinate) = GP_SPARSE_RANDOM
case('pivot')
this%sparse_method(i_coordinate) = GP_SPARSE_PIVOT
case('cluster')
this%sparse_method(i_coordinate) = GP_SPARSE_CLUSTER
case('uniform')
this%sparse_method(i_coordinate) = GP_SPARSE_UNIFORM
case('kmeans')
this%sparse_method(i_coordinate) = GP_SPARSE_KMEANS
case('covariance')
this%sparse_method(i_coordinate) = GP_SPARSE_COVARIANCE
case('uniq')
call system_abort("sparse method UNIQ is no longer in use. Use NONE instead." )
case('fuzzy')
this%sparse_method(i_coordinate) = GP_SPARSE_FUZZY
case('file')
this%sparse_method(i_coordinate) = GP_SPARSE_FILE
case('index_file')
this%sparse_method(i_coordinate) = GP_SPARSE_INDEX_FILE
case('cur_covariance')
this%sparse_method(i_coordinate) = GP_SPARSE_CUR_COVARIANCE
case('cur_points')
this%sparse_method(i_coordinate) = GP_SPARSE_CUR_POINTS
case('none')
this%sparse_method(i_coordinate) = GP_SPARSE_NONE
case default
call system_abort("unknown sparse method "//trim(sparse_method_str))
endselect
if( has_sparse_file ) then
if( this%sparse_method(i_coordinate) /= GP_SPARSE_FILE .and. &
this%sparse_method(i_coordinate) /= GP_SPARSE_INDEX_FILE ) then
call system_abort('"sparse_file" specified in command line, but sparse method not "file" or "index_file"')
endif
endif
select case(lower_case(trim(covariance_type_str)))
case('none')
call system_abort("covariance type cannot be"//trim(covariance_type_str))
this%covariance_type(i_coordinate) = COVARIANCE_NONE
case('gaussian')
this%covariance_type(i_coordinate) = COVARIANCE_ARD_SE
case('ard_se') ! backwards compatibility
this%covariance_type(i_coordinate) = COVARIANCE_ARD_SE
case('dot_product')
this%covariance_type(i_coordinate) = COVARIANCE_DOT_PRODUCT
case('bond_real_space')
this%covariance_type(i_coordinate) = COVARIANCE_BOND_REAL_SPACE
case('pp')
this%covariance_type(i_coordinate) = COVARIANCE_PP
case default
call system_abort("unknown covariance type"//trim(covariance_type_str)//". Available: Gaussian, DOT_PRODUCT, BOND_REAL_SPACE, PP (piecewise polynomial)")
endselect
enddo
call print('Descriptors have been parsed')
endsubroutine gap_fit_parse_gap_str
subroutine read_fit_xyz(this)
type(gap_fit), intent(inout) :: this
type(cinoutput) :: xyzfile
integer :: n_con
logical :: file_exists
if( allocated(this%at) ) then
do n_con = 1, this%n_frame
call finalise(this%at(n_con))
enddo
deallocate(this%at)
this%n_frame = 0
endif
inquire(file=this%at_file, exist=file_exists)
if( .not. file_exists ) then
call system_abort("read_fit_xyz: at_file "//this%at_file//" could not be found")
endif
call initialise(xyzfile,this%at_file,mpi=this%mpi_obj)
this%n_frame = xyzfile%n_frame
allocate(this%at(this%n_frame))
do n_con = 1, this%n_frame
call read(xyzfile,this%at(n_con),frame=n_con-1)
call set_cutoff(this%at(n_con), this%max_cutoff)
call calc_connect(this%at(n_con))
enddo
call finalise(xyzfile)
if(this%n_frame <= 0) then
call system_abort("read_fit_xyz: "//this%n_frame//" frames read from "//this%at_file//".")
endif
endsubroutine read_fit_xyz
subroutine read_descriptors(this)
type(gap_fit), intent(inout) :: this
integer :: i
this%max_cutoff = 0.0_dp
if(allocated(this%my_descriptor)) then
do i = 1, size(this%my_descriptor)
call finalise(this%my_descriptor(i))
enddo
deallocate(this%my_descriptor)
endif
allocate(this%my_descriptor(this%n_coordinate))
do i = 1, this%n_coordinate
call initialise(this%my_descriptor(i),this%gap_str(i))
if( this%max_cutoff < cutoff(this%my_descriptor(i)) ) this%max_cutoff = cutoff(this%my_descriptor(i))
enddo
endsubroutine read_descriptors
subroutine fit_n_from_xyz(this)
type(gap_fit), intent(inout) :: this
logical :: do_collect_tasks, do_filter_tasks
type(Atoms) :: at
integer :: n_con
logical :: has_ener, has_force, has_virial, has_stress_3_3, has_stress_voigt, has_hessian, has_local_property, has_force_mask, exclude_atom
real(dp) :: ener, virial(3,3), stress_3_3(3,3)
real(dp) :: stress_voigt(6)
real(dp), pointer, dimension(:,:) :: f, hessian_eigenvector_j
real(dp), pointer, dimension(:) :: local_property
logical, pointer, dimension(:) :: force_mask
integer :: i, j, k
integer :: n_descriptors, n_cross, n_hessian
integer :: n_current, n_last
do_collect_tasks = (this%task_manager%active .and. .not. this%task_manager%distributed)
do_filter_tasks = (this%task_manager%active .and. this%task_manager%distributed)
if (allocated(this%n_cross)) deallocate(this%n_cross)
if (allocated(this%n_descriptors)) deallocate(this%n_descriptors)
allocate(this%n_cross(this%n_coordinate))
allocate(this%n_descriptors(this%n_coordinate))
this%n_cross = 0
this%n_descriptors = 0
this%n_ener = 0
this%n_force = 0
this%n_virial = 0
this%n_hessian = 0
this%n_local_property = 0
n_last = 0
do n_con = 1, this%n_frame
if (do_filter_tasks) then
if (this%task_manager%tasks(n_con)%worker_id /= this%task_manager%my_worker_id) cycle
end if
has_ener = get_value(this%at(n_con)%params,this%energy_parameter_name,ener)
has_force = assign_pointer(this%at(n_con),this%force_parameter_name, f)
has_virial = get_value(this%at(n_con)%params,this%virial_parameter_name,virial)
has_stress_voigt = get_value(this%at(n_con)%params,this%stress_parameter_name,stress_voigt)
has_stress_3_3 = get_value(this%at(n_con)%params,this%stress_parameter_name,stress_3_3)
has_hessian = get_value(this%at(n_con)%params,"n_"//this%hessian_parameter_name,n_hessian)
has_local_property = assign_pointer(this%at(n_con),this%local_property_parameter_name, local_property)
has_force_mask = assign_pointer(this%at(n_con),trim(this%force_mask_parameter_name),force_mask)
if( has_ener ) then
this%n_ener = this%n_ener + 1
endif
if( has_force ) then
do i = 1, this%at(n_con)%N
exclude_atom = .false.
if(has_force_mask) exclude_atom = force_mask(i)
if( .not. exclude_atom ) this%n_force = this%n_force + 3
enddo
endif
if( has_stress_voigt .or. has_stress_3_3 ) then
if( has_stress_voigt .and. has_stress_3_3 ) then
call system_abort("fit_n_from_xyz: conflict in stress between 6-vector and 9-vector (really 3x3 matrix)")
endif
! if has_stress is true, virial is available whether or not virial
! field has been detected
has_virial = .true.
endif
if( has_virial ) then
this%n_virial = this%n_virial + 6
endif
if( has_hessian ) then
this%n_hessian = this%n_hessian + n_hessian
at = this%at(n_con)
endif
if( has_local_property ) then
this%n_local_property = this%n_local_property + this%at(n_con)%N
endif
if( has_local_property .and. ( has_ener .or. has_force .or. has_virial .or. has_hessian ) ) then
call system_abort("fit_n_from_xyz: local_property and (energy or force or virial or hessian) present in configuration, currently not allowed.")
endif
do i = 1, this%n_coordinate
call descriptor_sizes(this%my_descriptor(i),this%at(n_con),n_descriptors,n_cross)
if( has_force ) then
this%n_cross(i) = this%n_cross(i) + n_cross*3
endif
if( has_virial ) then
this%n_cross(i) = this%n_cross(i) + n_cross*6
endif
this%n_descriptors(i) = this%n_descriptors(i) + n_descriptors
if( has_hessian ) then
do j = 1, n_hessian
if( .not. assign_pointer(this%at(n_con),trim(this%hessian_parameter_name)//j, hessian_eigenvector_j) ) &
call system_abort("fit_n_from_xyz: could not find the "//j//"th of "//n_hessian//" hessian eigenvector")
hessian_eigenvector_j = hessian_eigenvector_j / sqrt( sum(hessian_eigenvector_j**2) )
do k = -1, 1, 2
at%pos = this%at(n_con)%pos + k * this%hessian_delta * hessian_eigenvector_j
call set_cutoff(at,this%max_cutoff)
call calc_connect(at)
call descriptor_sizes(this%my_descriptor(i),at,n_descriptors,n_cross)
this%n_descriptors(i) = this%n_descriptors(i) + n_descriptors
this%n_cross(i) = this%n_cross(i) + n_descriptors
enddo
enddo
endif
enddo
if (do_collect_tasks) then
n_current = this%n_ener + this%n_local_property + this%n_force + this%n_virial + this%n_hessian
call task_manager_add_task(this%task_manager, n_current - n_last)
n_last = n_current
end if
call finalise(at)
enddo
if (.not. do_filter_tasks) then
call print_title("Report on number of descriptors found")
do i = 1, this%n_coordinate
call print("---------------------------------------------------------------------")
call print("Descriptor "//i//": "//this%gap_str(i))
call print("Number of descriptors: "//this%n_descriptors(i))
call print("Number of partial derivatives of descriptors: "//this%n_cross(i))
enddo
call print_title("")
end if
end subroutine fit_n_from_xyz
subroutine fit_data_from_xyz(this,error)
type(gap_fit), intent(inout) :: this
integer, optional, intent(out) :: error
type(inoutput) :: theta_inout
type(descriptor_data) :: my_descriptor_data
type(Atoms) :: at
integer :: d, n_con
logical :: has_ener, has_force, has_virial, has_stress_voigt, has_stress_3_3, has_hessian, has_local_property, &
has_config_type, has_energy_sigma, has_force_sigma, has_virial_sigma, has_virial_component_sigma, has_hessian_sigma, &
has_force_atom_sigma, has_force_component_sigma, has_local_property_sigma, has_force_mask, exclude_atom
real(dp) :: ener, ener_core, my_cutoff, energy_sigma, force_sigma, virial_sigma, hessian_sigma, local_property_sigma, &
grad_covariance_cutoff, use_force_sigma
real(dp), dimension(3) :: pos
real(dp), dimension(3,3) :: virial, virial_core, stress_3_3, virial_component_sigma
real(dp), dimension(6) :: stress_voigt
real(dp), dimension(:), allocatable :: theta, theta_fac, hessian, hessian_core, grad_data
real(dp), dimension(:), pointer :: force_atom_sigma
real(dp), dimension(:,:), pointer :: f, hessian_eigenvector_i, f_hessian, force_component_sigma
real(dp), dimension(:), pointer :: local_property
logical, dimension(:), pointer :: force_mask
real(dp), dimension(:,:), allocatable :: f_core
integer, dimension(:,:), allocatable :: force_loc, permutations
integer :: ie, i, j, n, k, l, i_coordinate, n_hessian, n_energy_sigma, n_force_sigma, n_force_atom_sigma, &
n_force_component_sigma, n_hessian_sigma, n_virial_sigma, n_local_property_sigma, n_descriptors, n_virial_component_sigma
integer, dimension(:), allocatable :: xloc, hessian_loc, local_property_loc
integer, dimension(3,3) :: virial_loc
integer :: i_config_type, n_config_type, n_theta_fac
character(len=STRING_LENGTH) :: config_type
character(len=THETA_LENGTH) :: theta_string
character(len=STRING_LENGTH), dimension(:), allocatable :: theta_string_array
INIT_ERROR(error)
if (this%task_manager%active) then
if (.not. this%task_manager%distributed) then
call system_abort("fit_data_from_xyz: Tasks are not distributed.")
end if
do n_con = 1, this%n_frame
if (this%task_manager%tasks(n_con)%worker_id /= this%task_manager%my_worker_id) then
call finalise(this%at(n_con))
end if
end do
end if
this%my_gp%do_subY_subY = merge(gap_fit_is_root(this), .true., this%task_manager%active)
my_cutoff = 0.0_dp
call gp_setParameters(this%my_gp,this%n_coordinate,this%n_ener+this%n_local_property,this%n_force+this%n_virial+this%n_hessian,this%sparse_jitter)
do i_coordinate = 1, this%n_coordinate
d = descriptor_dimensions(this%my_descriptor(i_coordinate))
call gp_setParameters(this%my_gp,i_coordinate, d, this%n_descriptors(i_coordinate), this%n_cross(i_coordinate), this%delta(i_coordinate), this%f0(i_coordinate), &
covariance_type=this%covariance_type(i_coordinate) )
call gp_addDescriptor(this%my_gp,i_coordinate,trim(this%gap_str(i_coordinate)))
allocate(permutations(d,descriptor_n_permutations(this%my_descriptor(i_coordinate))))
call descriptor_permutations(this%my_descriptor(i_coordinate),permutations)
call gp_setPermutations(this%my_gp,i_coordinate,permutations)
deallocate(permutations)
my_cutoff = max(my_cutoff,cutoff(this%my_descriptor(i_coordinate)))
enddo
call print_title("Report on number of target properties found in training XYZ:")
call print("Number of target energies (property name: "//trim(this%energy_parameter_name)//") found: "//sum(this%task_manager%MPI_obj, this%n_ener))
call print("Number of target local_properties (property name: "//trim(this%local_property_parameter_name)//") found: "//sum(this%task_manager%MPI_obj, this%n_local_property))
call print("Number of target forces (property name: "//trim(this%force_parameter_name)//") found: "//sum(this%task_manager%MPI_obj, this%n_force))
call print("Number of target virials (property name: "//trim(this%virial_parameter_name)//") found: "//sum(this%task_manager%MPI_obj, this%n_virial))
call print("Number of target Hessian eigenvalues (property name: "//trim(this%hessian_parameter_name)//") found: "//sum(this%task_manager%MPI_obj, this%n_hessian))
call print_title("End of report")
if( this%do_core ) call Initialise(this%core_pot, args_str=this%core_ip_args, param_str=string(this%quip_string))
n_energy_sigma = 0
n_force_sigma = 0
n_force_atom_sigma = 0
n_force_component_sigma = 0
n_virial_component_sigma=0
n_hessian_sigma = 0
n_virial_sigma = 0
n_local_property_sigma = 0
do n_con = 1, this%n_frame
if (.not. is_initialised(this%at(n_con))) cycle
has_ener = get_value(this%at(n_con)%params,this%energy_parameter_name,ener)
has_force = assign_pointer(this%at(n_con),this%force_parameter_name, f)
has_virial = get_value(this%at(n_con)%params,this%virial_parameter_name,virial)
has_virial_component_sigma = get_value(this%at(n_con)%params,'virial_component_'//trim(this%sigma_parameter_name),virial_component_sigma)
has_stress_voigt = get_value(this%at(n_con)%params,this%stress_parameter_name,stress_voigt)
has_stress_3_3 = get_value(this%at(n_con)%params,this%stress_parameter_name,stress_3_3)
has_hessian = get_value(this%at(n_con)%params,"n_"//this%hessian_parameter_name,n_hessian)
has_config_type = get_value(this%at(n_con)%params,this%config_type_parameter_name,config_type)
has_local_property = assign_pointer(this%at(n_con),this%local_property_parameter_name,local_property)
has_energy_sigma = get_value(this%at(n_con)%params,'energy_'//trim(this%sigma_parameter_name),energy_sigma)
has_force_sigma = get_value(this%at(n_con)%params,'force_'//trim(this%sigma_parameter_name),force_sigma)
has_virial_sigma = get_value(this%at(n_con)%params,'virial_'//trim(this%sigma_parameter_name),virial_sigma)
has_hessian_sigma = get_value(this%at(n_con)%params,'hessian_'//trim(this%sigma_parameter_name),hessian_sigma)
has_force_atom_sigma = assign_pointer(this%at(n_con),'force_atom_'//trim(this%sigma_parameter_name),force_atom_sigma)
has_force_component_sigma = assign_pointer(this%at(n_con),'force_component_'//trim(this%sigma_parameter_name),force_component_sigma)
has_local_property_sigma = get_value(this%at(n_con)%params,'local_property_'//trim(this%sigma_parameter_name),local_property_sigma)
has_force_mask = assign_pointer(this%at(n_con),trim(this%force_mask_parameter_name),force_mask)
if ((.not. has_virial) .and. (has_stress_3_3 .or. has_stress_voigt)) then
if (has_stress_voigt) then
virial(1,1) = stress_voigt(1)