-
Notifications
You must be signed in to change notification settings - Fork 0
/
XNDiff.cpp
12907 lines (11064 loc) · 455 KB
/
XNDiff.cpp
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
/***************************************************************************
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program; if not, write to the *
* Free Software Foundation, Inc., *
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
***************************************************************************/
/*** History: ***/
/*
The development of a program for the evaluation of very complicated functions
was started in Oktober 2000 by T. Unruh. This program was called "func_calc".
The main goal was to write a program for the calculation of diffraction patterns
of very complicated structures numerically.
After first qualitatively successful tests in 2000 first quantitatively convincing
results could be achieved in 2007. At this stage the program was renamed to "XDiff"
and made available by Sourceforge.net ( http://sourceforge.net/projects/xdiff/).
15.10.2000: func_calc Version 0.02
av_corr applied to S not to E !
28.10.2000: Version 0.03
karthesian coordinates of a vectors corrected
hkl averaging (mode 1) bugs corrected
dispersion medium correction added
25.11.2000: func_calc Version 0.04
some speed optimizations
pid included in outputs
27.11.2000: func_calc Version 0.05
some changes of output data
log file name changed to <mfname>.log
3.11.2000: func_calc Version 0.06
V = V * n1 * n2 * n3!
some minor bugs
4.11.2000: func_calc Version 0.07
minor bug fix in output file
21.12.2000: func_calc Version 0.08
minor speed optimizations
11.01.2005: func_calc Version 0.09
background bug fixes
21.02.2006: func_calc Added stabilizer scattering
17.03.2006: func_calc Added new crystal (plus stabilizer) overlap test
26.06.2007: XDiff version 0.01 extracted from func_calc
(cf. J. Appl. Cryst. submitted)
The code is very preliminary and only intended to simulate very special
systems so far
26.06.2009: Begin of further developement of XNDiff
- add neutron diffraction support
- optimization for speed
- add crystallographic stuff such as space groups and symmetry operations
- add a tool for the generation of structure parameter files with complete molecules
from e.g. CSD (cif) or PDB file format
20140406 minor bug fix, for New crystal: the printed w3 and wt3 vectors were wrong w(t)3[0][k][:] was printed (i.e. always n3=1) but not w(t)3[n3[k]-1][k][:] for the correct n3[k] thickness
only printed results in the log-files were wrong but not the computations
*/
/**********************************************/
/* short description of the stackcpp function */
/**********************************************/
/* (to be added later) */
/*
folder structure:
XNDiff/ : XNDiff.cpp mtrand.cpp mtrand.h compile_XNDiff.sh Makefile
XNDiff/in/ : syminfo.lib (optional) f0_InterTables.dat DeBe_NeutronNews.dat atomweights.dat cif_core.dic
further user-provided input files like cif-files (less common pdb-files), symmetry operation files, parameter files (less common, parameters can be provided upon call of XNDiff on cmd)
XNDiff/out/ : XNDiff generated files for SAXS and SANS patterns, log-files, further output files for testing (distributions, (modified)cif-files, ...)
*/
/*
mandatory and optional external libraries / installations to run XNDiff on GNU Linux and (with minor code modifcations) on other POSIX systems (FreeBSD, ...):
-MTRand PRNG (MANDATORY, might be replaced with a PRNG from GSL later),
http://www.bedaux.net/mtrand/
just place mtrand.cpp and mtrand.h at the same location of XNDiff.cpp before compilation
-GSL (MANDATORY, tested with v1.15)
http://www.gnu.org/software/gsl/
comes along with most GNU Linux distros and is installed on most HPC systems by default
-CCP4 (OPTIONAL, tested with v6.1.13, when using -sym <ccp4 ...> option)
http://www.ccp4.ac.uk/
Note: CCP4 can be binary installation, since CCP4 is called via opening a shell
-Armadillo (OPTIONAL, tested with v3.0.3, mostly for test purposes, when using -bt <bt_twin, bt_double> option)
*/
/* compilation using g++, assuming everything is installed / put in the right place:
-without using optional Armadillo library, Optimization: DEBUG / O3-Level, no further architecture-specific optimization
e.g. ./compile_XNDiff.sh g++ DEBUG 0 NONE
e.g. ./compile_XNDiff.sh g++ 3 0 NONE
-using optional Armadillo library, Optimization: DEBUG / O3-Level, no further architecture-specific optimization
e.g. ./compile_XNDiff.sh g++ DEBUG 1 NONE
e.g. ./compile_XNDiff.sh g++ 3 1 NONE
*/
/*
further mandatory / optional parameter-files with physical constants and symmetry operations etc to run XNDiff:
-download from link below and save as syminfo.lib to in/
http://www.ccp4.ac.uk/cvs/viewvc.cgi/libccp4/data/syminfo.lib?view=co
(OPTIONAL, symmetry operations, required when using -sym <symopfile> or -sym <ccp4 ...> option without providing a symopfile, see source code for further details)
(further reading http://www.ccp4.ac.uk/html/symlib.html)
-download from link below and save as f0_InterTables.dat to in/
http://ftp.esrf.eu/pub/scisoft/xop2.3/DabaxFiles/f0_InterTables.dat
(MANDATORY, Cromer-Mann coefficients for X-ray scattering)
-download from link below and save as DeBe_NeutronNews.dat to in/
http://ftp.esrf.eu/pub/scisoft/DabaxFiles/DeBe_NeutronNews.dat
(MANDATORY, neutron coherent scattering lengths etc)
-save as atomweights.dat the ASCII table from link below into in/ and make some modifications (has been already done and included in your download of XNDiff!!!)
http://physics.nist.gov/cgi-bin/Compositions/stand_alone.pl?ele=&ascii=ascii&isotype=some
(MANDATORY, atom/isotope masses)
-download from link below the current version and save as cif_core.dic to in/
ftp://ftp.iucr.org/pub/cif_core.dic
(MANDATORY, cif-dictionary required to read cif-files, tested with v2.4.3 from 2012-05-16)
*/
#include <iostream>
#include <algorithm>
#include <string.h>
#include <sstream>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <sys/time.h>
#include <vector>
#include <ctype.h>
#include <time.h>
/* directory exists http://bytes.com/topic/c/answers/584434-check-directory-exists-c */
#include <dirent.h>
/* mkdir http://www.imb-jena.de/~gmueller/kurse/c_c++/c_file.html */
#include <sys/stat.h>
/* Unix-Default-Lib-Header (Windows C++ stdlib.h ) */
#include <unistd.h>
/* for typeid */
#include <typeinfo>
/**************************************/
/* OpenMP support for parallelization */
/**************************************/
#include <omp.h>
/* print variable name http://ubuntuforums.org/showthread.php?t=506402 */
#define var_name(var) #var
/**************************************************************************/
/* COMPLEX NUMBER HEADER OR CLASSES http://en.wikipedia.org/wiki/Complex.h*/
/**************************************************************************/
#include <complex>
/***************************************************************************/
/* PRNG Mersenne-Twister MTRand for MC mode (much better than C++'s rand())*/
/* http://www.bedaux.net/mtrand/ */
/***************************************************************************/
#include "mtrand.h"
/*****************************************************************/
/* GNU Scientific Library (GSL), PRNG for discrete distributions */
/* http://www.gnu.org/software/gsl/ */
/*****************************************************************/
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
using namespace std;
/****************************************/
/* ARMADILLO LIBRARY FOR LINEAR ALGEBRA */
/* http://arma.sourceforge.net/ */
/****************************************/
#if (ARMADILLO)
#include <armadillo>
#endif
#if (ARMADILLO)
using namespace arma;
#endif
typedef complex<double> dcmplx;
class XNDiff
{
public:
/* program information */
const char* XNDIFF_VERSION ; /* see constructor */
const char* XNDIFF_DATE ; /* see constructor */
const char* XNDIFF_AUTHOR ; /* see constructor */
double DR ; /* see constructor */
double RD ; /* see constructor */
/* List of available functions for XPPSA Method */
/* struct function, with func[i]->name and func[i].nop */
/* List of available functions for basis transformations */
static const unsigned int NOF = 1 ;
static const unsigned int NOBT = 2 ;
bool bt_flag ;
char bt[1024] ; /* bt function name */
struct function
{
const char *name ; /* Function name */
unsigned int nop ; /* Number of parameters */
} func_XPPSA[NOF], func_bt[NOBT] ;
struct par_entries
{
const char *THICKNESS_INLAY ;
const char *THICKNESS_OUTLAY ;
const char *XRAY_RHO_INLAY ;
const char *XRAY_RHO_OUTLAY ;
const char *XRAY_RHO_DM ;
const char *NEUT_SLD_INLAY ;
const char *NEUT_SLD_OUTLAY ;
const char *NEUT_SLD_DM ;
const char *NEUT_ICSD_INLAY ;
const char *NEUT_ICSD_OUTLAY ;
const char *NEUT_ICSD_DM ;
const char *CONC_CRY ;
const char *CONC_DM ;
} par_keyword ;
/* process ID http://www.tutorials.de/forum/linux-unix/158373-datentypen-pid_t-dev_t-ssize_t.html */
pid_t pid ;
unsigned int time_flag ; /* no protocolling(0), only at beginning/FFF/terminating(1) or detailed(2) */
unsigned int mem_flag ; /* no protocolling(0), only at beginning(1) or detailed(2) */
char mem_filename[64] ;
/* log file and output files */
bool singlefilesoutput_flag ;
bool extension_flag ;
bool intermediatesteps_flag ;
unsigned int imsteps ;
int log_flag;
FILE *logfile ;
char mfname[1024] ; /* Master file name of output files */
char **evarg ; /* Export function call parameters */
int *ecarg ; /* Export function call parameters */
/* XNDiff PRNG for particles, ... */
#define R1 16807
#define R2 2147483647
#define R3 (1.0/R2)
#define R4 127773
#define R5 2836
#define R6 32
#define R7 (1+(R2-1)/R6)
#define R8 1.1e-7
#define R9 (1.0-R8)
/* Parameters are contained in the structure par of type struct stack_parameters */
struct stack_parameters
{ /* shared parameters */
double n11, n12 ;
double n21, n22 ;
double n31, n32 ;
double d11, d12 ;
double d21, d22 ;
double D1, D2 ;
int td1 ;
int td2 ;
int stackmode ;
double rs1, rs2 ;
unsigned int np ;
unsigned int nr ;
double h, k, l ;
double hs, ks, ls ;
unsigned int nsp ; /* number of layers in single particles e.g. 10 */
unsigned int nms ; /* maximum number of particles in stacks e.g. 5 */
int av_mode ;
double av_pol[3] ;
double av_azi[3] ;
double r31, r32 ;
double av_r3[2] ; /* range of r3-rotation in [ av_r3[0] , av_r3[1] ] */
unsigned int n_r3 ; /* number of grid points for r3-rotation between -90° and +90°, e.g. 180 */
unsigned int n_MC ; /* number of points for MC mode e.g. 180*24 = 4320 */
} *par ;
static const int MAX_AV_MODE = 3 ; /* see +av flag */
unsigned int nav_pc, nav_ac, nav_r3, nav_MC ;
double av_r3d ;
double ds, s1 ;
vector <double> av_pol_ang ;
vector <double> av_azi_ang ;
char par_file[1024] ; /* filename of the .par-file which provides the necessary information about the stabilizer layer and dispersion media */
char ccp4_config_file[1024] ; /* configuration file with path information for an installed CCP4-library */
char ccp4_bash_file[1024] ; /* filename of the bash-script used for running the CCP4 (pdbset) */
const char *def_ccp4_bash_file ;
char cif_dic[1024] ; /* cif-dictionary filename */
char cif_file[1024] ; /* cif-file filename */
char cif_pdb_file[1024] ;
const char *def_cif_pdb_file ;
bool jmol_cif_flag ; /* write cif-file for molecular viewers like Jmol,PyMol,Rasmol,... */
bool sym_pdb_mode ;
char pdb_file[1024] ; /* pdb-filename either by default (def_pdb_file) generated by CCP4 from cif_pdb_file or user-supplied */
const char *def_pdb_file ;
bool pcr_flag ; /* flag for writing a file with atomic coordinates in the pcr-file format for FullProf */
char pcr_file[1024] ;
struct cif_dic_entries /* structure of entries in the cif-dictionary */
{
bool isnumeric;
char* units;
char* category;
char** names;
int num_names;
} *cde;
unsigned int num_cif_dic_entries; /* number of entries in the cif-dictionary */
struct cif_file_entries /* structure of entries in the cif-file */
{
bool isnumeric ;
char* units ;
char* name ;
unsigned int data_length ;
char** text_data ;
double* numeric_data ;
double* numeric_stderr_data ; /* default 0.0 */
} *cfe ;
unsigned int num_cif_file_entries ; /* number of entries in the cif-file */
bool sym_ccp4_mode ;
bool sym_cif_mode ;
bool symop_userdef ; /* flag whether to use a user-supplied symmetry operations file */
char symop_file[1024] ; /* filename of the user-supplied symmetry operations file */
char** symopxyz; /* string array of symmetry operations */
unsigned int num_symop; /* amount of symmetry operations */
char** cenopxyz; /* string array of center operations */
unsigned int num_cenop; /* amount of center operations */
vector<double> cis ; /* c_i's fitted from diluted systems (\sum_i c_i will be normalized to 1 after read in) */
vector<double> nis ; /* n_i = (c_i/i)/[\sum_i(c_i/i)] in case cis are provided, else if nis provided normaliyed to 1 */
bool cis_isdef ;
bool nis_isdef ;
vector< vector<unsigned int> > bond_site12 ; /* connectivity table derived from cif-file */
vector< vector<string> > bond_symmetry12 ;
vector<double> bond_distance ;
/* Standard Atomic Weight for elements and Relative Atomic Mass for isotopes */
struct StAtWtandRelAtMa
{
char* nuclei_name;
double mass; /* Standard Atomic Weight for elements and Relative Atomic Mass for isotopes */
} *atom_mass;
int num_atom_mass_entries;
/* NEutron SCattering LEngth and CRoss SECtions */
struct nesclecrsec
{
char* nuclei_name;
double nat_abd; /* natural abundance [0-100 %] (for isotopes the half-life is given instead) */
dcmplx coh_sc_length; /* complex bound coherent scattering length [fm] */
dcmplx incoh_sc_length; /* complex bound incoherent scattering length [fm] */
double coh_sc_cross; /* bound coherent scattering cross section [barn] */
double incoh_sc_cross; /* bound incoherent scattering cross section [barn] */
double total_sc_cross; /* total bound scattering cross section [barn] */
double abs_cross; /* absorption cross section for 2200 m/s neutrons [barn] */
double mass; /* Standard Atomic Weight for elements and Relative Atomic Mass for isotopes */
} *n_atom_sc;
int num_n_atom_sc_entries;
double cell_par[6] ; /* Unit cell parameters, edge lengths and angles of the triclinic unit cell in [nm] */
struct atom_entries /* Atom parameters */
{
char label[10] ; /* Atom label */
char type[5] ; /* Atom type (cf. get_atomic_scattering_factor */
double coord[3] ; /* Atom fractional coordinates */
double mult ; /* Multiplicity */
double temp[7] ; /* 1 isotropic and 6 anisotropic temperature factors (optional) */
dcmplx coh_sc_length ; /* coherent scattering length */
double incoh_sc_cross ;/* incoherent scattering cross section */
int f0ind ; /* index to a table with the pre-computed scattering amplitudes as a function of s, f_0(s) */
} *atom ;
unsigned int noa ; /* Number of atoms in the unit cell */
double** f0table ; /* numf0 x par->np table containing f_0^j(s), j=0,...,numf0 */
int numf0 ;
bool par_num_def ;
double *rho_dm ; /* electron density of dispersion medium [electrons/nm^3] */
double *sld_dm ; /* neutron SLD of dispersion medium [1e-6 A^(-2)] */
double *icsd_dm ; /* neutron ICSD of the dispersion medium [1/cm] */
double thickness_isl ; /* Thickness of inner stabilizer layer [nm] */
double *rho_isl ; /* electron density of inner stabilizer layer [electrons/nm^3] */
double *sld_isl ; /* neutron SLD of inner stabilizer layer [1e-6 A^(-2)] */
double *icsd_isl ; /* neutron ICSD of inner stabilizer layer [1/cm] */
double thickness_osl ; /* Thickness of outer stabilizer layer [nm] */
double *rho_osl ; /* electron density of outer stabilizer layer [electrons/nm^3] */
double *sld_osl ; /* neutron SLD of outer stabilizer layer [1e-6 A^(-2)] */
double *icsd_osl ; /* neutron ICSD of outer stabilizer layer [1/cm] */
double rho_cry ; /* X-ray electron density for the crystal [electrons/nm^3] */
double icsd_cry ; /* neutron incoherent cross section density for the crystal [1/cm] */
double sld_cry ; /* neutron scattering length density for the crystal [10^(-6) 1/A^2] */
bool conc_userdef ; /* flag whether to use default concentrations for conc_dm and conc_cry ( false, 1.0 and 0.01) or those specified in the par-file (true) */
double conc_dm ; /* vol% concentration of dispersion medium in the nanosuspension */
double conc_cry ; /* vol% concentration of crystals in the nanosuspension */
struct rho_multi_struct
{
bool multiq ;
unsigned int n ;
char* XRAY_RHO_INLAY_file ;
unsigned int XRAY_RHO_INLAY_col ;
char* XRAY_RHO_OUTLAY_file ;
unsigned int XRAY_RHO_OUTLAY_col ;
char* XRAY_RHO_DM_file ;
unsigned int XRAY_RHO_DM_col ;
} rho_multi ;
struct sld_multi_struct
{
bool multiq ;
unsigned int n ;
char* NEUT_SLD_INLAY_file ;
unsigned int NEUT_SLD_INLAY_col ;
char* NEUT_SLD_OUTLAY_file ;
unsigned int NEUT_SLD_OUTLAY_col ;
char* NEUT_SLD_DM_file ;
unsigned int NEUT_SLD_DM_col ;
char* NEUT_ICSD_INLAY_file ;
unsigned int NEUT_ICSD_INLAY_col ;
char* NEUT_ICSD_OUTLAY_file ;
unsigned int NEUT_ICSD_OUTLAY_col ;
char* NEUT_ICSD_DM_file ;
unsigned int NEUT_ICSD_DM_col ;
} sld_multi ;
double r_e ; /* see constructor */
double sc_scl ; /* see constructor */
double sc_rho ; /* see constructor */
double sc_sld ; /* see constructor */
double sc_icsd ; /* see constructor */
double sc_irr ; /* see constructor */
const char *unit_S_X, *unit_S_n ; /* see constructor */
const char *unit_rho, *unit_sld, *unit_icsd ; /* see constructor */
const char *unit_F_X, *unit_F_n ; /* see constructor */
const char *unit_Yc_X[10], *unit_Yc_n[10], *unit_Yi_n[4] ; /* see constructor */
const char *unit_V_CRY, *unit_V_ISL, *unit_V_OSL, *unit_V_DM, *unit_V_IRR ; /* see constructor */
int Xn ; /* flag, that decides whether X-ray (Xn=-1), neutron (Xn=+1) or both patterns (average of both, Xn=0) should be computed */
dcmplx ***FFF_X ; /* structure amplitude for X-rays [nm] */
dcmplx ***FFF_n ; /* structure amplitude for neutrons [nm] */
bool orav_flag ;
bool openmp_flag ;
bool Fr_flag ;
int num_Fr_flag ;
char Fr_file[2][1024] ;
bool Fw_flag ;
bool FNoMem_flag ;
bool distr_flag ;
bool weight_flag ;
int init_n_rand ;
bool init_n_rand_userdef ;
long unsigned int init_d_rand ;
bool init_d_rand_userdef ;
unsigned int init_MTRand ;
bool init_MTRand_userdef ;
double ***eQvj ; /* scalar products for all three a_j's with all the unit vectors of Q depending on theta,phi/r3 */
double a1[3], a2[3], a3[3] ; /* Lattice parameters in a Cartesian coordinate system */
double b1[3], b2[3], b3[3] ; /* Reciprocal unit cell vectors in kartesian coordinate system */
double G0[3] ; /* Reciprocal lattice unit vector parallel to (s-s0) */
double kg1[3], kg2[3] ; /* Two vectors in (hkl)-plane */
double Gs[3] ; /* Reciprocal lattice unit vector of stacking direction */
double k1[3], k2[3] ; /* Two vectors in (hs ks ls)-plane */
double V_uc ; /* Volume of unit cell = det|a_1 a_2 a_3| [nm^3] */
char startingtime[64];
const char *def_ifo, *def_ofo ;
char ifo[1024], ofo[1024] ;
bool overwrite_flag ;
bool silent_flag ;
bool test_XnEquiv_flag ;
bool test_I_at_s_equal_zero_flag ;
public:
/* constructor */
XNDiff()
{
/* program informations */
this->func_XPPSA[0].name = "stackcpp" ;
this->func_XPPSA[0].nop = 27 ;
this->func_bt[0].name = "bt_twin" ;
this->func_bt[0].nop = 0 ;
this->func_bt[1].name = "bt_double" ;
this->func_bt[1].nop = 0 ;
this->XNDIFF_VERSION = "0.4.1" ;
this->XNDIFF_DATE = "8.2000-04.2014" ;
this->XNDIFF_AUTHOR = "Tobias Unruh and Martin Schmiele" ;
/* class constants */
this->DR = M_PI/180.0 ;
this->RD = (180.0/M_PI) ;
this->r_e = 2.8179402894 ; /* classical electron radius in [fm] */
this->sc_scl = 1.0e-6 ; /* scattering lengths / r_e in [fm] -> [nm] for F_X/n and sc_rho */
this->sc_rho = r_e * sc_scl ; /* [ 1/nm^3 * ( r_e = 2.82 fm ) ] -> [1/nm^2] */
this->sc_sld = 1.0e-4 ; /* [10^(-6) 1/A^2 = 10^10 1/cm^2] -> [1/nm^2] */
this->sc_icsd = 1.0e-7 ; /* incoherent cross section density (ICSD) in [1/cm], volume in [nm^3], aim is [nm^2] for B_n */
this->sc_irr = 1.0e-7 ; /* [nm^3] for V -> [nm^2 * cm] for Vtot_irr */
/* units */
this->unit_S_X = "1/cm" ; /* absolute units */
this->unit_S_n = "1/cm" ; /* absolute units */
this->unit_rho = "1/nm^3" ; /* user-friendly unit for rho, note that r_e * rho can be expressed as for neutrons in [10^10 1/cm^2] */
this->unit_sld = "10^(-6) 1/A^2" ; /* user-friendly unit for sld, equivalent to [10^10 1/cm^2] */
this->unit_icsd = "1/cm" ; /* absolute units */
this->unit_F_X = "nm" ; /* atom form factor F */
this->unit_F_n = "nm" ;
this->unit_V_CRY = "nm^3" ;
this->unit_V_ISL = "nm^3" ;
this->unit_V_OSL = "nm^3" ;
this->unit_V_DM = "nm^3" ;
this->unit_V_IRR = "nm^2 cm" ;
this->unit_Yc_X[0] = "1/cm" ; /* absolute units */
for ( int q=1; q<7; ++q) { this->unit_Yc_X[q] = "nm^6/cm" ; } /* absolute units */
for ( int q=7; q<10; ++q) { this->unit_Yc_X[q] = "nm^3/cm" ; } /* absolute units */
this->unit_Yc_n[0] = "1/cm" ; /* absolute units */
for ( int q=1; q<7; ++q) { this->unit_Yc_n[q] = "nm^4/cm" ; } /* absolute units */
for ( int q=7; q<10; ++q) { this->unit_Yc_n[q] = "nm^2/cm" ; } /* absolute units */
this->unit_Yi_n[0] = "1/cm" ; /* absolute units */
for ( int q=1; q<4; ++q) { this->unit_Yi_n[q] = "1" ; } /* absolute units */
this->par_keyword.THICKNESS_INLAY = "THICKNESS_INLAY" ;
this->par_keyword.THICKNESS_OUTLAY = "THICKNESS_OUTLAY" ;
this->par_keyword.XRAY_RHO_INLAY = "XRAY_RHO_INLAY" ;
this->par_keyword.XRAY_RHO_OUTLAY = "XRAY_RHO_OUTLAY" ;
this->par_keyword.XRAY_RHO_DM = "XRAY_RHO_DM" ;
this->par_keyword.NEUT_SLD_INLAY = "NEUT_SLD_INLAY" ;
this->par_keyword.NEUT_SLD_OUTLAY = "NEUT_SLD_OUTLAY" ;
this->par_keyword.NEUT_SLD_DM = "NEUT_SLD_DM" ;
this->par_keyword.NEUT_ICSD_INLAY = "NEUT_ICSD_INLAY" ;
this->par_keyword.NEUT_ICSD_OUTLAY = "NEUT_ICSD_OUTLAY" ;
this->par_keyword.NEUT_ICSD_DM = "NEUT_ICSD_DM" ;
this->par_keyword.CONC_CRY = "CONC_CRY" ;
this->par_keyword.CONC_DM = "CONC_DM" ;
this->num_symop = 0 ;
this->num_cenop = 0 ;
this->noa = 0 ;
}
/* check if a directory path exists */
bool directory_exists( const char* path)
{
if ( path == NULL) return false ;
DIR *dir ;
bool ret = false ;
dir = opendir(path) ;
if ( dir != NULL )
{
ret = true ;
(void) closedir(dir) ;
}
return ret ;
}
/* Test computer code if it returns the same results for neutrons as for X-rays when replacing neutron parameters with X-ray ones
replaced will be neutron slds by rhos and scaling factors
FFF_n will be defacto copied from FFF_X elsewhere in the code, see:
compute_structure_amplitude()
read_structure_amplitude()
stackcpp()
*/
void test_XnEquiv()
{
/* test works only for X-ray AND neutron scattering */
if ( Xn != 0 ) { XNDIFF_ERROR(82) ; }
/* test should better not be applied when using multiple rho/sld/icsd's '*/
if ( rho_multi.multiq || sld_multi.multiq ) { XNDIFF_ERROR(83) ; }
/* coherent neutron scattering -> X */
this->sc_sld = this->sc_rho ;
sld_isl[0] = rho_isl[0] ;
sld_osl[0] = rho_osl[0] ;
sld_dm[0] = rho_dm[0] ;
/* no incoherent scattering */
this->sc_icsd = 0.0 ;
icsd_isl[0] = 0.0 ;
icsd_osl[0] = 0.0 ;
icsd_dm[0] = 0.0 ;
/* units cannot be replaced because the of const char* */
}
/* Test computer code if it returns the correct intensity at I(s=0)
Use :
conc_cry = 1.0 (and optionally conc_dm = 0.0)
to normalize only by volume of crystals (without stabilizer layer)
Note that in XNDiff Vtot_irr ~ Vtot_cry / conc_cry, hence only conc_cry counts,
but not conc_dm, which is only relevant for Vtot_dm (via VN_dm) ~ Vtot_cry * conc_dm / conc_cry
and by this for B_n and Yi_n
Hence with conc_cry = 1.0 it holds Vtot_irr = Vtot_cry
Generally normalization is done in XNDiff via :
S_X/n -> S_X/n / MULT(coh) / SUM_WEIGHT_CORR / Vtot_irr
Yc_X/n -> Yc_X/n / MULT(coh) / SUM_WEIGHT_CORR / Vtot_irr
and
B_n -> B_n / MULT(inc) / Vtot_irr
Yi_n -> Yi_n / MULT(inc) / Vtot_irr
In exported files for S_X/n and Yc_X/n etc
V_CRY ~ Vtot_cry ~ volume of crystal
V_ISL ~ Vtot_isl ~ volume of crystal AND inner stabilizer layer
V_OSL ~ Vtot_osl ~ volume of crystal AND inner AND outer stabilizer layers
V_DM and V_IRR correspond to Vtot_dm and Vtot_irr as given from above via Vtot_cry
V_CRY corresponds physically to the volume of crystal and therefore for XNDiff or later the correct conc_cry should be applied !!!
( V_ISL - V_CRY )
and
( V_OSL - V_ISL )
correspond to volumes of stabilizer layer covering the parallelepipedal particles. Ideally they should be also related to the concentration of the
stabilizer layer (e.g. phospholipid + sodium glycocholate etc ). However they must not, cause the stabilizer layer must not be densely packed and dispersion
medium and crystal substance may contribute. Their contributions might be later estimated by the rho/sld and from the difference between the expected volume
concentration of the stabilizer (sample composition) and the one from the simulation/fit.
*/
void test_I_at_s_equal_zero()
{
/* NOT IMPLEMENTED YET */
}
#if (ARMADILLO)
/* basis transformation for twinned triglyceride crystals */
void bt_twin()
{
char sdummy[1024] ;
int idummy ;
vec a1_vec = zeros<vec>(3) ;
vec a2_vec = zeros<vec>(3) ;
vec a3_vec = zeros<vec>(3) ;
vec delta_vec = zeros<vec>(3) ;
vec ap1_vec = zeros<vec>(3) ;
vec ap2_vec = zeros<vec>(3) ;
vec ap3_vec = zeros<vec>(3) ;
mat A = zeros<mat>(3,3) ; /* A=(a1,a2,a3) */
mat Ap = zeros<mat>(3,3); /* Ap=(ap1,ap2,ap3) */
mat T = zeros<mat>(3,3); /* inv(A) */
vec x_vec = zeros<vec>(3) ; /* x_vec=(x,y,z) */
vec xi_vec = zeros<vec>(3) ; /* xi_vec=(xi,eta,zeta) */
/* setup a1_vec,...,a3_vec, delta_vec and ap1_vec,...,ap3_vec */
for ( unsigned int i = 0; i<3; ++i)
{
a1_vec(i) = a1[i] ;
a2_vec(i) = a2[i] ;
a3_vec(i) = a3[i] ;
}
delta_vec(0) = 0.0 ;
delta_vec(1) = 0.0 ;
delta_vec(2) = 2.0 * a3_vec(2) ;
ap1_vec = a1_vec ;
ap2_vec = a2_vec ;
ap3_vec = delta_vec ;
if ( log_flag ) { fprintf( logfile, "\t a1_vec = ( %15.10G, %15.10G, %15.10G)\n", a1_vec[0], a1_vec[1], a1_vec[2]) ; fflush (logfile) ; }
if ( log_flag ) { fprintf( logfile, "\t a2_vec = ( %15.10G, %15.10G, %15.10G)\n", a2_vec[0], a2_vec[1], a2_vec[2]) ; fflush (logfile) ; }
if ( log_flag ) { fprintf( logfile, "\t a3_vec = ( %15.10G, %15.10G, %15.10G)\n", a3_vec[0], a3_vec[1], a3_vec[2]) ; fflush (logfile) ; }
if ( log_flag ) { fprintf( logfile, "\n") ; }
if ( log_flag ) { fprintf( logfile, "\tap1_vec = ( %15.10G, %15.10G, %15.10G)\n", ap1_vec[0], ap1_vec[1], ap1_vec[2]) ; fflush (logfile) ; }
if ( log_flag ) { fprintf( logfile, "\tap2_vec = ( %15.10G, %15.10G, %15.10G)\n", ap2_vec[0], ap2_vec[1], ap2_vec[2]) ; fflush (logfile) ; }
if ( log_flag ) { fprintf( logfile, "\tap3_vec = ( %15.10G, %15.10G, %15.10G)\n", ap3_vec[0], ap3_vec[1], ap3_vec[2]) ; fflush (logfile) ; }
/* setup matrix A=(a1,a2,a3) and Ap=(ap1,ap2,ap3) */
A.col(0) = a1_vec ;
A.col(1) = a2_vec ;
A.col(2) = a3_vec ;
Ap.col(0) = ap1_vec ;
Ap.col(1) = ap2_vec ;
Ap.col(2) = ap3_vec ;
/* set up transformation matrix from old coordinates/basis to new ones, T = inv(Ap) */
T = Ap.i(true) ;
/**********************/
/* update cif-entries */
/**********************/
/* update symmetry */
strcpy( sdummy, "_symmetry_space_group_name_H-M") ;
idummy = find_index_in_cfe(sdummy) ;
if ( idummy >= 0 )
{
strcpy( sdummy, "Unknown") ;
cfe[idummy].text_data[0] = (char *) realloc ( cfe[idummy].text_data[0], ( strlen(sdummy) + 1 ) * sizeof(char)) ;
strcpy( cfe[idummy].text_data[0], sdummy);
}
strcpy( sdummy, "_symmetry_space_group_name_Hall") ;
idummy = find_index_in_cfe(sdummy) ;
if ( idummy >= 0 )
{
strcpy( sdummy, "Unknown") ;
cfe[idummy].text_data[0] = (char *) realloc ( cfe[idummy].text_data[0], ( strlen(sdummy) + 1 ) * sizeof(char)) ;
strcpy( cfe[idummy].text_data[0], sdummy);
}
strcpy( sdummy, "_symmetry_equiv_pos_as_xyz") ;
idummy = find_index_in_cfe(sdummy) ;
if ( idummy >= 0 )
{
for ( unsigned int i=0; i<cfe[idummy].data_length; ++i)
{
free( cfe[idummy].text_data[i] ) ;
}
free( cfe[idummy].text_data ) ;
cfe[idummy].data_length = 1 ;
cfe[idummy].text_data = (char **) calloc( 1, sizeof(char *)) ;
strcpy( sdummy, "+x,+y,+z") ;
cfe[idummy].text_data[0] = (char *) calloc( strlen(sdummy)+1, sizeof(char)) ;
strcpy( cfe[idummy].text_data[0], sdummy);
}
/* cell, update c and alpha, beta */
strcpy( sdummy, "_cell_length_c") ;
idummy = find_index_in_cfe(sdummy) ;
if ( idummy >= 0 ) { cfe[idummy].numeric_data[0] = 10.0 * norm( ap3_vec, 2) ; }
strcpy( sdummy, "_cell_angle_alpha") ;
idummy = find_index_in_cfe(sdummy) ;
if ( idummy >= 0 ) { cfe[idummy].numeric_data[0] = 90.0 ; }
strcpy( sdummy, "_cell_angle_beta") ;
idummy = find_index_in_cfe(sdummy) ;
if ( idummy >= 0 ) { cfe[idummy].numeric_data[0] = 90.0 ; }
/* update atom with coordinates */
atom = (atom_entries *) realloc (atom, 2*noa * sizeof(atom_entries)) ;
for ( unsigned int i = 0; i<noa; ++i)
{
/* first unit cell remains the same */
x_vec(0) = atom[i].coord[0] ;
x_vec(1) = atom[i].coord[1] ;
x_vec(2) = atom[i].coord[2] ;
/* define position vector with old basis and coordinates x, y, z */
x_vec = A * x_vec ;
/* transform to new basis coordinates xi, eta, zeta */
xi_vec = T * x_vec ;
atom[i].coord[0] = xi_vec(0) ;
atom[i].coord[1] = xi_vec(1) ;
atom[i].coord[2] = xi_vec(2) ;
/* (001)-mirrored coordinates for second unit cell, use x_vec from above,
define position vector with old basis and coordinates x, y, z */
/* x_vec hasn't changed, use it again */
x_vec(2) *= -1 ;
x_vec += delta_vec ;
/* transform to new basis coordinates xi, eta, zeta */
xi_vec = T * x_vec ;
atom[i+noa].coord[0] = xi_vec(0) ;
atom[i+noa].coord[1] = xi_vec(1) ;
atom[i+noa].coord[2] = xi_vec(2) ;
}
/* update noa and atom labels and types */
noa *= 2 ;
derive_atom_labels_and_types_from_cif() ;
}
#endif
#if (ARMADILLO)
/* basis transformation, test scenario with a doubled unit cell along a3-direction for triglyceride crystals */
void bt_double()
{
char sdummy[1024] ;
int idummy ;
vec a1_vec = zeros<vec>(3) ;
vec a2_vec = zeros<vec>(3) ;
vec a3_vec = zeros<vec>(3) ;
vec ap1_vec = zeros<vec>(3) ;
vec ap2_vec = zeros<vec>(3) ;
vec ap3_vec = zeros<vec>(3) ;
mat A = zeros<mat>(3,3) ; /* A=(a1,a2,a3) */
mat Ap = zeros<mat>(3,3); /* Ap=(ap1,ap2,ap3) */
mat T = zeros<mat>(3,3); /* inv(A) */
vec x_vec = zeros<vec>(3) ; /* x_vec=(x,y,z) */
vec xi_vec = zeros<vec>(3) ; /* xi_vec=(xi,eta,zeta) */
/* setup a1_vec,...,a3_vec, delta_vec and ap1_vec,...,ap3_vec */
for ( unsigned int i = 0; i<3; ++i)
{
a1_vec(i) = a1[i] ;
a2_vec(i) = a2[i] ;
a3_vec(i) = a3[i] ;
}
ap1_vec = a1_vec ;
ap2_vec = a2_vec ;
ap3_vec = 2.0 * a3_vec ;
if ( log_flag ) { fprintf( logfile, "\t a1_vec = ( %15.10G, %15.10G, %15.10G)\n", a1_vec[0], a1_vec[1], a1_vec[2]) ; fflush (logfile) ; }
if ( log_flag ) { fprintf( logfile, "\t a2_vec = ( %15.10G, %15.10G, %15.10G)\n", a2_vec[0], a2_vec[1], a2_vec[2]) ; fflush (logfile) ; }
if ( log_flag ) { fprintf( logfile, "\t a3_vec = ( %15.10G, %15.10G, %15.10G)\n", a3_vec[0], a3_vec[1], a3_vec[2]) ; fflush (logfile) ; }
if ( log_flag ) { fprintf( logfile, "\n") ; }
if ( log_flag ) { fprintf( logfile, "\tap1_vec = ( %15.10G, %15.10G, %15.10G)\n", ap1_vec[0], ap1_vec[1], ap1_vec[2]) ; fflush (logfile) ; }
if ( log_flag ) { fprintf( logfile, "\tap2_vec = ( %15.10G, %15.10G, %15.10G)\n", ap2_vec[0], ap2_vec[1], ap2_vec[2]) ; fflush (logfile) ; }
if ( log_flag ) { fprintf( logfile, "\tap3_vec = ( %15.10G, %15.10G, %15.10G)\n", ap3_vec[0], ap3_vec[1], ap3_vec[2]) ; fflush (logfile) ; }
/* setup matrix A=(a1,a2,a3) and Ap=(ap1,ap2,ap3) */
A.col(0) = a1_vec ;
A.col(1) = a2_vec ;
A.col(2) = a3_vec ;
Ap.col(0) = ap1_vec ;
Ap.col(1) = ap2_vec ;
Ap.col(2) = ap3_vec ;
/* set up transformation matrix from old coordinates/basis to new ones, T = inv(Ap) */
T = Ap.i(true) ;
/**********************/
/* update cif-entries */
/**********************/
/* update symmetry */
strcpy( sdummy, "_symmetry_space_group_name_H-M") ;
idummy = find_index_in_cfe(sdummy) ;
if ( idummy >= 0 )
{
strcpy( sdummy, "Unknown") ;
cfe[idummy].text_data[0] = (char *) realloc ( cfe[idummy].text_data[0], ( strlen(sdummy) + 1 ) * sizeof(char)) ;
strcpy( cfe[idummy].text_data[0], sdummy);
}
strcpy( sdummy, "_symmetry_space_group_name_Hall") ;
idummy = find_index_in_cfe(sdummy) ;
if ( idummy >= 0 )
{
strcpy( sdummy, "Unknown") ;
cfe[idummy].text_data[0] = (char *) realloc ( cfe[idummy].text_data[0], ( strlen(sdummy) + 1 ) * sizeof(char)) ;
strcpy( cfe[idummy].text_data[0], sdummy);
}
strcpy( sdummy, "_symmetry_equiv_pos_as_xyz") ;
idummy = find_index_in_cfe(sdummy) ;
if ( idummy >= 0 )
{
for ( unsigned int i=0; i<cfe[idummy].data_length; ++i)
{
free( cfe[idummy].text_data[i] ) ;
}
free( cfe[idummy].text_data ) ;
cfe[idummy].data_length = 1 ;
cfe[idummy].text_data = (char **) calloc( 1, sizeof(char *)) ;
strcpy( sdummy, "+x,+y,+z") ;
cfe[idummy].text_data[0] = (char *) calloc( strlen(sdummy)+1, sizeof(char)) ;
strcpy( cfe[idummy].text_data[0], sdummy);
}
/* cell, update c */
strcpy( sdummy, "_cell_length_c") ;
idummy = find_index_in_cfe(sdummy) ;
if ( idummy >= 0 ) { cfe[idummy].numeric_data[0] = 10.0 * norm( ap3_vec, 2) ; }
/* update atom with coordinates */
atom = (atom_entries *) realloc (atom, 2*noa * sizeof(atom_entries)) ;
for ( unsigned int i = 0; i<noa; ++i)
{
/* first unit cell remains the same */
x_vec(0) = atom[i].coord[0] ;
x_vec(1) = atom[i].coord[1] ;
x_vec(2) = atom[i].coord[2] ;
/* define position vector with old basis and coordinates x, y, z */
x_vec = A * x_vec ;
/* transform to new basis coordinates xi, eta, zeta */
xi_vec = T * x_vec ;
atom[i].coord[0] = xi_vec(0) ;
atom[i].coord[1] = xi_vec(1) ;
atom[i].coord[2] = xi_vec(2) ;
/* second unit cell is the same as the first one, but shifted by a3, use x_vec from above */
/* define position vector with old basis and coordinates x, y, z */
/* x_vec' = A * x_vec hasn't changed, use it again */
/* x_vec = x_vec' + a3_vec */
x_vec += a3_vec ;
/* transform to new basis coordinates xi, eta, zeta */
xi_vec = T * x_vec ;
atom[i+noa].coord[0] = xi_vec(0) ;
atom[i+noa].coord[1] = xi_vec(1) ;
atom[i+noa].coord[2] = xi_vec(2) ;
}
/* update noa and atom labels and types */
noa *= 2 ;