C, 208 205 175 169 bytes

argv[1] : cation
argv[2] : anion
Charges of ions are taken on stdin.

#define z(b)for(i=*++v;*++i>95;);printf(b>1?*i?"(%s)_%d":"%s_%d":"%s",*v,b); main(c,v,d,e,g,h,i)char**v,*i;{scanf("%d%d",&c,&d);for(e=c,h=-d;g=h;e=g)h=e%g;z(-d/e)z(c/e)} 

Retina, 86 80 bytes

Thanks to Neil for saving 6 bytes.

^ ( \^.+ )_$&$* (1+)(\1|¶.+)+_(\1)+$ $#3$2_$#2 _1$ m)T`()``.+\)$|\(.[a-z]?\) ¶ 

Try it online!

Input is linefeed-separated (test suite uses comma-separation for convenience).

Explanation

^ ( 

We start by prepending a ( to each molecule. The ^ matches on line beginnings because the m) towards the end of the program sets multiline mode for all preceding stages.

\^.+ )_$&$* 

We replace the ^[-+]n part with )_, followed by n copies of 1 (i.e. we convert the charges to unary, dropping the signs).

(1+)(\1|¶.+)+_(\1)+$ $#3$2_$#2 

This stage does three things: it divides both charges by their GCD, converts them back to decimal and swaps them. The GCD can be found quite easily in regex, by matching the longest 1+ that lets us match both charges using only the backreference \1. To divide by this, we make use of Retina's "capture count" feature, which tells us how often a group has been used. So $#2 is the first charge divided by the GCD and $#3 is the second charge divided by the GCD (both in decimal).

_1$ 

We remove _1s from the ends of both parts.

m)T`()``.+\)$|\(.[a-z]?\) 

And we drop the parentheses from lines which end in a ) (i.e. those that had a _1 there), as well as lines that only contain a single atom.

¶ 

Finally, we concatenate the two molecules by dropping the linefeed.

Haskell, 101 97 bytes

(s#t)n m|let x#i|j<-div(lcm n m)i=l$x:[l(x:['(':x++")"|l x<':'])++'_':show j|j>1]=s#n++t#m l=last 

Try it online! Example usage: Fe^+3, O^-2 is taken as ("Fe"#"O")3 2.

RPL (HP48 S/SX), 294.5 bytes

Yes, ridiculously large submission, not sure how competitive it will be...

DIR M « S ROT S SWAP ABS 4 PICK ABS DUP2 WHILE DUP2 REPEAT MOD SWAP END DROP2 SWAP OVER / 4 ROLL J 3 ROLLD / ROT J ROT 0 > IF THEN SWAP END + » S « DUP "^" POS DUP2 1 + OVER SIZE SUB OBJ🡢 3 ROLLD 1 - 1 SWAP SUB » J IF OVER 1 == THEN SWAP DROP ELSE DUP SIZE DUP2 DUP SUB "Z" > - 1 > IF THEN "(" SWAP + ")" + END "_" + SWAP + END END 

3 routines packaged neatly in a directory. M is the main one. It expects 2 strings on the stack formatted as ions and pushes a molecule string onto the stack.

S splits the ion into charge as a number and the element formula as a string. For example, "PO_4^-3" would be taken from the stack and -3 and "PO_4" pushed onto the stack.

J joins the number of ions with the formula and decides whether to wrap the formula into brackets. The bit before ELSE deals with 1 ion, leaving the string as it is. For example, if 1 and "PO_4" are on the stack, they are replaced by "PO_4". 1 and "H" gives "H".

The rest deals with multiple ions; if it's a single atom it's not in brackets, otherwise it is. To decide whether it is, I check the length of the string and check if the last character is >"Z". Boolean expressions return 1 for true and 0 for false. By subtracting the result of this comparison from the length of the string, I get 1 or less when it's single atom, otherwise more: length 1 is a single atom; length 2 will have a letter as the last character; for a single atom it's lower case, so >"Z", making result 1, otherwise 2; length 3 or more means more than 1 atom and with 0 or 1 subtracted from the length the result will be at least 2. For example, 3 and "PO_4" gives "(PO_4)_3". 3 and "Al" gives "Al_3".

M first splits each ion using S. After the first line the level 5 of the stack (so the deepest buries object) contains charge of the second ion, level 4 second ion's formula, level 3 first ion's formula, level 2 absolute value of the first ion's charge and level 1 absolute value of the second ion's charge again. For example, if given ions on the stack are "Al^+3", "SO_4^-2", we get -2, "SO_4", "Al", 3, 2.

Second line calculates gcd of the 2 charges (leaving the charges intact).

Third line divides each charge by the gcd (to calculate the multiples) and joins it with the ion's formula using J. So we have two strings each with one given ion with charge removed (or a multiple of it) and a charge of the second one buried behind them. For example, -2, "Al_2", "(SO_4)_3" (-2 is the charge of SO_4).

Fourth line unburies the charge and if it's positive it swaps the two strings (so that the cation comes first) before joining them. So in example above, because it's negative, they are joined in the order as they are: "Al_2(SO_4)_3".

Python 3, 131 bytes

lambda c,a,C,A:g(c,A/gcd(C,A))+g(a,C/gcd(C,A)) from math import* g=lambda s,n:[s,[s,'(%s)'%s][sum(c<'`'for c in s)>1]+'_%d'%n][n>1] 

Try it online!

Python 2, 196 174 170 155 149 140 136 bytes

lambda c,a,C,A:g(c,A/gcd(C,A))+g(a,C/gcd(C,A)) from fractions import* g=lambda s,n:[s,[s,'(%s)'%s][sum(c<'`'for c in s)>1]+'_'+`n`][n>1] 

Try it online!

Python 3, 129 bytes

lambda E,e,I,i,m=lambda e,i:(len(e)>2<=i)*"("+e+(len(e)>2<=i)*")"+"_%d"%i*(i>1):m(E,i/gcd(i,I))+m(e,I/gcd(i,I)) from math import* 

Try it online!

If we need to handle negative charges of the anions, then 153 bytes:

lambda E,e,I,i,a=abs,m=lambda e,i:(len(e)>2<=i)*"("+e+(len(e)>2<=i)*")"+"_%d"%i*(i>1):m(E,a(i)/gcd(a(i),a(I)))+m(e,a(I)/gcd(a(i),a(I))) from math import* 

Try it online!

JavaScript, 145 bytes

(c,a,p,n,g=(a,b)=>b?g(b,a%b):a)=>`(${c})_${n/-g(p,-n)}(${a})_${p/g(p,-n)}`.replace(/\(((\w+)\)(?=_1)|([A-Z][a-z]*)\))/g,"$2$3").replace(/_1/g,"") 

Takes arguments c is cation, a is anion, p is positive charge, n is negative charge.