Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Gauss star gonio #985

Open
wants to merge 13 commits into
base: master
Choose a base branch
from
Open

Gauss star gonio #985

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
b53c320
gauss star and goniometer support
dermen Apr 21, 2024
230f535
type gauss star
dermen Apr 21, 2024
93b023f
device wrapper
dermen Apr 21, 2024
4965e1b
forward models
dermen Apr 21, 2024
87db641
mos samples phil
dermen Apr 21, 2024
8ae7d71
ma on detector only
dermen Apr 21, 2024
bf5be16
deletes moved test
dermen Apr 21, 2024
6b608a6
script to test gaussian star fix
dermen Apr 22, 2024
eea05fa
plot fix for test
dermen Apr 22, 2024
c0cba77
remove unused import
dermen Apr 22, 2024
b928fbf
remove option
dermen Apr 22, 2024
31a1840
Merge branch 'master' into gauss_star_gonio
dermen May 4, 2024
25ddd67
update run_tests and fudge in forward model
dermen May 4, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions simtbx/diffBragg/device.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,7 @@ def __init__(self, gpu_id, force=False):
"""
self.gpu_id = gpu_id
self.force_kokkos_init = force
self.is_using_kokkos_gpu = USE_KOKKOS_GPU

def __enter__(self):
if USE_KOKKOS_GPU or self.force_kokkos_init:
Expand Down
3 changes: 3 additions & 0 deletions simtbx/diffBragg/phil.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1319,6 +1319,9 @@
.help = diffbragg offers CUDA support via the DIFFBRAGG_USE_CUDA=1 environment variable specification
.help = or openmp support using the OMP_NUM_THREADS flag
.help = The exascale only uses CUDA (will raise error if CUDA is not confugured)
mosaic_samples_override = None
.type = int
.help = Specify the number of mosaic spread samples
}
"""

Expand Down
2 changes: 1 addition & 1 deletion simtbx/diffBragg/src/diffBragg.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1967,7 +1967,7 @@ void diffBragg::add_diffBragg_spots(const af::shared<size_t>& panels_fasts_slows
db_cryst.Nd = Nd;
db_cryst.Ne = Ne;
db_cryst.Nf = Nf;

db_cryst.xtal_shape = xtal_shape;

db_beam.number_of_sources = sources;
db_beam.source_X = source_X;
Expand Down
11 changes: 3 additions & 8 deletions simtbx/diffBragg/src/diffBraggCUDA.cu
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,13 +39,6 @@ void diffBragg_sum_over_steps_cuda(
diffBragg_cudaPointers& cp,
timer_variables& TIMERS){


if (db_cryst.phisteps > 1){
printf("PHI (goniometer rotations) not supported in GPU code: phi0=%f phisteps=%d, phistep=%f\n", db_cryst.phi0, db_cryst.phisteps, db_cryst.phistep);
printf("SPINDLE VEC: %f %f %f\n",db_cryst.spindle_vec[0], db_cryst.spindle_vec[1], db_cryst.spindle_vec[2]);
exit(-1);
}

int numblocks;
int blocksize;
char* diffBragg_blocks = getenv("DIFFBRAGG_NUM_BLOCKS");
Expand Down Expand Up @@ -513,7 +506,9 @@ void diffBragg_sum_over_steps_cuda(
cp.data_freq, cp.data_trusted,
cp.FhklLinear_ASUid,
cp.Fhkl_channels,
cp.Fhkl_scale, cp.Fhkl_scale_deriv
cp.Fhkl_scale, cp.Fhkl_scale_deriv,
db_cryst.xtal_shape==GAUSS_STAR,
db_cryst.xtal_shape==SQUARE
);

error_msg(cudaGetLastError(), "after kernel call");
Expand Down
15 changes: 7 additions & 8 deletions simtbx/diffBragg/src/diffBraggKOKKOS.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -57,12 +57,6 @@ void diffBraggKOKKOS::diffBragg_sum_over_steps_kokkos(
cuda_flags& db_cu_flags,
// diffBragg_kokkosPointers& kp,
timer_variables& TIMERS) {
if (db_cryst.phi0 != 0 || db_cryst.phisteps > 1) {
printf(
"PHI (goniometer position) not supported in GPU code: phi0=%f phisteps=%d, phistep=%f\n",
db_cryst.phi0, db_cryst.phisteps, db_cryst.phistep);
exit(-1);
}

Kokkos::Tools::pushRegion("diffBragg_sum_over_steps_kokkos");
Kokkos::Tools::pushRegion("local detector, beam and crystal");
Expand Down Expand Up @@ -529,7 +523,10 @@ void diffBraggKOKKOS::diffBragg_sum_over_steps_kokkos(
m_data_freq, m_data_trusted,
m_FhklLinear_ASUid,
m_Fhkl_channels,
m_Fhkl_scale, m_Fhkl_scale_deriv);
m_Fhkl_scale, m_Fhkl_scale_deriv,
db_cryst.xtal_shape==GAUSS_STAR,
db_cryst.xtal_shape==SQUARE
);
} else {
kokkos_sum_over_steps(
Npix_to_model, m_panels_fasts_slows, m_floatimage, m_wavelenimage, m_d_Umat_images,
Expand Down Expand Up @@ -568,7 +565,9 @@ void diffBraggKOKKOS::diffBragg_sum_over_steps_kokkos(
m_data_freq, m_data_trusted,
m_FhklLinear_ASUid,
m_Fhkl_channels,
m_Fhkl_scale, m_Fhkl_scale_deriv
m_Fhkl_scale, m_Fhkl_scale_deriv,
db_cryst.xtal_shape==GAUSS_STAR,
db_cryst.xtal_shape==SQUARE
);
}

Expand Down
1 change: 1 addition & 0 deletions simtbx/diffBragg/src/diffBraggKOKKOS.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@
#include "simtbx/diffBragg/src/util.h"
#include "simtbx/diffBragg/src/util_kokkos.h"
#include "simtbx/diffBragg/src/diffBragg_refine_flag.h"
#include "simtbx/nanoBragg/nanotypes.h"

using vector_vec3_t = view_1d_t<KOKKOS_VEC3>;
using vector_mat3_t = view_1d_t<KOKKOS_MAT3>;
Expand Down
95 changes: 67 additions & 28 deletions simtbx/diffBragg/src/diffBragg_cpu_kernel.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -401,6 +401,13 @@ void diffBragg_sum_over_steps(
double fp_fdp_manager_dI[2] = {0,0};
double dI_latt_diffuse[6] = {0,0,0,0,0,0};

// compute unit cell volume
Eigen::Vector3d ap_vec(1e10*db_cryst.eig_B(0,0), 1e10*db_cryst.eig_B(1,0), 1e10*db_cryst.eig_B(2,0));
Eigen::Vector3d bp_vec(1e10*db_cryst.eig_B(0,1), 1e10*db_cryst.eig_B(1,1), 1e10*db_cryst.eig_B(2,1));
Eigen::Vector3d cp_vec(1e10*db_cryst.eig_B(0,2), 1e10*db_cryst.eig_B(1,2), 1e10*db_cryst.eig_B(2,2));
double cell_vol = ap_vec.dot(bp_vec.cross(cp_vec)); // cubic Angstrom
double xtal_size = pow(db_cryst.Na*db_cryst.Nb*db_cryst.Nc*cell_vol, (double)1/double(3)); // Angstrom

for (int i_step=0; i_step < db_steps.Nsteps; i_step++){

int subS = db_steps.subS_pos[i_step];
Expand Down Expand Up @@ -474,25 +481,21 @@ void diffBragg_sum_over_steps(
}
}

double phi = db_cryst.phi0 + db_cryst.phistep*phi_tic;
Eigen::Matrix3d Bmat_realspace = db_cryst.eig_B;
if( phi != 0.0 )
{
double cosphi = cos(phi);
double sinphi = sin(phi);
Eigen::Vector3d ap_vec(db_cryst.eig_B(0,0), db_cryst.eig_B(1,0), db_cryst.eig_B(2,0));
Eigen::Vector3d bp_vec(db_cryst.eig_B(0,1), db_cryst.eig_B(1,1), db_cryst.eig_B(2,1));
Eigen::Vector3d cp_vec(db_cryst.eig_B(0,2), db_cryst.eig_B(1,2), db_cryst.eig_B(2,2));

ap_vec = ap_vec*cosphi + db_cryst.spindle_vec.cross(ap_vec)*sinphi + db_cryst.spindle_vec*(db_cryst.spindle_vec.dot(ap_vec))*(1-cosphi);
bp_vec = bp_vec*cosphi + db_cryst.spindle_vec.cross(bp_vec)*sinphi + db_cryst.spindle_vec*(db_cryst.spindle_vec.dot(bp_vec))*(1-cosphi);
cp_vec = cp_vec*cosphi + db_cryst.spindle_vec.cross(cp_vec)*sinphi + db_cryst.spindle_vec*(db_cryst.spindle_vec.dot(cp_vec))*(1-cosphi);

Bmat_realspace << ap_vec[0], bp_vec[0], cp_vec[0],
ap_vec[1], bp_vec[1], cp_vec[1],
ap_vec[2], bp_vec[2], cp_vec[2];
}
Bmat_realspace *= 1e10;
Eigen::Matrix3d Bmat_realspace = 1e10*db_cryst.eig_B;
//if( phi != 0.0 )
//{
// double cosphi = cos(phi);
// double sinphi = sin(phi);

// ap_vec = ap_vec*cosphi + db_cryst.spindle_vec.cross(ap_vec)*sinphi + db_cryst.spindle_vec*(db_cryst.spindle_vec.dot(ap_vec))*(1-cosphi);
// bp_vec = bp_vec*cosphi + db_cryst.spindle_vec.cross(bp_vec)*sinphi + db_cryst.spindle_vec*(db_cryst.spindle_vec.dot(bp_vec))*(1-cosphi);
// cp_vec = cp_vec*cosphi + db_cryst.spindle_vec.cross(cp_vec)*sinphi + db_cryst.spindle_vec*(db_cryst.spindle_vec.dot(cp_vec))*(1-cosphi);

// Bmat_realspace << ap_vec[0], bp_vec[0], cp_vec[0],
// ap_vec[1], bp_vec[1], cp_vec[1],
// ap_vec[2], bp_vec[2], cp_vec[2];
//}
//Bmat_realspace *= 1e10;
if (db_flags.use_diffuse && db_flags.gamma_miller_units){
anisoG_local = anisoG_local * Bmat_realspace;
for (int i_gam=0; i_gam<3; i_gam++){
Expand All @@ -501,9 +504,28 @@ void diffBragg_sum_over_steps(
}

Eigen::Matrix3d U = db_cryst.eig_U;
Eigen::Matrix3d UBO = (db_cryst.UMATS_RXYZ[mos_tic] * U*Bmat_realspace*(db_cryst.eig_O.transpose())).transpose();

Eigen::Matrix3d Ainv = U*(Bmat_realspace.transpose().inverse())* (db_cryst.eig_O.inverse());
Eigen::Matrix3d UBOt=U*Bmat_realspace*(db_cryst.eig_O.transpose());
double phi = db_cryst.phi0 + db_cryst.phistep*phi_tic;
if (phi != 0){
if (i_pix==0){
printf("phistep=%f, phi0=%f, phi=%f\n", db_cryst.phistep, db_cryst.phi0, phi);
}

double c = cos(phi);
double omc = 1-c;
double s = sin(phi);
Eigen::Matrix3d Rphi;
double gx = db_cryst.spindle_vec[0];
double gy = db_cryst.spindle_vec[1];
double gz = db_cryst.spindle_vec[2];
Rphi << c + gx*gx*omc, gx*gy*omc-gz*s, gx*gz*omc+gy*s,
gy*gx*omc + gz*s, c + gy*gy*omc, gy*gz*omc - gx*s,
gz*gx*omc - gy*s, gz*gy*omc + gx*s, c + gz*gz*omc;
UBOt = Rphi*UBOt;
}
Eigen::Matrix3d UBO = (db_cryst.UMATS_RXYZ[mos_tic]*UBOt).transpose();
Eigen::Matrix3d Ainv = UBO.inverse();
Eigen::Vector3d q_vec(scattering[0], scattering[1], scattering[2]);
q_vec *= 1e-10;
Eigen::Vector3d H_vec = UBO*q_vec;
Expand All @@ -528,12 +550,30 @@ void diffBragg_sum_over_steps(
Eigen::Vector3d V = NABC*delta_H;
double hrad_sqr = V.dot(V);
double NABC_determ = NABC.determinant();
double F_latt;
if (db_flags.no_Nabc_scale)
F_latt = exp(-( hrad_sqr / 0.63 * db_cryst.fudge ));
else
//F_latt = db_cryst.Na*db_cryst.Nb*db_cryst.Nc*exp(-( hrad_sqr / 0.63 * db_cryst.fudge ));
F_latt = NABC_determ*exp(-( hrad_sqr / 0.63 * db_cryst.fudge ));
double F_latt = 0;
if (db_cryst.xtal_shape==SQUARE){
F_latt = 1;
if (db_cryst.Na > 1)
F_latt *= sincg(M_PI * h, db_cryst.Na);
if (db_cryst.Nb > 1)
F_latt *= sincg(M_PI * k, db_cryst.Nb);
if (db_cryst.Nc > 1)
F_latt *= sincg(M_PI * l, db_cryst.Nc);
}
else { // gaussian relp model
double exparg;
if (db_cryst.xtal_shape==GAUSS_STAR){
Eigen::Vector3d delta_Q = UBO.inverse()*delta_H;
double rad_star_sqr = delta_Q.dot(delta_Q)*xtal_size*xtal_size;
exparg = rad_star_sqr *1.9* db_cryst.fudge ;
}
else
exparg = hrad_sqr / 0.63 * db_cryst.fudge;

F_latt = exp(-exparg);
if (! db_flags.no_Nabc_scale)
F_latt *= NABC_determ;
}

/* convert amplitudes into intensity (photons per steradian) */
if (!db_flags.oversample_omega && !db_flags.Fhkl_gradient_mode)
Expand Down Expand Up @@ -749,7 +789,6 @@ void diffBragg_sum_over_steps(
// printf("hkl= %f %f %f hkl1= %d %d %d Fcell=%f\n", h,k,l,h0,k0,l0, F_cell);

double two_C = 2*C;
Eigen::Matrix3d UBOt = U*Bmat_realspace*(db_cryst.eig_O.transpose());
if (db_flags.refine_Umat[0]){
Eigen::Matrix3d RyRzUBOt = db_cryst.RotMats[1]*db_cryst.RotMats[2]*UBOt;
Eigen::Vector3d delta_H_prime = (db_cryst.UMATS[mos_tic]*db_cryst.dRotMats[0]*RyRzUBOt).transpose()*q_vec;
Expand Down
Loading