Can we optimize a molecular structure using Moller-Plesset MP4 method on GAUSSIAN software? If yes what keyword do we have to add?
Gaussian can do all the way up to MP5 for single-point energies, but analytic gradients are only available for MP2, MP3, MP4(DQ) and MP4(SDQ), the latter including only single, double and quadruple substitutions.
MP4(SDTQ) which is like MP4(SDQ) but also with triple substitutions, and MP5, are listed in the manual as being able to do numerical gradients. Keep in mind that the Gaussian documentation already linked above says "MP4: Fourth-order MP theory correction [Raghavachari78], which defaults to full MP4 with single, double, triple and quadruple substitutions [Raghavachari78, Raghavachari80] (MP4(SDTQ))", meaning that MP4(SDTQ) and MP4 mean the same thing in Gaussian's terminology.
All of what I wrote above applies only to RHF and/or UHF reference states, since Gaussian only provides single-point energies for MP methods with an ROHF reference (you can see some discussion about that at the Chemistry SE here: Is it possible to perform ROMP2 numerical optimizations in Gaussian09?).
Beware that MP4 calculations are quite a lot slower than MP2 ones, which can make geometry optimization quite prohibitive depending on the size of your system, so you would probably want to use analytic gradients with RHF or UHF (instead of ROHF) and with MP4(DQ) or MP4(SDQ) rather than trying to do numericaly gradients in Gaussian with MP4 or numerical gradients manually on your own using ROMP4.
Furthermore, the MP series (MP2, MP3, MP4, MP5, etc.) is rarely used beyond MP2, largely because methods like CCSD can be far more accurate, and many DFT functionals may be more accurate while also being less computationally expensive (especially double-hybrid functionals which combine DFT with MP2-style calculations and can be astonishingly close to CCSD(T) in accuracy while being only a bit more expensive than MP2, see my answer to this: What are some recent developments in density functional theory?). So if you're about to launch 100s of geometry optimization calculations at the MP3 or MP4 level, I would encourage you to consider carefully the purpose of the project.
Complementing the @Nike answer:
Can we optimize a molecular structure using Moller-Plesset MP4 method on GAUSSIAN software?
If yes what keyword do we have to add?
You need to add the corresponding keywords like below:
# opt mp4(sdtq)/6-31++g(d,p) geom=connectivity
This keyword requests that a geometry optimization be performed. The geometry will be adjusted until a stationary point on the potential surface is found. Analytic gradients will be used if available. For the Hartree-Fock, CIS, MP2, MP3, MP4(SDQ), CID, CISD, CCD, CCSD, QCISD, BD, CASSCF, and all DFT and semi-empirical methods, the default algorithm for both minimizations (optimizations to a local minimum) and optimizations to transition states and higher-order saddle points is the Berny algorithm using GEDIIS [Li06] in redundant internal coordinates [Pulay79, Fogarasi92, Pulay92, Baker93, Peng93, Peng96] (corresponding to the Redundant option).
The MPn method keywords request a Hartree-Fock calculation (by default, RHF for singlets, UHF for higher multiplicities) followed by a Møller-Plesset correlation energy correction (...) MP4: Fourth-order MP theory correction, which defaults to full MP4 with single, double, triple and quadruple substitutions (MP4(SDTQ)).
6-31++g(d,p): Is the chosen basis set. You need to select the basis set that better represent your system.
The Geom keyword specifies the source of the molecule specification input, options related to coordinate definitions, and geometry-related output. By default, it is read from the input stream, as described previously. Geom may be used to specify an alternate input source. It also controls what geometry-related information is printed and use of internal consistency checks on the Z-matrix. The Geom keyword is not meaningful without at least one item selection option.