function autoRELAX = rhoSim(pdbNumCodeReduced,pdbNumCode,diagDecay,conditions,user, work)

fprintf('\nCaculation autorelaxation values')

if user.rhoSupValues == 1
    fprintf('\nand using experimental autorelaxation values if available')
end

deleteScreenMsg                                 = [];


length_pdbNumCodeReduced                        = length(pdbNumCodeReduced.Model.Atom);
autoRELAX                                       = zeros(length_pdbNumCodeReduced,4);

w                                               = user.B0field *1e6;
tauM                                            = user.tauC;
J0                                              = J(tauM, w, 0);
J1                                              = J(tauM, w, 1);
J2                                              = J(tauM, w, 2);


for i = 1:length_pdbNumCodeReduced
    
    
    if (mod(i,2) == 0) || i == length_pdbNumCodeReduced || i == 1
        screenMsg = sprintf(['\n', num2str(length_pdbNumCodeReduced + 1 - i), ' remaining.\n']);
        fprintf([deleteScreenMsg, screenMsg])
        deleteScreenMsg = repmat(sprintf('\b'), 1, length(screenMsg));
    
        if i == length_pdbNumCodeReduced
            fprintf(deleteScreenMsg)
            fprintf('\nDONE\n')
        end
    end
    
    
    ATOMi_resSeq                                = pdbNumCodeReduced.Model.Atom(i).resSeq;
    ATOMi_NumCode                               = pdbNumCodeReduced.Model.Atom(i).NumCode;
  	ATOMi_ser_num                               = pdbNumCodeReduced.Model.Atom(i).AtomSerNo;
    
    supflag = 0;
    
    if user.rhoSupValues == 1
        
        % use the calculated autorelaxation rate if available
        for j=1:length(diagDecay)
           if supflag == 0  
                if diagDecay(j).numcode.ATOMi_pdbSerNo  == ATOMi_ser_num
                
                    autoRELAX(i,1:3)                = [ATOMi_resSeq ATOMi_NumCode ATOMi_ser_num];
                    autoRELAX(i,4)                  = diagDecay(j).rho;
                    supflag                         = 1;
                
                end
           end
            
        end
        
        if supflag == 0 
            
        % reasign single proton of methyl groups ser number of pseudo
        % methyl to find corresponding relaxation rate (same for methyl
        % pseudo as single atoms
            
        for a= 1:length(conditions.methyl)
            if supflag == 0
                if length(conditions.methyl(a).atoms) == 3
                 	if size(find(ismember(ATOMi_NumCode, conditions.methyl(a).atoms)), 2) == 1
        
                            pseudo_numcode = conditions.methyl(a).grpNum;
                            
                            for s = 1:length(pdbNumCode.Model.Atom)
                            
                                if supflag == 0
                                    
                                    if ATOMi_resSeq == pdbNumCode.Model.Atom(s).resSeq && pseudo_numcode == pdbNumCode.Model.Atom(s).NumCode
                                    
                                    pseudo_ser_num  = pdbNumCode.Model.Atom(s).AtomSerNo;
                                
                                        for j=1:length(diagDecay)
            
                                            if diagDecay(j).numcode.ATOMi_pdbSerNo  == pseudo_ser_num
                
                                                autoRELAX(i,1:3)                = [ATOMi_resSeq ATOMi_NumCode ATOMi_ser_num];
                                                autoRELAX(i,4)                  = diagDecay(j).rho;
                                                supflag                         = 1;
                
                                            end
            
                                        end
                                    end
                                end
                
                            end
                    end
                end
            end
            
        end

        end
        
        % use the average autorelaxation rate provided by the user
        
        if supflag == 0 
            
                autoRELAX(i,1:3)                = [ATOMi_resSeq ATOMi_NumCode ATOMi_ser_num];
                autoRELAX(i,4)                  = user.rhoAverage;
                supflag                         = 1;
                
        end
        
    end
    
    
    if user.rhoSupValues == 0 || supflag == 0
        
    ATOMi_coor                                  = [pdbNumCodeReduced.Model.Atom(i).X pdbNumCodeReduced.Model.Atom(i).Y pdbNumCodeReduced.Model.Atom(i).Z];
    
    %************************************************************
    %        Get all protons involved           
    %************************************************************
    
    count_ATOMx                                 = 1; 
    gamma_ATOMx                                 = [];
    
    for j = 1:length_pdbNumCodeReduced
        
        skipFlag                                = 0;    
        ser_num                                 = pdbNumCodeReduced.Model.Atom(j).AtomSerNo;
        
        if ATOMi_ser_num ~= ser_num
        
            resSeq                              = pdbNumCodeReduced.Model.Atom(j).resSeq;
            NumCode                             = pdbNumCodeReduced.Model.Atom(j).NumCode;
            resName                             = pdbNumCodeReduced.Model.Atom(j).resName;
            
            % skip atoms belonging to same methylgroups, for these the formalism fails!
                            
          	if ATOMi_resSeq == resSeq
                
              	bothAtoms = [ATOMi_NumCode NumCode];
               	for a= 1:length(conditions.methyl)
        
                 	if length(conditions.methyl(a).atoms) == 3
                     	if size(find(ismember(bothAtoms, conditions.methyl(a).atoms)), 2) == 2
        
                            skipFlag            = 1;
                            
                            break
                
                        end
                    end
                end
                                
            end
            

            if skipFlag == 0

                coor_H = [ pdbNumCodeReduced.Model.Atom(j).X pdbNumCodeReduced.Model.Atom(j).Y pdbNumCodeReduced.Model.Atom(j).Z ];
                                
                    if norm( coor_H - ATOMi_coor ) <= user.upperlimit       	% is H close enough to ATOMi
                      
                        coor                        = [coor_H ATOMi_coor];
                        r                           = norm(coor(1:3) - coor(4:6));
                        
                        gamma                       = R(r,J0,J1,J2);
                        
                    	if gamma == Inf
        
                        	gamma  = 0;
        
                        end
                                        
                    	gamma_ATOMx(count_ATOMx)    = scaleRho2labSchemeFun(gamma,NumCode,resName,conditions,user);
                       	count_ATOMx                 = count_ATOMx + 1;            
                                        
                    end             
            end 
        
        end
    end
    
         	autoRELAX(i,1:3)                        = [ATOMi_resSeq ATOMi_NumCode ATOMi_ser_num];
            autoRELAX(i,4)                          = sum(gamma_ATOMx(:));
            
            
    end

end

autoRELAXav                                         = sum(autoRELAX(:,4))/length(autoRELAX(:,4));

for i = 1 : size(autoRELAX,1)
    if autoRELAX(i,4) == 0
       autoRELAX(i,4)                               =  autoRELAXav;
    end 
end

file = fopen (work.rhosimOutFileName, 'w'); 


fprintf (file,'%6d%6d%6d%12.3f\n', autoRELAX');

fclose(file);

while (~exist(work.rhosimOutFileName, 'file'))
  sleep(0.5);
end
      
