Density Functional Tight Binding Method is 2-3 orders of magnitude faster than typical DFT but with accuracy approaching that of DFT for many systems including graphene. You can do all things with DFTB which you can with DFT including electronic structure, geometry optimization, MD, electronic transport, etc. Besides accuracy the advantage w.r.t. other semiempirical methods is that it needs minimal empirical parameters which you would not need to fit yourself as they already exist for carbon. One of the best implementations of DFTB is DFTB+ software which you can find here: http://www.dftb-plus.info/

Hi. Obviously, molecular dynamics (MD) simulation is quite suitable for your case.

You can find a lot of publications on it.

Actually, coronene molecule can serve as the smallest possible graphene model in many calculations. It consists of 24 carbon atoms. Also, there are many works in which graphene models consist of 100 carbon atoms or fewer were studied. If you wish to be on the "safe side" of calculations, you may use periodic boundary calculations using VASP or CP2K programs. Using the programs you can involve in your study infinite sheets of graphene.

You might want to take a look at the SIESTA code http://icmab.cat/leem/siesta/ to see if it fits your needs for this project. It uses localized basis sets resulting in linear-scaling algorithms. You can effectively tune accuracy-vs-cost with SIESTA to facilitate treating very large systems composed of thousands of atoms. - On a side-note: we have worked with interaction of biomolecules (DNA/RNA nucleobases and amino acids) with 2-D materials like graphene and h-BN in the recent (http://dx.doi.org/10.1021/jp402403f) and not so recent (http://dx.doi.org/10.1103/PhysRevB.76.033401) past . But our simulated systems were very small compared to what you have apparently in mind.

I'm curious about your need for greater than 1000 atoms - I would imagine that a smaller system on the order of 100 atoms could capture the important interactions with organic molecules unless you are trying to treat disorder in your systems. Both SIESTA and Quantum Espresso are free codes that could treat 100 atom systems using DFT. Otherwise, free molecular dynamics codes such as DLPOLY or Lammps would work for a larger 1000 atom system.

With 1000 atoms you need use the DFTB+ http://www.dftb-plus.info/

Density Functional Tight Binding Method is 2-3 orders of magnitude faster than typical DFT but with accuracy approaching that of DFT for many systems including graphene. You can do all things with DFTB which you can with DFT including electronic structure, geometry optimization, MD, electronic transport, etc. Besides accuracy the advantage w.r.t. other semiempirical methods is that it needs minimal empirical parameters which you would not need to fit yourself as they already exist for carbon. One of the best implementations of DFTB is DFTB+ software which you can find here: http://www.dftb-plus.info/

Couldn't you try to understand the behaviour of graphene (don't know what you are working on) and use a modified aproach of the methods listed above by the others?

Any molecular dynamics code like Gromacs, LAMMPS or NAMD should do the work. For what you want. NAMD in particular (see here for details http://www.ks.uiuc.edu/Research/namd/ ) is probably a good choice as it has been designed for large molecules. Nevertheless I would suggest you try to start from a small model and use DFT-like methods as Jonathan has mentioned earlier. If the scalling is not accurate enough for your system then Molecular Dynamics is probably your only option.

Out of curiosity what sort of system are you trying to study? How big are the organic molecules so that you need so many atoms?

I think that coronene is not necessarily a model of graphene surfaces. Before using Molecular Dynamics make a graph of the electronic chemical potential versus several sizes of graohene models. Please, see “Minimal molecular models for the study of nanostructures” and "Short nanotubes or big molecules?" posted in my profile, Thxs.

Another thing: experience has shown that when you reach an enough size you can use molecular mechanics (OPLS level) to fully optimize the geometry because nothing new will appear in the final geometry. (This is for the graph).

If your organic molecules are rather small, I think you can use DFT as mentioned above. For inspiration you can see "Adsorption of Small Organic Molecules on Graphene" ( http://pubs.acs.org/doi/abs/10.1021/ja403162r ). We used both finite model (coronene; "benchmark" calculations - CCSD(T), DFT-SAPT, etc. - allow us to choose suitable density functional (vdW-DF)) and infinite model of graphene (to perform ab initio molecular dynamics). Obtained adsorption enthalpies were in excellent agreement with the experimental data. OPLS-AA force field was able to recognize weakly and strongly bound molecules, and thus may be suitable for semi-quantitative estimates of the interaction energies of large molecules with graphene.

As many previous replies have pointed out, there are many options depending on what properties you want to calculate. Typical methods are: for macro length scale, there is Finite Element Method (FEM); for micro to nano length scale, there are Molecular Dynamics (MD) and Molecular Mechanics (MM), etc.; for nano to angstrom length scale, there are Tight Binding (TB) or Density Functional Theory (DFT). For your system of ~1000 atoms, I would suggest MD for mechanical properties and TB for electronic properties. Lammps, Gromacs (good for organic molecule compuation), VASP and Material studios are good software choices. Siesta is DFT package and 1000 atom system is bit too large for it. DFT is accurate but slow, MD usually ignores or approximates electrons but relatively fast.

FYI, for MD and TB examples, you can check some of my publications at http://scholar.google.com/citations?user=xcfH90cAAAAJ&hl=en

One of my papers (Nano Lett., 2013, 13 (6), pp 2692–2697) combined MM and TB methods for a system of ~6000 atoms using Lammps and own TB code. Hope it helps.

You'll probably need to combine a Tight-binding method with a dft calculation for your purpose. In order to obtain spectral function I suggest you to use a recursive green function approach or the kernel polynomial expansion method, which are the commonly used for big graphene. At least to my knowledge. Regards

You did not mention what is the initial state of graphene - pristine or defective? Are you going to apply periodic boundary conditions? Is your model of interaction includes orientational dependence? What is the ultimate goal of your computations?

The answers to these questions will navigate your choice of a code or a combination of codes. Most of the well known codes have already been listed. However, 10 cpus is a rather limited resource for a system with 1000 and more atoms, unless you have all the time of the Universe ahead :) Regards

Thanks for the suggestions! DFTB+ seems to be a good choice, as far as the system fits available parametrization. I have already used MM (OPLS-AA) and it worked very nice in many cases.

We have found that semiempirical PM6-DH2 method is also reliable for graphene systems.

Well if you have unit cells between 500 and 1000 atoms, you can use the fireball code, where there is a nice implementation for vdW calculations (LCAO-S^2+vdW)

We wrote codes by hands using numerical Kubo-Greenwood-formalism-based calculations, which are especially suited to treat large graphene samples containing millions of atoms, i.e. with dimensions approaching realistic systems. For example, in our recent "graphene" papers we calculated transport properties (diffusivity, conductivity) of graphene sheets containing 1700000 and 6800000 atoms, corresponding to the sizes of 210*210 and 420*420 nm^2, respectively.

Our experience of modelling graphene systems with PM6-DH2, perhaps might be usefull.

http://dx.doi.org/10.1039/c3cp53189a

FYI: DFTB+ is implemented into Materials Studio Program. It has a nice graphical interface and allows you create input files easily.

I've just found by chance the following paper (2013) - you might be interested in their DFT results combined with experiments:

Adsorption of Small Organic Molecules on Graphene; J. Am. Chem. Soc., 2013, 135 (16), pp 6372–6377

The strength of interaction between graphene and the organic molecules was estimated by density functional theory (PBE, B97D, M06-2X, and optB88-vdW), wave function theory (MP2, SCS(MI)-MP2, MP2.5, MP2.X, and CCSD(T)), and empirical calculations (OPLS-AA) using two graphene models: coronene and infinite graphene.