I remember that Pearson correlation works for continuous variables and also if one is continuous and the other a dummy. For the latter, it does not provide a correlation but provides a portion. What is the best approach, when we have 2 dummy variables. Would you suggest to calculate correlation for them or what is a good way to look at the relationship between 2 dichotomous variables. Ex: gender and religious status.
You're looking for estimates of association for contingency tables. You've got a few options. Here's a great summary of two of the most popular measures, the phi coefficient and Cramer's V: http://www.people.vcu.edu/~pdattalo/702SuppRead/MeasAssoc/NominalAssoc.html
Both measure the strength of the association between categorical variables. The chi coefficient only works for 2x2 tables (like the example you give: "2 dichotomous variables. Ex: gender and religious status"), but Cramer's V works for larger tables as well. For a 2x2 table, phi and V will be equal.
Phi is a chi-square based measure of association. The chi-square coefficient depends on the strength of the relationship and sample size. Phi eliminates sample size by dividing chi-square by n, the sample size, and taking the square root.
Cramer's V is the most popular of the chi-square-based measures of nominal association because it gives good norming from 0 to 1 regardless of table size, when row marginals equal column marginals. V equals the square root of chi-square divided by sample size, n, times m, which is the smaller of (rows - 1) or (columns - 1): V = SQRT(X2/nm).
Another measure similar to Cramer's V is Cohen's w, explained on this page: http://stats.idre.ucla.edu/other/mult-pkg/faq/general/effect-size-power/faqhow-is-effect-size-used-in-power-analysis/
Effect size w is the square root of the standardized chi-square statistic.