Help regarding parallelizing code

for i=1:(length(water)/3)-1
            
            for j=i+1:length(water)/3
                
                donor_mol_index=[i j]; %between 2 water molecules either can be donor/acceptor
                for p=donor_mol_index
                    donor_h_index=[1 3];% water atom data 1st and last atoms are hydrogen atoms
                    for k=donor_h_index
                        hb_ww_count=hb_ww_count+1;
                        
                        %di-donor index, ai-acceptor index, donor is the molecule which gives the h-atom 
                        if p==i
                            di=i; ai=j; %if donor molecule index equals i
                        else
                            di=j;ai=i;%if donor molecule index equals j
                        end
                        %dist_ww(cnt_steps,hb_ww_count)=sqrt(sum((water_mol(di).atom_data(k,3:5)-water_mol(ai).atom_data(2,3:5)).^2));
                        v1=water_mol(di).atom_data(k,3:5)-water_mol(di).atom_data(2,3:5);%vector connecting donor hydrogen-donor oxygen
                        v2=water_mol(di).atom_data(k,3:5)-water_mol(ai).atom_data(2,3:5);%vector connecting donor hydrogen-acceptor oxygen
                        angle=acos(dot(v1, v2) / (norm(v1) * norm(v2)));
                        angle_ww(cnt_steps,hb_ww_count)=radtodeg(angle);
                        if sum((water_mol(di).atom_data(k,3:5)-water_mol(ai).atom_data(2,3:5)).^2)<=dist_hb^2 %distance between donor hydrogen and acceptor oxygen
                            v1=water_mol(di).atom_data(k,3:5)-water_mol(di).atom_data(2,3:5);%vector connecting donor hydrogen-donor oxygen
                            v2=water_mol(di).atom_data(k,3:5)-water_mol(ai).atom_data(2,3:5);%vector connecting donor hydrogen-acceptor oxygen
                            angle=acos(dot(v1, v2) / (norm(v1) * norm(v2)));
                            angle=radtodeg(angle);
                            if angle<=angle_hb
                                hb_ww(cnt_steps,hb_ww_count)=1;
                            end
                        else
                            hb_ww(cnt_steps,hb_ww_count)=0;
                        end
                    end
                end
            end
        end

I would like to parallelize the outer loop. Can I put all the code inside the first for loop in a function and supply the loop variable i and the structures to the function ? hb_ww_count is a reduction variable. I use this variable for indexing another array within my code.

Please provide some suggestions. I am sort of lost.

[Jan]

Start with simplifying your code: why do you calculate v1 and v2 and the angle twice? Change:

v1=water_mol(di).atom_data(k,3:5)-water_mol(di).atom_data(2,3:5);%vector connecting donor hydrogen-donor oxygen
v2=water_mol(di).atom_data(k,3:5)-water_mol(ai).atom_data(2,3:5);%vector connecting donor hydrogen-acceptor oxygen
angle=acos(dot(v1, v2) / (norm(v1) * norm(v2)));
angle_ww(cnt_steps,hb_ww_count)=radtodeg(angle);
if sum((water_mol(di).atom_data(k,3:5)-water_mol(ai).atom_data(2,3:5)).^2)<=dist_hb^2 %distance between donor hydrogen and acceptor oxygen
    v1=water_mol(di).atom_data(k,3:5)-water_mol(di).atom_data(2,3:5);%vector connecting donor hydrogen-donor oxygen
    v2=water_mol(di).atom_data(k,3:5)-water_mol(ai).atom_data(2,3:5);%vector connecting donor hydrogen-acceptor oxygen
    angle=acos(dot(v1, v2) / (norm(v1) * norm(v2)));

to:

% Before the loops:
dist_hb2 = dist_hb ^ 2;  % Once only
len = floor(length(water)/3);
n = 4 * (len - 1) * len / 2;
hb_ww = zeros(cnt_steps, n);   % Pre-allocate!!!
angle_ww = zeros(cnt_steps, n);
...
    v1 = water_mol(di).atom_data(k,3:5) - water_mol(di).atom_data(2,3:5);%vector connecting donor hydrogen-donor oxygen
    v2 = water_mol(di).atom_data(k,3:5) - water_mol(ai).atom_data(2,3:5);%vector connecting donor hydrogen-acceptor oxygen
    angle = acos(dot(v1, v2) / (norm(v1) * norm(v2))); % Not accurate, prefer atan2
    angle_ww(cnt_steps, hb_ww_count) = angle * (180 / pi);
    if sum(v2 .^ 2) <= dist_hb2 %distance between donor hydrogen and acceptor oxygen
       if angle<=angle_hb
         hb_ww(cnt_steps, hb_ww_count)=1;
       end
    end
end

With avoiding repeated calculation and a proper pre-allocation the code is at least efficient without parallelizing.

[Soumyadipta Sengupta]

preallocation of all the arrays and structures have been done but not shown here.

Answer 1

[Soumyadipta Sengupta]

Hello Jan,

I have preallocated those variables as well as the structures . I didn't show the complete code.

[Jan]

Please post comments inthe section for comments. Not showing the relevant part of the code, especially the pre-allocation, is not helpful for the readers.

[Soumyadipta Sengupta]

Even without the double angle calculation, each inner loop takes around 0.5 secs on a 2.6 ghz core. i would have to run approx 2000 such inner loops and 2000 such outer loops

Answer 2

[Soumyadipta Sengupta]

i have parallelized it.