From idealised to predictive models of liquid crystals

ABSTRACT A brief review and personal perspective of modelling and simulations of liquid crystals and of how this has dramatically changed over the years, from idealised lattice to molecular level and now to predictive atomistic models. Graphical Abstract


Introduction
Upon approaching the Italian mandatory retirement age [1] for university professors it is difficult to resist the temptation of looking back at how the field of liquid crystals modelling and simulations I have been involved for over 40 years has developed, and radically changed, during this time span.
In the seventies, when the work about modelling the properties of liquid crystals started, the systems studied with computer simulations were normally rare gases and other so called simple fluids made of spherical particles [2,3], while the first pioneering molecular dynamics (MDs) simulations of water only started to appear in those years [4]. It is maybe worth recalling that the computational power at the time was embarrassingly weak for today standards. The top machine in the seventies was the CDC7600, with 36 Mflops peak power, hardly comparing to that of a current smartphone, which is of the order of a few hundred Mflops instead. Input to the machine had evolved from teletype paper tape to decks of punched cards and storage was also amazingly meagre for today standards. A RP07 removable disk unit of the size of a washing machine (and equally noisy) could store 500 MB of data, a feast of bytes for the time, while a current micro-SD card the size of a fingernail can approach 1 TB. Last, but not less annoying, printouts of results from individual runs had often to be transcribed by hand for subsequent averaging and analysis.
Facing this situation, the task of bottom up modelling and of trying to reproduce the properties of liquid crystals of mesogenic molecules as complex as those shown as examples in Figure 1, rather than 'simple liquids', seemed just impossible.
Indeed, computer simulations of liquid crystals started with a simple lattice model proposed by Lebwohl and Lasher (LL) [6], described later on, along the lines of the somewhat similar Ising and Heisenberg models for magnetic systems [7], more or less when I started my Ph.D. thesis with Geoffrey Luckhurst [8]. My thesis was actually on statistical theories of liquid crystals and on the simulations of Electron Spin Resonance (ESR) spectra for radical probes in liquid crystals, but I also started doing Monte Carlo (MC) simulations on the Lebwohl-Lasher system, and these first results were published in the wonderful 'green' book resulting from the lecture notes of the NATO ASI organised in Cambridge by Geoffrey Luckhurst and George Gray in 1977 [9]. I was still under 30 when I delivered those lectures [10,11] to an impressive audience (see Figure 2) and I still regard that School as the most important and stimulating I ever attended and remain grateful to Geoffrey for that daring invitation.
The models upgraded over the years from lattice to off-lattice ones, where molecules are replaced by relatively simple objects either taken as purely hard repulsive [12] or endowed with attractive as well as soft repulsive interactions, like the Gay-Berne (GB) [13], and, only recently, to rather realistic fully atomistic models [14].
In the next Sections I shall try to offer some personal reflections on the development of these models and on future perspectives.

Lattice models
Choosing a lattice model means treating a liquid crystal as a sort of plastic crystal, where molecular positions are restricted to a lattice and only orientations change, a choice that may seem crazy at the outset, looking at how very different this is from real liquid crystals, where fluidity is a defining feature. The solution of the paradox is that, when choosing such a lattice model, the aim is not to try and reproduce all of the properties of a real liquid crystal, but rather its orientational ones alone, particularly around the nematic-isotropic transition, where the behaviour is expected not to depend too much on molecular details. Some of the universal characteristics of the nematic-isotropic phase transition are: being of weak first order, with small entropy of transition and small volume change, both roughly 10% of the solid-isotropic liquid one [15], with an orientational scalar order parameter hP 2 i ¼ h3ðu Á nÞ 2 À 1i=2 (with u the molecule axis and n the director) which at the isotropic transition is hP 2 i NI % 0:3 À 0:4: hP 2 i is just the most important example, of rank L ¼ 2; of a complete set of order parameters defined as Legendre polynomial averages ðP L Þ [10]. The temperature variation of the order parameter in reduced units (T Ã ;T=T NI ) is also, to a first approximation, quite similar for different nematic, if not universal. The NI transition also shows strong pretransitional effects in Kerr or Cotton-Mouton experiments [16,17], similar to those expected from a second order transition, diverging at a temperature T y NI slightly below T NI , although differences exist for different materials [18], so that the phenomenon does not appear to be truly universal.
Going into some details, in the LL model molecules (or rather a tightly packed clusters of molecules) are represented by unit vectors, u i , placed at the sites of a simple cubic lattice and interacting with a pair potential where ij is a positive constant, , for neighbouring sites i and j, and zero otherwise. β ij is the relative orientation of the two particles (see Figure 3), i.e. cos β ij ¼ u i Á u j . This simple model has never been solved exactly, and chances are it will not be in the foreseeable future, like no other three dimensional lattice model has [19], but it has certainly been investigated by a number of authors using a great variety of theoretical techniques (see [20][21][22][23]). In 1985 we performed MC simulations of a 30 Â 30 Â 30 LL lattice with Periodic Boundary Conditions (PBC), and showed, by an analysis of energy and order parameter histograms, that it has a weak first order transition at T Ã NI ;k B T NI == 1.1232 AE 0.0006 [20]. When we did the simulations, the computational effort was significant and we could do it with Umberto Fabbri, at CINECA, the major Italian computer centre, only using some weeks of the burn-in time on the newly installed Cray XMP, the top supercomputer of those U ij ij u u Figure 3. (Colour online) The Lebwohl -Lasher interaction potential U ij between two particles i and j against the relative orientation angle β ij for ε ij = 1 (left) and a cubic lattice showing nearest neighbour spins u i , u j (right). days, with its 235 Mflops peak power (!). The transition properties, including the transition temperature, depend on sample size N; as shown in Figure 4 for smaller lattices, but the T Ã NI value has essentially been confirmed by other groups [24] also using larger, 60 Â 60 Â 60, lattices [25] and by other recent simulations [22]. Interestingly the model shows pretransitional effects diverging about one degree below the first order transition temperature, a behaviour entirely consistent with real experiments. Thus, albeit simple, the model possesses the fundamental aspects of the nematic ordering. It also features, although with the limitation of having only one elastic constant (like most often assumed in continuum type theories [16] anyway), the essentials of the director field and of its topological defects [26], allowing the predictions of defects in systems like droplets or thin films with specific boundary conditions.
When applied to simulating nematics confined to thin films [28] or to microdroplets embedded in some medium [35], i.e. situations where topological defects have to emerge, lattice model simulations provide on one hand results analogous to those obtained from continuum theory in relatively simple cases where both can be applied, while on the other being able of treating systems with complicated geometries or additional features. As an example, Nelson has showed in a seminal paper [39] that a thin shell of nematic spread over a sphere of some material implementing tangential anchoring of the nematic should show four defects, each with a topological charge 1/2. MC simulations of a lattice model of the system [40] do reproduce the effect, but also show that different defect configurations can be obtained in the presence of appropriate external fields ( Figure 5) [40] or distorting the shape of the supporting sphere to that of an ellipsoid [41]. MC simulations are indeed very useful to treat director distributions in a variety of situations applicable to real experiments for confined nematics, e.g. for nanostructured wave guides [38].   [39] and MC simulations [40] in the absence of external field (centre) or the two defects obtained by an applied uniaxial field (right) [40].
Even though lattice models have been and are enormously successful for orientational order parameter, director field and defects, they have a number of obvious limitations, e.g. in that they cannot handle smectics or discotics or rheological properties (see Figure 1).

Molecular models
The introduction of generic molecular models, where a molecule is represented with a rigid particle of anisotropic shape approaching that of the real mesogen, e.g. rodlike or disklike (see Fig. 6), and is endowed with translational and rotational freedom, was a major achievement. It allowed to study not only nematics and their orientational order, but also smectics and columnar discotics, and to observe their spontaneous formation upon changing temperature or density. However, even with the introduction of off-lattice models, theories and computer simulations of liquid crystals have made drastic and often contrasting assumptions on the essential features of models representing constituent mesogens and on their interactions.
Similarly to what happened for simple fluid made of spherical particles (rare gases for instance), two main classes of models have been developed. The first focuses on steric repulsions, using hard particles somewhat similar in shape to real mesogens but with no attractive forces [12]. The real surprising and important result, mainly due to Daan Frenkel and collaborators, was that nematic phases could be formed even in the absence of attractive interactions, just by increasing the particle concentration. As internal energy is absent for purely hard particles, this implies that the entropy of a nematic formed by hard spherocylinders [42] or ellipsoids [43,44] is higher than that of the corresponding isotropic phase. These examples can be very reminescent of (and applicable to) suspensions of elongated colloidal particles, even though one would imagine that in real cases the inevitable shape dependent dispersive attractions would be difficult to discount [45]. It also illustrates, like the simpler case of the crystallisation of systems of hard spheres, the difficulty in finding the microscopic organisation with the highest entropy (in the case of hard spheres at sufficiently high concentrations, the crystal [46]). This difficulty is also manifested for the further transition to smectics, occurring for hard spherocylinders, but not for hard ellipsoids. Simulations have now been performed for a variety of hard particle polyhedral shapes, particularly by Sharon Glotzer and her group [47].
Attractive-repulsive models for mesogens, perhaps more akin to chemical intuition and more appropriate to treat thermotropics, were pioneered by Bruce Berne and collaborators and in particular his GB model [13], providing an anisotropic version of the Lennard-Jones (LJ) model ubiquitous in describing simple fluids, has become a reference model for liquid crystals. In Figure 7 we see the GB potential between two uniaxial ellipsoidal particles for three important approach configurations, each with a LJ type behaviour but, of course, each with a different contact distance and well depth.
We started work on GB models with Andy Emerson, then a postdoc from Southampton, who afterward married a colleague and started working at CINECA, thus remaining to date in Bologna and with Roberto Berardi with whom, as well as, later, with younger colleagues like Silvia Orlandi, Matteo Ricci, Luca Muccioli, Lara Querciagrossa, I have continued collaborating to date. At that time we proposed a GB parameterisation that favours side-side interactions [48], thus extending the liquid crystal range (see Figures 7,8). With Roberto we also developed a biaxial version of the GB pair potential [49] and we showed that a parameterisation contrasting the strong face to face attraction leading to smectics and crystals could unequivocally allow the formation of biaxial nematics [50,51]. This confirmed that the 'elusive biaxial nematic' [52] could exist in a system of particles that had the choice to become smectic rather than nematic and not just biaxial rather than uniaxial nematic as already indicated long times before by mean field theory [53] and lattice MC simulations [27]. The formation of an off-lattice biaxial nematic was also demonstrated by Mike Allen [54] in a system of hard ellipsoids but, as already mentioned, this particular system does not form smectics [44], so one of the main competing phases is in that case absent.
Apart from predicting thermodynamic properties, molecular resolution models provide a playground for virtual experiments. For instance we simulated the switching of a pixel in a twisted nematic LC display [55] and the fast relaxation of the secondary director expected for a biaxial nematic device [56].
We also actively pursued two additional ways of generalising the GB model.
The first was to allow for different particle shapes and in particular, apart from biaxial [57], bowl like mesogens [58] with the aim of modelling pyramidic liquid crystals [59] and polyphilic [60] tapered shapes to try to obtain ferroelectric nematics [61].  The other route was to decorate the GB with reactive sites and put forward a procedure for obtaining polymers as chains of GB particles connected by flexible springs [62]. More recently, together with Gregor Skačej, we have also prepared main chain LC elastomers (LCE), using our polymerisation procedure to implement, as far as possible the experimental protocol of Finkelmann to obtain LCE. This has proved quite difficult and was only practically possible by employing a softer version of the GB pair potential that we had developed in the meantime with Mark Wilson, Roberto Berardi and Juho Lintuvuori [63]. The LCE work allowed us to perform virtual stress-strain MC experiments obtaining Young moduli and shedding some light on the mechanism of elongation-contraction leading to actuation driven thermally by going through the NI transition [64] or by the application of an electric field [65] in the swollen LCE systems experimentally prepared and studied by Urayama [66,67]. We performed stress-strain virtual experiments to examine the relation between mechanical properties and molecular structure and we arrived at proposing a mechanism for the phenomenon of supersoft elasticity [68], one of the most fascinating properties of some LCE, particularly the main chain ones polymerised in the isotropic phase .
Also for off-lattice models I had the privilege to collaborate over the years with many colleagues, apart from those already mentioned: Martin Bates, Reiner Memmer, Joachim Stelzer, Demetri Photinos, Alekos Vanakaras and Doug Cleaver.

Predictive atomistic models
While molecular models can be of great help in understanding generic trends and to answer questions concerning the effect of a specific feature, e.g. of the aspect ratio or of an electric dipole at a certain position and orientation [69,70] on the molecular organisation, or the shift in transition temperatures, the prediction of the existence or not of a liquid crystal phase is an entirely different and more demanding challenge. A major problem in this respect is the subtle, but possibly large, effect that changing details of molecular structure has on phase transition temperatures. This is well illustrated by the huge change in T NI in the cinnamates homologous series in Figure 9, where the introduction of 1, 2, . . ., n apparently innocuous methylene groups shows a major odd-even effect, with a decrease of some 250°C when introducing the first CH 2 group and an increase of over 100°C with the second, and so on. This is quite different from the typical effect of adding methylene groups in normal alkanes C n H 2n+ 2 where a regular stabilisation quantified by the linear increase of the standard vaporisation enthalpy, ΔH vap , with n occurs. The gloomy feeling on the predictive possibilities of atomistic MD simulations, until the end of the nineties is well described by this statement from [71]: 'Contemporary computer power is insufficient to reproduce unambiguously phase transitions and order parameters with realistic potentials'. The situation was, however, destined to improve, and we took up the challenge of trying to reproduce, at least approximately, the odd-even effect in cinnamates.
In Figure 10 we compare experimental [72] and simulated [73] values of TNI for the first three members of the series and we can see that the accord with experiment is reasonable, with the odd-even effect reproduced with fair accuracy, demonstrating the level of realism that could be reached with atomistic simulations. Simulations also indicated the origin of the odd even effect in the shape changes due to non°°°°°°F igure 9. (Colour online) The odd-even alternation of the N I transition temperature in the phenyl alkyl-4-(4ʹcyanobenzylidene) amino cinnamates series [72] and the first three homologues n = 0, 1 and 2 showing the odd-even shape change in the fully stretched conformation, We indicate with thick lines the inertial molecular axis (light grey) and the terminal phenyl para axis (dark grey). The torsional angles ϕ 1 , and ϕ 2 are also shown [73]. Adapted with permission from ref. [73]. rigidity of the molecules leading to a bimodal distribution of aspect ratios.
This first successful result preceded a vast effort [14,74] to extend simulations to the cyanobiphenyls series (nCB), possibly the most important and wellstudied family of liquid crystals. Indeed, following the bibliometrics fashion of our times, according to Web of Science 5CB is a topic for over 2300 papers, has nearly 35,000 citations and an h-index of 71 from 1985 to 2017, but also 8CB is doing well, with around 700 papers, 11,500 citations and h-index of 47. In Figure 11 we show our results for the nematic-isotropic transition temperatures of the n = 5-8 nCB, which differ from experiment by less than 4 degrees. This required a suitable tuning of the Force Field (FF), validated by comparison of a large series of predicted results for observables like birefringence, NMR, Raman, X-ray diffraction. In particular the most important order parameters hP 2 i and hP 4 i are plotted in Figure 12. It is worth noting that the results for hP 2 i are in good agreement with experiments, in the sense that they are not more different from typical experimental data than discrepancies observed between experimental data obtained, even from the same physical observables (e.g. birefringence) by different groups. Moreover, the predicted order parameter hP 4 i that appeared to be rather far from the much lower and even slightly negative values obtained a few years before in [75], is found to be in much better agreement with Raman work from Giesselmann group [76] appeared after publication of our simulations (see Figure 12). Considering 8CB instead, the first of the nCB series to yield a smectic phase, we obtained layer spacing d= 32.4 AE 0.2 Å [77], in good agreement with experimental X-ray studies from literature that give: d= 31.6Å [78]; d= 31.432Å [79], d= 31.Å [80].
While the agreement with experiment is exciting and promising, it is worth noting that Force Fields still represent the weak point, if not the fly in the ointment of the procedure. Indeed, tuning the FF is of crucial importance and using a standard off the shelf one, we have found, like others [85], that one can easily give T NI 70-100 K off from experiment. Various attempts at defining a general procedure to tune a FF  Figure 9right [73]. Adapted with permission from [43]. Figure 11. (Colour online) Experimental [9] (blue squares) and simulated (red squares) [14] nematic-isotropic transition temperatures for the nCB series. The error on the assignment of the simulated T NI is estimated to be 2 À 4 K at most. Adapted with permission from [14].  [14] as a function of T À T NI compared with experimental values from various techniques: birefringence: (a) [81], (b) [82] and (c) [83]; Depolarized Raman: hP 2 i (d) [75] and (e) [84] and hP 4 i (f) [75], (g) [76]. The blue lines are a guide for the eye, error bars reported in [14] are omitted here to reduce clatter. with permission from [14].
for liquid crystals have been put forward [86,87], but the goal still seems not so close.
To start with, it is important to stress that all types of intermolecular contributions should be included: so, e.g. allowing only for atomistically detailed shape and dispersions is not sufficient. For instance McDonald and Hanna [88], used atomistic MD to study the phase behaviour of 8CB with the Amber united-atom force field without electrostatics contributions and found T NI more than 100 K higher than experiment. Even including all relevant interactions, parameters have to be tuned. This has been observed and discussed for a series of alkenic fluoroterphenyls by Yan and Earl [89] who also provided a nice summary of other cases.
It may be worth mentioning that the use of different FF and different MD codes can lead to relatively high systematic errors in simulations. As shown by a recent round-robin study [90], results on simple alkanes can produce, e.g. differences in densities from 0.620 to 0.760 kg/m 3 for butane. Such density differences, perhaps acceptable for normal liquids, would lead to major errors in ordering for LC phases [91].
A general issue that atomistic simulations can tackle is to what extent real molecules can be approximated by simple rigid objects, like the molecular resolution models of the previous Section. A simple example is provided by quinquephenyl (P5) [91] which at least has the advantage of not having the flexible chains nearly universally present in mesogens (see Figure 1(a)). Ed Samulski and collaborators studied P5 experimentally and suggested that: If ever there were a calamitic mesogen that corresponded to the approximations used to derive S, the rod-like thermotropic LC PPPPP is one among them [92]. However, atomistic simulations [91] show that even P5 is not so rigid and straight as could be imagined, and that it shows a non negligible distribution of the P5 internal bending angle. Comparison with a hard spherocylinder of similar aspect ratios also shows considerable differences [91].
As a natural evolution of the study of bulk liquid crystals, in the last few years I have been very interested in the problem of predicting the alignment of LC on a given surface. Since most, if not all, LC devices are based on thin films, the direction of molecular alignment and the strength of anchoring of the LC on a given material is of obvious importance. However, this knowledge has been, to date, mainly empirical, rather than being predicted bottom-up from molecular principles. Having prepared and validated some in silico nematics, and in particular 5CB in the bulk, the natural question is thus: can we predict the surface anchoring of a given mesogen on a certain support surface by atomistic MDs simulations? And moreover, given a certain LC material, how important is the chemical nature of the surface? Having fixed a certain chemical composition of the surface, how important is the morphology of the surface and its roughness? To try and make predictions we need a well defined surface (chemical composition, morphology, roughness, . . .) and ideally we would like to try to build a library of different materials. With this in mind, we have now studied various types of support surfaces, as shown in Figure 13, either solid (left column) or 'soft' (right column).
As an example of the type of detailed information that can be obtained, we show in Figure 14 the scalar order parameter across two thin 5CB films (11 and 22 nm thick, respectively), either in the isotropic or nematic phase, when deposited on the (001) surface of a slab of H-terminated crystalline silicon [95]. From the snapshot in the lower part of Figure 14 we also see that the anchoring direction of the 5CB is planar at the silicon interface and homeotropic at the free surface, where we also have chains tending to point outside the film. We also see from Figure 14 (upper part) that the scalar order at the crystalline surface is very similar in all cases and higher than that inside the film. We have found a similar behaviour for other crystalline substrates and in particular for 5CB on cristobalite quartz silica [96]. Quite differently, 5CB on silica, but in the form of glass has an order at the solid surface overlayer lower than that in the middle of the film, showing the importance of surface morphology [96]. This suggests that changing by some treatment the crystallinity of the solid support materials could help in varying the order at the surface itself and possibly improve technological properties connected to it, e.g. the contrast ratio in some LC devices. A detailed analysis of the alignment on self-assembled monolayers (SAM) of different types grafted on a flat silica surface ( Figure 13) is also useful in understanding the conditions that lead to homeotropic alignment of the nematic, which we find favoured if the SAM is a mixture of long and short chain amphiphiles [97].
A few years ago, with the development of fairly reliable atomistic simulations, we started a stimulating collaboration, continuing to date, on ordered phases of interest for Organic Electronic with, amongst others, Yves Geerts at ULB, Brussels, David Beljonne and Jerome Cornil in Mons, Alison Walker at University of Bath. Our contribution has been and still is that of trying to obtain morphologies for organic functional materials, some of which liquid crystalline, e.g. phthalocyanines [98] or perylene derivatives [99] or triphenylenes [100] forming columnar molecular organisations which can be used as input for calculation of charge transport or, in the case of donor-acceptor (D-A) Figure 14. (Colour online) Local hP 2 i order parameter of 5CB as a function of distance z from a silicon Si(001):H support surface for a film % 11 nm thick (at T = 300 K, red, and T= 315 K, green) and for a % 22 nm thick one (black at 300 K and blue at 315 K). The local director changes from planar at the solid interface to homeotropic at the free surface, represented by the colour change from grey to blue as shown by the snapshot cut-out at the bottom [95]. Adapted from [95] by permission from Royal Society of Chemistry. Figure 13. (Colour online) On the left some 'hard' solid surfaces we have studied: hydrogen terminated silicon [95] (a), silica in a cristobalite quartz morpholgy (b) silica glass (c) [96], fullerene [93] (d) and on the right some 'soft' surfaces: self assembled monolayer (SAM) on: an alkylsilane (OTS) on silica glass [97] (e), a fluoroalkysilane (FDTS) SAM [97](f), and two polymers: PMMA [94] (g) and polystyrene [94](h). systems for estimating the chance of charge dissociation and recombination at the D-A interface in a model photovoltaic cell. In this context I should certainly mention that atomistic work has only been possible with the essential collaboration of Luca Muccioli, Roberto Berardi, Gabriele D'Avino, Giustiniano Tiberio, Antonio Pizzirusso, Adriana Pietropaolo, Yoann Olivier, Domenico Summa, Lara Querciagrossa, Otello Roscioni, Julien Idé and Mattia Palermo.

Conclusions and perspectives
On closing this brief overview of model development and of our contributions, it might be worth, as a true sign of old age, to close with some perspective views, if not words of wisdom.
(i) The increasing success of realistic atomistic simulations shows that the outlook for predictive modelling has to be clearly optimistic if computer performance, that have already increased by a factor of the order of 10 4 in the last 20 years [101], will continue to evolve from current petascale to hexascale and beyond. (ii) Increase of computer performance, while to be thankfully acknowledged, is not enough and, beyond the Force Field developments that we have already mentioned, software improvements are much in demand, in particular better techniques for time integration of Newton's equation of motion, that would allow increasing the elementary time step [102] and thus the overall time window for observation of a system. (iii) We have discussed models at different scales as separate, independent approaches, but a coarse graining procedure that allows going in a systematic way from the atomistic to a molecular level where molecules are represented, e.g. by suitably connected sets of GB beads would be highly desirable in perspective and starts to be feasible at least for certain classes of organic functional materials employed in organic electronics [103]. (iv) While we have been concerned with performing simulations of equilibrium phases, guided by a free energy minimisation principle, many practical fabrication techniques involve nonequilibrium processes, e.g. evaporation, vapour deposition, casting and so on. A huge amount of work remains to be done in this area, both from the methodological and predictive point of views.
A persisting problem in predictive simulation is the need for extensive supercomputing resources, which even though more widespread may be not always available when needed. Thus the suggestion reported by Löwdin back in 1978 still holds: If you do not have access to such big computers, your only chance is to follow the advice from Peter Debye's 1955: 'Beat them' i.e. develop theoretical methods so powerful that you do not need so much computational hardware [104], even if this might be rather difficult to do.

Acknowledgments
I would like to thank at least some of the many friends and colleagues in Bologna that I have not already cited in the text, and that, over the years, have shared my passion for liquid crystals, computer simulations and the spectroscopic techniques, like ESR, NMR, Fluorescence Depolarization, used to study them, and in particular: Alberto Arcioni, Corrado Bacchiocchi, Sigismondo Boschi, Carlo Fava, Renato Ferreira, Manuele Lamarra, Isabella Miglioli, Franco Semeria, Francesco Spinozzi, Sai Preeti Gouripeddi, Riccardo Tarroni, Ilaria Vecchi.
I am grateful to EU for funding a number of collaborative projects with leading European Groups and in particular for the support from grant EXTMOS (EXTended Model of Organic Semiconductors) EU-H2020-646176. I also wish to thank Italian MIUR for PRIN 2015XJA9NT (Molecular Organization in Organic Thin Films via Computer Simulation of their Fabrication Processes) for financial support. Last but not least, I would like to take this opportunity to thank my wife Nicoletta for her patience and support while I was playing with molecules and liquid crystals.

Disclosure statement
No potential conflict of interest was reported by the author.