Dear Geant4 users, I face a problem in simulation of Nanoparticles in my code, would any one can give me a clue how I can tackle this problem? the position of the Nanoparticles in the Geant4 code is as shown in snapshot below.
Are you sure “NP” (the number of nanoparticles) i supposed to be a double, rather than an integer?
You don’t show any code where that variable is being used, which makes it hard for anyone other than you to know how it’s supposed to be used.
Please don’t use screenshots for code. Copy and paste the code as text, using a “block” delimiter of three backquotes at the beginning and end:
```
void myFunction(G4int someArg) {
G4cout << someArg << " Hello, world!" << G4endl;
}
```
which will be displayed as:
void myFunction(G4int someArg) {
G4cout << someArg << " Hello, world!" << G4endl;
}
Thank you very much for response and Sorry, for the snapshot post. Yes you are true for NPs G4int or G4intLong is the true representation, how ever, the number of of the NPs are beyond it. When I run it gives error, but I changed to double it works and generate the nano particles. My challenge is when i increase the number beyond 10*7 it takes weeks and even doesn’t show any result. So is there any mechanism that makes my simulation fast? here is the code. To distribute the NPs I used std::vector and G4UniformRand
// Cylindrical Water Phantom
G4double pRmin = 0.cm;
G4double pRmax = 3.25cm;
G4double pHeight = 7.5cm;
G4double pSPhi = 0.deg;
G4double pDPhi = 360.deg;
G4Tubs solidTube = new G4Tubs(“solidTube”, pRmin, pRmax, pHeight, pSPhi, pDPhi);
G4LogicalVolume logicTube = new G4LogicalVolume(solidTube, water,“logicTube”);
G4VPhysicalVolume physTube = new G4PVPlacement(0, G4ThreeVector(0, 0, 0),
logicTube, “physTube”, logicPhantom, false, 0, checkOverlaps);
G4VisAttributes* solidTubeVisAtt = new G4VisAttributes(brown);
solidTubeVisAtt → SetVisibility(true);
solidTubeVisAtt → SetForceWireframe(true);
logicTube → SetVisAttributes(solidTubeVisAtt);
// generating positions inside the cylindrical phantom
G4double NP = 1.8e10; // Number of Sherical gold Nanoparticles in (18mg/g)
G4double r = 5e-6*cm; //radius of the Sherical Nanoparticles
G4double DMIN = 2*r; // Minimum distance between two Sherical Nanoparticles
G4ThreeVector fmin, fmax;
solidTube->BoundingLimits(fmin,fmax); // get bounding box of the cylinder
std::vector position;
while(position.size() < NP)
{
if you place more than G4INT_MAX (~2e9) volumes, the instanceID of type G4int for the physical volume will overflow. don’t know if that necessarily causes problems, but apparently it very likely does so 
also consider the probability to generate a free spot for the next NP to place. if the volume is already critically filled, it will take significantly longer and longer to randomly generate a non-overlapping location.
Such a huge number of objects as 1.8*10^10 nanoparticles cannot be generated in one volume. Just filling a std::vector with nanoparticle positions would require half a terabyte of operating memory and hundreds of years of CPU time:
One position (G4ThreeVector) takes 24 bytes, 1.8*10^10 positions will take 432 Gb
Current algorithm to generate NP positions performs more than NP^2/2 comparisons. Assuming that NP/2 = 0.9*10^10 comparisons can be done in 1 sec (actually more), it means that all comparisons will take 1.8*10^10 secs = 570 years
So, you should use much smaller volumes. For example, reduce the size of solidTube:
pRmax = 3.25*cm / 25; // 0.13*cm
pHeight = 7.5*cm / 25; // 0.3*cm
NP = 1.152e6;
Alternatively, you can try to construct the original volume as a set of sub-volumes (multiple replication of a very small volume).
Surprising, Thank you very much evc for your clarification. I will try both of your ways. In the 1st case when you reduce the pRmax and pHeight of the cylinder by 25, you put the NP = 1.152e6; Would you please explain how you determine the NP value = 1.152e6?. I mean from where comes from this 1.152e6 number?
Both, the radius and height of the cylinder are reduced in 25 times, it means that the volume of the cylinder is reduced in 25*25*25 = 15625 times. To keep the density of nanoparticles in the volume unchanged, the number of nanoparticles should also be reduced in 15625 times:
1.8*10^10 / (25*25*25) = 1.152*10^6
G4Box* tumor = new G4Box(“tumor”, 0.5cm, 0.5cm, 0.5cm);
G4LogicalVolume logicTumor = new G4LogicalVolume(tumor, softTissue, “logicTumor”);
G4VPhysicalVolume* physTumor = new G4PVPlacement(0, G4ThreeVector(0, 1*cm, 0), logicTumor, “physTumor”, logicTube, false, 0, checkOverlaps);
Here, I need to divide the tumor volume of 1cc in 125000 parts, that means in each part we will have 1.44*10^4 NPs. i.e each side of the cube is divided by 50.
The minimum and maximum value of the cube are[-0.5, 0.5) = G4ThreeVector fmin, fmax;
To distribute the NPs along the x-axis without dividing the cube we can use
G4double x = fmin.x() + (fmax.x() -fmin.x())*G4UniformRand();
but after dividing the into 125000 parts how we can distribute the NPs?
There are three logical volumes: big box, small box, and nanoparticle.
You fill the small box with randomly distributed 1.44*10^4 nanoparticles using an algorithm very similar to the one that was used for the cylinder.
Then, in three nested for loops, you fill the big box with 50x50x50 small boxes using corresponding translation vectors (G4ThreeVector).
Dear evc really appreciate for all you thoughtful suggestions. Based on the quoted above recommendation I wrote the three nested loops as below but I couldnt succeeded. Would you please suggest me what mistake I made?
//=========================================================================
//// Voxel
G4Box* voxel = new G4Box(“voxel”, 0.1mm, 0.1mm, 0.1mm);
G4LogicalVolume logicVoxel = new G4LogicalVolume(voxel, softTissue, “logicVoxel”);
G4int N =50;
for(G4int j = 0; j < N; j++)
{
for(G4int m = 0; m < N-2; m++)
{
for(G4int n = 0; n < N; n++)
G4VPhysicalVolume* physVoxel = new G4PVPlacement(0, G4ThreeVector(-0.5*cm+(j+0.5*cm)/N), 1*cm+(m+2)/N),-0.5*cm +(n+0.5*cm)/N)), 0), logicVoxel, "physVoxel", logicTumor, false, 0, checkOverlaps);
}
}
*The minimum and maximum value of the cube are[-0.5, 0.5)
Try the code below in visual mode. If the image looks Ok, then you can try to replace the parameters in the code with the ones you need.
// Main parameters
G4double sizeTumor = 10*mm;
G4double R = 0.05*mm; // radius of particle
G4int Nv = 7; // number of voxels = Nv * Nv * Nv
G4int Np = 100; // number of particles in voxel
G4double sizeVoxel = sizeTumor / Nv;
// Materials
G4Material* tumor_mat = nist->FindOrBuildMaterial("G4_WATER");
G4Material* voxel_mat = tumor_mat;
G4Material* nano_mat = nist->FindOrBuildMaterial("G4_Au");
// Make logical volumes
auto solidTumor = new G4Box("Tumor", 0.5*sizeTumor, 0.5*sizeTumor, 0.5*sizeTumor);
auto logicTumor = new G4LogicalVolume(solidTumor, tumor_mat, "Tumor_LV");
auto solidVoxel = new G4Box("Voxel", 0.5*sizeVoxel, 0.5*sizeVoxel, 0.5*sizeVoxel);
auto logicVoxel = new G4LogicalVolume(solidVoxel, voxel_mat, "Voxel_LV");
auto solidNano = new G4Orb("Nano", R);
auto logicNano = new G4LogicalVolume(solidNano, nano_mat, "Nano_LV");
// Fill Voxel with nanoparticles
std::vector<G4ThreeVector> pos;
G4double DMIN = 2*R;
while (pos.size() < Np) {
// generate position
G4double x = (0.5*sizeVoxel - R)*(2.*G4UniformRand() - 1.);
G4double y = (0.5*sizeVoxel - R)*(2.*G4UniformRand() - 1.);
G4double z = (0.5*sizeVoxel - R)*(2.*G4UniformRand() - 1.);
G4ThreeVector p(x, y, z);
// check position
G4int size = pos.size();
G4int k = 0;
for ( ; k < size; ++k) { if ((pos[k] - p).mag() < DMIN) break; }
// add position to the list
if (k == size) pos.push_back(p);
}
// Place nanoparticles at generated positions
for (G4int i = 0; i < Np; ++i) {
new G4PVPlacement(nullptr,pos[i],logicNano, "Nano_PV", logicVoxel, false,0,false);
}
// Fill Tumor with Voxels
G4double dd = 0.5*sizeVoxel;
for (G4int ix = 0; ix < Nv; ++ix) {
G4double x = -(0.5*sizeTumor - dd) + 2*ix*dd;
for (G4int iy = 0; iy < Nv; ++iy) {
G4double y = -(0.5*sizeTumor - dd) + 2*iy*dd;
for (G4int iz = 0; iz < Nv; ++iz) {
G4double z = -(0.5*sizeTumor - dd) + 2*iz*dd;
new G4PVPlacement(nullptr,G4ThreeVector(x,y,z),logicVoxel,"Voxel_PV",logicTumor, false,0,false);
}
}
}
This topic was automatically closed 7 days after the last reply. New replies are no longer allowed.