Protein Docking: discretization

discretize_molecules.jl

volume_dim = 128
receptor_volume, ligand_volume = zeros(volume_dim, volume_dim, volume_dim), zeros(volume_dim, volume_dim, volume_dim)
vols = Dict('l' => ligand_volume, 'r' => receptor_volume)

asa_min = 1.0

# receptor -- iteration 1: [core: sqrt(1.5)*r, surface: sqrt(0.8)*r], iteration 2: [core: 0.0 (skip), surface: r + 3.4]
rr_iter1 = [receptor_atoms_asa[i] < asa_min ? sqrt(1.5)*receptor_radius_list[i] : sqrt(0.8)*receptor_radius_list[i] for i in 1:length(receptor_atoms_asa)]
rr_iter2 = [receptor_atoms_asa[i] < asa_min ? 0.0 : receptor_radius_list[i]+3.4 for i in 1:length(receptor_atoms_asa)]

# ligand -- iteration 1: [core: sqrt(1.5)*r, surface: 0.0 (skip)], iteration 2: [core: 0.0 (skip), surface: 1.0 * r]
lr_iter1 = [ligand_atoms_asa[i] < asa_min ? sqrt(1.5)*ligand_radius_list[i] : 0.0 for i in 1:length(ligand_atoms_asa)]
lr_iter2 = [ligand_atoms_asa[i] < asa_min ? 0.0 : ligand_radius_list[i] for i in 1:length(ligand_atoms_asa)]

for (iter, designation, radius_list) in [(1, 'r', rr_iter1), (2, 'r', rr_iter2), (1, 'l', lr_iter1), (2, 'l', lr_iter2)]
    atoms = if designation == 'r' receptor_coords else ligand_coords end
    set_val = if iter == 1 2.0 else 1.0 end

    for atom_idx in 1:length(atoms)
        atom, atom_r = atoms[atom_idx], radius_list[atom_idx]

        vol_idxs = [1 + floor( ( atom[i] - vol_start ) / voxel_size ) for i in 1:3]
        n_neighbor_voxels = ceil(atom_r / voxel_size)
        start_end_idx = [(max(1, vol_idxs[i] - n_neighbor_voxels), min(vol_idxs[i] + n_neighbor_voxels, volume_dim)) for i in 1:3]

        for x in start_end_idx[1][1]:start_end_idx[1][2]
            for y in start_end_idx[2][1]:start_end_idx[2][2]
                for z in start_end_idx[3][1]:start_end_idx[3][2]
                    ( (iter == 2) && (designation == 'r') && (vols[designation][Int(x),Int(y),Int(z)] == 2.0) ) && continue

                    voxel_coords = [x, y, z]
                    voxel_center = [vol_start + voxel_size * (voxel_coords[i] - 0.5) for i in 1:3]
                    center_dist = norm(voxel_center - atoms[atom_idx])

                    (center_dist <= atom_r) && (vols[designation][Int(x),Int(y),Int(z)] = set_val)
                end
            end
        end
    end
end