Axelle Vatrican

Evidentiality and Epistemic Modality in the Rumor / Journalistic Conditional in Spanish

Keywords: conditional, epistemic, journalistic, rumor, evidential

The aim in this paper is to explore the functions of the so-called rumor / journalistic conditional in Spanish. In order to do this, I will try to account for the fact that rumor conditionals are epistemic as well as evidential whereas conjecture conditionals in Spanish are epistemic but not evidential, as they convey uncertainty and do not encode the source of information. I will claim that the morphological marker -ría is a modal epistemic operator of possibility in both cases. In a rumor conditional, the epistemic operator quantifies over the illocutionary force of an embedded proposition p (enunciation / truth of the information), which means "maybe the information about p is true". The situation p is anchored in the present or in the future. In a conjecture conditional, the modal epistemic operator quantifies over the realization of an embedded proposition p (fact), which means "maybe the realization of p is true". The situation p is anchored in the past.

Introduction

The Spanish language has four modal uses of the conditional: a hypothetical conditional as in (1), a conditional of politeness or mitigation as shown in (2), a conjecture conditional as in (3) and a rumor or journalistic conditional as in (4):

(1) Si fuera rica, viajaría todo el tiempo.

'If I were rich, I would travel all the time'.

(2) Querría pedirte un favor.

'I would like to ask you a favor'.

(3) Juan tendría unos 40 años en aquella época.

'Juan must have been in his forties at that time'.

(4) Todos los españoles estarían a salvo menos uno, el mallorquín Guillermo Gual del 68. ("Abandonados y enfadados", El País, 15/01/12) 'It is said/ it seems that all the Spaniards are safe, except for one person the Majorcan Gual of the 68'.

(pre-print) Axelle Vatrican, (2015) Belgian Journal of Linguistics, 29, p. 83-100

The Rumor Conditional, in (4), which is also called "journalistic" conditional [de información no asegurada [START_REF] Lapesa | Tendencias y problemas actuales de la lingüística española[END_REF], periodístico [START_REF] Veiga | Condicionales, concesivas y modo verbal en español (Verba Anexo 34)[END_REF], del rumor [START_REF] Casado | El lenguaje de los medios de comunicación[END_REF]NGLE 2009)] has two basic features: the first one is evidentiality and the second one is epistemicity. On the one hand, the rumor conditional is evidential because it expresses indirect evidentiality: a non-firsthand knowledge, a type of reportative evidence or "hearsay" [START_REF] Willet | A cross-linguistic survey of the grammaticalization of evidentiality[END_REF]Cruschina & Remberger 2008: 97;Frawley 1992: 413) as shown in ( 5) and ( 6):

(5) On the other hand, a rumor conditional is epistemic as it expresses probability / uncertainty, as shown in the following examples, in which (b) is a possible equivalent of (a): the conditional can be understood as a future form combined with the adverb supposedly (Sp.

supuestamente). The speaker introduces a new and doubtful information ("Sevilla tendrá carrera de psicología"). As we will show, this conditional form conveys uncertainty with respect to the truth of the information and not with respect to the fact that Sevilla offers a degree in psychology.

(7) a. Sevilla tendría próximamente carrera de psicología para los jóvenes que deseen estudiar esta carrera en el municipio. (Canal Vivavisión, 08/07/2014) b. Sevilla supuestamente tendrá carrera de psicología para los jóvenes que deseen estudiar esta carrera en el municipio.

'Supposedly, Sevilla will offer a degree in psychology'.

(8) a. Two comments need to be made. Authors like Lapesa (1977: 227), Veiga (1991: 51;140) (11) Juan tendría entonces 40 años.

[Conjecture Conditional] 'Juan probably was in his forties / Juan must have been in his forties at that time'.

Although the same verb tendría is used in ( 10) and ( 11), a rumor conditional is evidential whereas a conjecture conditional is not, since in (12), tendría cannot be substituted by it is said that.

(12) a. Juan tendría entonces 40 años.

b. # 2 Según dicen, Juan tenía entonces 40 años.

# 'It is said that Juan was in his forties.' # 'Someone told me that Juan was in his forties'.

Since we have shed light on this contrast, two questions arise: why are rumor conditionals evidential whereas conjecture conditionals are not? Our aim is to account for the difference between the two conditionals. Our proposal is that the morphology -ría is a modal epistemic operator which operates on a situation p; this situation p is anchored in the present or in the future when we use a journalistic conditional; the situation p is anchored in the past when we use a conjecture conditional.

Description

In this section, we will underscore two features of rumor conditional in Spanish. The first one is a cross-linguistic one: the conditional is an evidentiality strategy and the second one is typical of Spanish: rumor conditional is a controversial form.

Rumor Conditional as an "evidentiality strategy"

The rumor or journalistic conditional exists in many Romance Languages: Spanish, French, Italian, Portuguese [START_REF] Squartini | The internal structure of evidentiality in Romance[END_REF][START_REF] Kronning | El condicional epistémico de atribución en francés, italiano y español : aspectos diafásicos, diatópicos y diacrónicos[END_REF]. It is usually described as an "evidentiality strategy" (Aikhenvald 2004, 105 & 276). This means that Spanish uses a verbal form to express the source of information ("evidential"). Thus, in a certain context, the conditional acquires an evidential meaning, although it is not always evidential (unlike dizque in American Spanish, for example, which is always evidential: Demonte & Soriano 2013). To express non-firsthand knowledge, it is well known that many languages use tenses like future ("unwitnessed and of subsequent lack of certainty"), non-indicative moods and modalities (South Estonian marker: the potential mood marker -ne), modal verbs (sollen. Germ. deber.

Sp.) whereas "declarative and assertive markers may give rise to direct and visual evidentials" (Aikhenvald 2004, 106-107, 277 ;Wiemer 2010, 65;81). The conditional is one of these "strategies", but the rumor conditional in Spanish has for a long time beenbeen considered as a controversial form. "The use of the conditional in that kind of sentences is strictly forbidden in the newspaper. In addition to the ungrammaticality, it minimizes the credibility of the information". This is a rumor, not an information, although it is not an ungrammatical form according to Grijelmo (2012):

Rumor Conditional in

Y esos "habría" "estaría" o "sería" no hacen sino narrar algo que el periodista no puede dar por seguro. Es decir, algo que no se ha comprobado o que no se puede atribuir a una fuente que lo asuma como cierto. Un rumor. (Grijelmo, 29/07/2012 El País, "el rumor enmascarado") "And those words as "should" "would be" or "could" do only narrate what the journalist cannot assure. That is, something that has not been proven or that cannot be attributed to a source that assumes this information is true. A rumor."

and Fundeu (2007):

El condicional de rumor no es gramaticalmente incorrecto, pero los manuales de estilo periodístico lo consideran inadecuado porque presenta un rumor como noticia.

[…] como según fuentes consultadas, posiblemente, cabe la posibilidad de que es posible, es posible, se cree que.

"The rumor Conditional is not grammatically incorrect, but journalistic style manuals consider it inappropriate because it introduces a rumor as news.

[…] and according to sources, possibly, it is possible that it is possible, it is believed that".

(pre-print) Axelle Vatrican, (2015) Belgian Journal of Linguistics, 29, p. 83-100

However, the Nueva Gramática de la Lengua Española -ASALE ( 2009) allows it as a form which introduces information in a very careful way:

[…] la variante del condicional de conjetura que se usa a menudo para presentar las informaciones de forma cautelosa (NGLE 2009, §23.15m).

" […] the variation of the conjecture conditional which is often used to introduce the information in a very careful way"

Rumor Conditional in Spanish: Epistemicity and Evidentiality

In this section, we will expose why a rumor conditional can be considered as an evidential / epistemic form in Spanish and which problems this gives rise to.

Previous analysis

The Rumor Conditional is usually described as following [START_REF] Abouda | Les emplois journalistique, polémique et atténuatif du conditionnel. Un traitement unitaire[END_REF]): 1. Firstly, it expresses uncertainty concerning the information conveyed. 2. Secondly, it indicates that the information has been taken from some other source. 3. Thirdly, the speaker / writer takes no responsibility for the information (Lapesa 1977: 227;Gili Gaya 1943: 168;[START_REF] Veiga | Condicionales, concesivas y modo verbal en español (Verba Anexo 34)[END_REF]Dendale & Tasmowski 2001: 339;Abouda 2001: 279;Gosselin 2001: 46;[START_REF] Haillet | Le conditionel en francais: une approche polyphonique[END_REF]Liddicoat 1997). Let us now look in more detail at these three dimensions, which as will become clear, cannot easily be combined with each other. 'Vicky Xipolitakis supposedly will travel to Dubai'.

The speaker introduces an information viewed as uncertain. This means that there is uncertainty because the Speaker repeats some information, acquired through another source without knowing if it is true or not. The conditional is considered epistemic because the possibility operates on a proposition [START_REF] Lyons | Sémantique linguistique[END_REF]NGLE 2009;Faller 2006, 29) The Rumor Conditional differs from other evidential forms in that it conveys both an epistemic and an evidential meaning. By contrast, it is well known that some evidential forms do not convey a hypothetical ("doubt") meaning, e.g. "que" in (18), which is purely evidential:

(18) Oye, que Manolo puso la bandera.

= 'Listen, (that) Manolo put the flag' # 'Manolo probably put the flag' (Demonte & Soriano 2013: 2).

2. It indicates that the information has been taken from some another source. In the example (19), we need to distinguish three elements: a situation p, the speaker 1 and the speaker 2:

(19) Vicky Xipolitakis viajaría a Dubai para ver a Maradona.

The situation p, i.e. the embedded proposition, is "Vicky Xipolitakis viajará a Dubai".

Speaker 1 says ( 19), which means that he assumes the hypothesis about this situation p, that is, "V.X. viajará a Dubai" is uncertain information. The epistemic hypothetical point of view comes from the Speaker 1. The speaker 1 is only saying: maybe what another speaker (the original source, speaker 2) has told me is true. However the source of the situation p, is another speaker, Speaker 2. Speaker 2 only asserts p: "V.X. viajará a Dubai", because as we said it is a reported speech, as if the speaker would have deleted "alguien me dijo que Vicky / someone told me that Vicky". We can give the following schema:

Speaker 1 (current speaker): [POSSIB [p]] (asserts a possibility about the truth of the information about a situation p)

Speaker 2 (original source) = p (asserts a situation p)

3. The speaker takes no responsibility for the information. It means that the Speaker 1 does not assert the embedded proposition p 'Vicky Xipolitakis viajará a dubai'. As a consequence, the speaker 1 does not commit to the truth value of p. He can say (20) if he is convinced that "V.X. will travel to Dubai" is true, as well as if he is convinced that it is false ("modalisation zéro", Kronning 2002, 567).

(20) Vicky Xipolitakis viajaría a Dubai para ver a Maradona.

(pre-print) Axelle Vatrican, (2015) Belgian Journal of Linguistics, 29, p. 83-100 = 'V.X. supposedly will travel to Dubai to meet with Maradona {and I think it is true / but I think it is false'}.

Problems

One of the major problems is the relation between epistemic modality and evidentiality. These are the most frequent questions: is the probability, i.e. epistemic modality, an inference derived from the evidentiality or is the epistemic modality a primary semantic feature of the conditional? Does the reported speech (hearsay) of the situation p prevent Speaker 1 from committing to the truth value of p? Does the fact that a speaker does not commit to the truth value of the situation p entail that we make an inference: someone else reported the situation p? According to the majority of authors, evidentiality is a primary feature. In the literature, one can find several elements to support this view: "the kind of evidence a speaker has will often determine the degree of certainty with which she believes a proposition" (Faller 2006, 2). Hearsay or inference convey uncertainty (Wiemer 2010, 85). Conditionals are basically evidential [START_REF] Dendale | Les problèmes linguistiques du conditionnel français[END_REF][START_REF] Dendale | Il serait à Paris en ce moment. Serait-il à Paris ? A propos de deux emplois épistémiques du conditionnel. Grammaire, syntaxe, sémantique[END_REF]. All tenses in Spanish have a primary evidential meaning [START_REF] Bermúdez | Evidencialidad. La codificación lingüística del punto de vista[END_REF]). Yet, the problem in Spanish is how to explain, as we said above (example ( 12)) that the morphology -ría only expresses epistemic modality in a conjecture conditional (21) whereas the same morphology -ría expresses epistemic modality together with evidentiality in rumor conditional ( 22) and ( 23).

(21) Juan tendría unos 40 años.

(22) Sevilla tendría próximamente carrera de psicología para los jóvenes que deseen estudiar esta carrera en el municipio.

(23) Todos los españoles estarían a salvo.

A very intuitive answer consists in making a difference between the Possible World and the Actual World (modal conditionals). On one hand, the possible World is the world where the situation p is located [Juan tener 40 años], [Sevilla tener carrera] or [V. X. viajar a Dubai];

if Speaker 2 asserts p, he will also be probably situated in this Possible W. The source of the information-Evidentiality-could be this speaker 2 situated in a possible world. On the other hand, the Actual W is the world where Speaker 1 is located: what Speaker 1 thinks at the present speech time (actual world) about a situation asserted by another speaker which is Speaker 2 (possible world). We will try to account for it in a more formal way.

(pre-print) Axelle Vatrican, (2015) Belgian Journal of Linguistics, 29, p. 83-100

Our proposal: [-ría] is a modal epistemic operator

In this section, we will propose a model of a semantic analysis of -ría morphology as an epistemic modal operator. First, we will present the schema of our proposal; then we will try to demonstrate that the embedded proposition p needs to be anchored in the present or in the future; we will explain that the modal epistemic operator operates on the illocutionary force of p. At the end, we will observe that rumor conditional, most of the time, appears with a third person and that using a first person would make difficult to interpret it as a rumor conditional.

Finally, to complete our analysis, we will compare the rumor conditional with the politeness conditional.

Modal conditional -ría as a modal operator

Our hypothesis is that [-ría] is an epistemic operator which operates on a situation p expressed in the present or in the future in a rumor conditional (Kearns, 2011: 185;[START_REF] Laca | Tiempo, aspecto y la interpretación de los verbos modales en español[END_REF]Escandell-Vidal 2004, 244;251;[START_REF] Mccready | Modal subordination in Japanese: Dynamics and evidentiality[END_REF]. We will give an example and make some comment. When we assert It's raining, this assertion can be true in two ways: it means that what I have said is true, i.e. the information of being raining is true: we are talking about the illocutionary force; or it means that at the speech time, it is raining actually: in that case we are talking about the fact of being raining. This difference will allow us to distinguish the rumor conditional from the conjecture conditional. Let us explore the first one. We will claim that a rumor conditional contains a modal operator -ría [POSS] which operates on the illocutionary force of a situation p; this situation p must be anchored in the present or in the future. In such a case we are talking about the possible truth of the information, i.e, of what someone said about a present or a future situation. We will give the following schema:

[POSS [to be true (p) [PRESENT / FUTURE]]]

The Tense focus of the situation p is the present or the future

The tense focus of the situation p is the present in (24)(25) On the opposite, the tense focus is the past with a conjecture conditional but neither a present nor a future as shows the contrast in (29):

(29) a.En aquella época, Juan tendría unos 40 años.

b.Juan {tenía / *tiene / *va a tener} probablemente unos 40 años.

'Juan probably {was/ *is / *will be} in his forties'.

(pre-print) Axelle Vatrican, (2015) Belgian Journal of Linguistics, 29, p. 83-100

The epistemic possibility operator and the situation p

When we talk about epistemic possibility, we make the assumption that a situation p is possible or probable according to the speaker's knowledge. However, no correlation exists between the presence of reportative and the degree of certainty with which the speaker believes the embedded proposition. Speaker 1, i.e. the speaker who asserts the whole sentence, can be convinced that the situation p is true or that the situation p is false. This is a feature which is shared by rumor conditional and conjecture conditional, since both examples are compatible with I am not sure of that as shown in ( 30) and ( 31).

(30) Vicky Xipolitakis viajaría a Dubai para ver a Maradona.

'It is said that V.X. will travel to Dubai but I am not sure of that [the information].

(31) Juan tendría en aquella época unos 40 años.

'Juan was probably in his forties but I am not sure of that

[the fact].'
The question is what we are talking about when we say "I am not sure of that". When we are talking about the truth-value of p, we need to distinguish between on the one hand, the commitment to the truth value of the rumor of p, i.e. the commitment to what has been said (the information), and on the other hand, the commitment of the truth value of the realization (the fact) of p, the commitment to the fact of p. In a rumor conditional, the speaker cannot commit to the truth of the 'rumor' of p, he cannot assert that what has been said is true.

Hence, in a rumor conditional the operator of possibility operates on the rumor that p: maybe the rumor that p is true as in (32), where (32a) is an equivalent of (32b) but not exactly an equivalent of (32c):

(32) a.Sevilla tendría próximamente carrera de psicología para los jóvenes que deseen estudiar esta carrera en el municipio. (Canal Vivavisión, 08/07/2014) b. = Maybe the rumor that in the future, Sevilla will supposedly offer a degree in psychology, is true. [information] c. # Maybe in the future (it is true that) Sevilla will offer a degree.

[fact] (pre-print) Axelle Vatrican, (2015) Belgian Journal of Linguistics,29, On the other hand, in a conjecture conditional, the speaker asserts that maybe the situation of p is true, which means that to the speaker, the situation p may have occurred already. We are talking about the realization of the embedded proposition p: in this case, p is a fact. (33a) is an equivalent of (33c) but not an equivalent of (33b).

(33) a. En aquella época, Juan tendría unos 40 años. b. #Maybe the rumor that at the time Juan was in his forties is true.

[information] c. =Maybe (it is true that) he was in his forties. [fact] In short, in a rumor conditional, we are talking about the possibility of the truth of what has been said. In a conjecture conditional, we are talking about the possibility of the realization of a situation. This could be schematized as following:

The rumor (information) that p = Rumor Conditional (p= present or future)

The realization (fact) of p = Conjecture Conditional (p= past) Another question could be: is the rumor conditional a kind of conditional construction <if p, then q>? In a certain sense, it is. It is a conditional construction but a 'logical' one or 'epistemic' as defined by Sweetser (1990: 128): "if I know that p is true, then I conclude that q is true". We have to distinguish the truth-value of the proposition (true/false) and the commitment with respect to the information (asserted by the speaker or by another person: "modalisation zéro", [START_REF] Kronning | Le conditionnel 'journalistique': médiation et modalisation épistémiques[END_REF]. A rumor conditional then depends on a protasis which is linked to the truth of what the speaker said:

(34) Sevilla tendría próximamente carrera de psicología para los jóvenes que deseen estudiar esta carrera en el municipio, si {fuese/es} cierta esta información 'Supposedly in the future, Sevilla would offer a degree in psychology, if this information is confirmed'.

A rumor conditional cannot depend on a protasis which expresses a situation which would have to be realized. With a protasis which refers to a fact, the rumor meaning disappears as shown in ( 35):

(35) Sevilla tendría próximamente carrera de psicología para los jóvenes que deseen estudiar esta carrera en el municipio, si lo permitiera el (36) Sevilla tendría próximamente carrera de psicología para los jóvenes que deseen estudiar esta carrera en el municipio, si fuese/es cierta esta información ''Supposedly, in the future, Sevilla would offer a degree in psychology if this information is/were confirmed'.

In a conjecture conditional, the conditional construction corresponds to: <q is true if p is true>, i.e. <q is true, if p is true>. The speaker is not sure that p occurred but he is sure that the unique source of information is him or herself. It is hypothetical because of him or herself, not because of another speaker's assertion.

(37) En aquella época, Juan tendría unos 40, {si no me equivoco/*si fuese/es cierta esta información}.

'At that time, John was probably in his forties {if I am not mistaken/ *if this information is / were true'. This can be explained by the fact that a double reading of conditional in Spanish is allowed.

Conditional actually has two temporal anchors and two interpretations: we can read it with respect to the situation time: the situation p is situated in a possible world, and the realization of p is located in a possible world as in <if p, then q>; we need to situate p in a possible world in a conditional construction; but we can also read it with respect to the present speech time.

This happens with a politeness conditional whose meaning is derived from a present implicature as in I would be happy if you came, which means in my present time, "I am asking you to come". In a rumor conditional, q is true if the information is confirmed; in a conjecture conditional, q is true if I am not mistaken. The meaning comes from the present of speech time since the modal operator [POSS] is anchored in the present. 'Someone else told me that I will travel/ am traveling to Dubai'.

Difference with other modal conditionals: the politeness conditional

Our last point will consist in making some differences with other conditionals. Unlike rumor conditional or conjecture conditional, politeness conditional in Spanish does not convey any epistemic meaning of possibility:

(46) Querría pedirte un favor. # 'All the Spanish people are safe".

Conclusions

In this paper we have investigated the factors that allow us to read the conditional as a rumor conditional. Our claim is that -ría is an epistemic modal operator. It operates on the illocutionary force of a situation p which has not been asserted by the Speaker himself (the rumor that p). It means "maybe the information about p is true". This situation p needs to be anchored in the present or in the future but not in the past. To make it clearer, we have compared it with the conjecture conditional which expresses present assessment of probability about a past event. In this case, -ría is an epistemic modal operator which operates on the realization of p as p refers to a fact. It means "maybe the realization of p is true". One of the remaining problems is why evidentiality only arises in Spanish with present and future situations, and never with past situations.

(

  pre-print) Axelle Vatrican, (2015) Belgian Journal of Linguistics, 29, p. 83-100 4.4. A third person Our last point deals with the use of the grammatical person. Most of the times, a rumor conditional appears with a third person, as we can observe in our examples: (38) Todos los españoles estarían a salvo, menos uno, el mallorquín Guillermo Gual del 68. ("Abandonados y enfadados", El País, 15/01/12) 'Supposedly / It is said that all the Spanish people are safe' (39) Sevilla tendría próximamente carrera de psicología para los jóvenes que deseen estudiar esta carrera en el municipio. (Canal Vivavisión, 08/07/2014) 'Supposedly / It is said that, in the future, Sevilla will offer a degree in psychology'. When using the first person I, we change the meaning of the conditional 3 as we can observe in the examples (40) and (41): (40) Vicky Xipolitakis viajaría a Dubai para ver a Maradona. (41) Yo viajaría a Dubai para ver a Maradona. With the first person yo, we obtain a hypothetical conditional adding a protasis si yo pudiese 'if I could': (42) Yo viajaría a Dubai para ver a Maradona, si yo pudiese. 'I would travel to Dubai to meet with Maradona, if I could'. We can also obtain a conjecture conditional adding a temporal complement like en 1995 'in 1995': (43) Yo viajaría a Dubai en 1995. 'I must have traveled to Dubai in 1995' 3 I thank Patrick Dendale to remind me that there exists in French a conjecture conditional in interrogative sentences: Serais-je un imbécile? This use is impossible in Spanish with a conditional. We would use a future: ¿Seré imbécil? (pre-print) Axelle Vatrican, (2015) Belgian Journal of Linguistics, 29, p. 83-100 But it is not as natural to obtain the rumor conditional meaning; if it is possible, only the future reading is possible as in (45): (44) # Según me ha informado, estoy yo viajando a Dubai. # 'Someone else told me that I am traveling to Dubai'. (45) Según me han informado, estaré yo viajando a Dubai mañana.

'

  I would like to ask you a favor.' (47) # Es probable que quiera pedirte un favor. # '{Probably/ it is possible} that I want to ask you a favor.' However an equivalence with the present tense is allowed with politeness conditional, (48) a. Querría pedirte un favor. 'I would like to ask you a favor'. b. = Quiero pedirte un favor. 'I want to ask you a favor'. but not with the rumor conditional as the contrasts between (49a) and (49b), as well as (50a) and (50b) show: (49) a. Vicky Xipolitakis viajaría a Dubai para ver a Maradona. 'It is said that V.X. is traveling to Dubai'. b. # V. X. está viajando a Dubai. # 'V.X. is traveling to Dubai'. (pre-print) Axelle Vatrican, (2015) Belgian Journal of Linguistics, 29, p. 83-100 (50) a. Todos los españoles estarían a salvo, menos uno, el mallorquín Guillermo Gual del 68. 'It is said/ It seems that all the Spaniards are safe'. b.# Todos los españoles están a salvo.

  Todos los españoles estarían a salvo, menos uno, el mallorquín Guillermo Gual del 68. ("Abandonados y enfadados", El País, 15/01/12) 'It is said/ It seems that all the Spaniards are safe. [Someone told me that… but I didn't see it']

	(6) Sevilla tendría próximamente carrera de psicología para los jóvenes que
	deseen estudiar esta carrera en el municipio. (Canal Vivavisión,
	08/07/2014)

'It seems that / it is said that, Sevilla will offer a degree in psychology.

[Someone told me that… but I didn't see it']

Spanish: a controversial form

	(pre-print) Axelle Vatrican, (2015) Belgian Journal of Linguistics, 29, p. 83-100
	"In the current century, news written this way appear with some frequency
	in the newspapers of all the Hispanic countries [...]. Obviously these are
	telegrams translated from English or French."
	Es puro francés. X rencontrerait Y prochainement. […] Pero el idioma
	cuenta, para advertir que algo no está comprobado, con propios y
	acreditados recursos: se dice…, Parece ser que…, Aseguran…, Es
	probable o posible que… (Lázaro Carreter 2001: 386)
	"It is pure French. X rencontrerait and prochainement. […] But the
	language has its accredited and own resources to advise that something is
	not checked... it is said that, it seems to be, they say, it is probable or
	possible that..."
	According to El libro de estilo de El País, it minimizes the credibility of the information:
	El uso del condicional en ese tipo de frases queda terminantemente
	prohibido en el periódico. Además de incorrecto gramaticalmente, resta
	credibilidad a la información. (El Libro de Estilo, El País, 1991: 124)
	Within the group of Romance languages, the French Rumor Conditional has been fully
	studied: it is called journalistique (Gosselin 2001), de citation (Kozen & Nølke 1990),
	d'altérité énonciative (Haillet 2002), de non prise en charge (Abouda 2001), de l'information
	incertaine (Dendale 2001). By contrast, many authors are reluctant to accept it in Spanish, at
	least until 2009, as shown in these examples. For Lázaro Carreter (2001) and Gili Gaya
	(1943), this conditional is a translation from French or English:
	En el siglo actual, aparecen con alguna frecuencia en los periódicos de
	todos los países hispánicos noticias redactadas de este modo […].
	Evidentemente se trata de telegramas traducidos del inglés o del francés.
	(Gili Gaya 1943: §129, n. 6)

  pre-print) Axelle Vatrican, (2015) Belgian Journal of Linguistics, 29, p. 83-100 'Supposedly, in the future, Sevilla will offer a degree in psychology.' (15) a.Cerca de un millón y medio de extranjeros podrían vivir con o sin papeles, en Grecia, un país de poco más de diez millones de habitantes (El

	Mundo, 04/05/2013)
	b. Cerca de un millón y medio de extranjeros supuestamente pueden
	vivir con o sin papeles.
	'Supposedly a million and a half foreigners can live with or without
	identification papers'
	(16) a.Junker aceptaría limitar la libertad de los trabajadores en la UE para
	que Londres le apoye. (El Mundo, 05/06/2014)
	b.Junker aceptará limitar supuestamente la libertad de los trabajadores
	en la UE.
	(13) a.Todos los españoles estarían a salvo, menos uno, el mallorquín
	Guillermo Gual del 68. ("Abandonados y enfadados", El País, 15/01/12)
	b. 'Todos los españoles supuestamente están a salvo'
	'Supposedly all the Spaniards are safe.'
	(14) a.Sevilla tendría próximamente carrera de psicología para los jóvenes
	que deseen estudiar esta carrera en el municipio. (Canal Vivavisión,
	08/07/2014)

1. It expresses uncertainty concerning the information conveyed; this means that a rumor conditional conveys an epistemic meaning, i.e. a uncertainty meaning as shown in the following examples with the equivalence with supposedly (supuestamente in Spanish): b. Sevilla tendrá supuestamente carrera de psicología para los jóvenes que deseen estudiar esta carrera en el municipio. ('Junker will supposedly accept freedom of workers in the EU'. (17) a.Vicky Xipolitakis viajaría a Dubai para ver a Maradona. La vedette partiría rumbo al emirato para "consolar" al recientemente separado "Diez". (La Gaceta, 8 de marzo de 2014, Argentina) b.Vicky Xipolitakis viajará supuestamente a Dubai para ver a Maradona.

There is no Conjecture Conditional neither in French nor in Italian. In French, a very restrictive case: see Dendale

2 # means 'not an equivalent'.


Benjamin Stamm

Louis Lagardère

Giovanni Scalmani

Paolo Gatto

Eric Cancès

Jean-Philip Piquemal

Yvon Maday

Benedetta Mennucci

Filippo Lipparini

| Louis Lagardère

Domain Decomposition

Keywords: Continuum Solvation, Linear Scaling, Domain Decomposition 1

How to make continuum solvation incredibly fast in a few simple steps: a practical guide to the domain decomposition paradigm for the Conductor-like Screening Model Continuum Solvation, Linear Scaling,

| INTRODUCTION

Continuum Solvation models (CSMs) [START_REF] Tomasi | Quantum Mechanical Continuum Solvation Models[END_REF][START_REF] Mennucci | Continuum Solvation Models: What Else Can We Learn from Them?[END_REF][START_REF] Cramer | Implicit Solvation Models: Equilibria, Structure, Spectra and Dynamics[END_REF][START_REF] Mennucci | Polarizable continuum model[END_REF][START_REF] Lipparini | Perspective: Polarizable continuum models for quantum-mechanical descriptions[END_REF] are nowadays part of the standard toolbox of computational chemists.

Historically, in the quantum chemistry community, polarizable CSMs such as the Polarizable Continuum Model (PCM) [START_REF] Tomasi | Quantum Mechanical Continuum Solvation Models[END_REF] or the Conductor-like Screening Model (COSMO) [START_REF] Klamt | COSMO: a new approach to dielectric screening in solvents with explicit expressions for the screening energy and its gradient[END_REF] have been developed as a cheap, but physically sound way, to include solvation effects in the quantum mechanical (QM) description of a molecule and its properties. As the computational cost is usually dominated by the solution of the QM equations, the computational performance of CSM have not been historically taken into much consideration, as the setup and solution of the CSM equations, which are equivalent in some way to solving Poisson's equation in a heterogeneous dielectric medium, has always been assumed to be a negligible additional computational cost. Due to advances in hardware, more efficient implementations and the spread of linear scaling techniques within quantum chemistry, such an assumption started to be less and less true in the last decade. Further, the diffusion of multiscale methods such as quantum mechanics/molecular mechanics (QM/MM) has made large to very large systems accessible to computational chemists. For such systems, the computational cost associated with continuum solvation, which scales as the second or third power of the size of the system, can easily become the real bottleneck of the calculation [START_REF] Lipparini | Quantum calculations in solution for large to very large molecules: A new linear scaling QM/continuum approach[END_REF]. Nevertheless, polarizable CSM have been used successfully in a large number of different applications and are a standard feature of most quantum chemistry codes. Polarizable CSM have not had the same success in Molecular Dynamics (MD) applications, as their computational cost makes them incompatible with the usual time frame of these simulations. For this reason, the CSMs used in classical MD are usually additive approximation to Poisson's equation such as the Generalized Born model [START_REF] Still | Semianalytical treatment of solvation for molecular mechanics and dynamics[END_REF].

Recently, a new paradigm has been introduced for the solution of the COSMO equations [START_REF] Cancès | Domain decomposition for implicit solvation models[END_REF][START_REF] Lipparini | Fast domain decomposition algorithm for continuum solvation models: Energy and first derivatives[END_REF]. The new algorithm, called ddCOSMO, is based on domain decomposition (dd) and exhibits linear scaling and an overall very limited computational cost [START_REF] Lipparini | Fast domain decomposition algorithm for continuum solvation models: Energy and first derivatives[END_REF][START_REF] Lipparini | Quantum calculations in solution for large to very large molecules: A new linear scaling QM/continuum approach[END_REF][START_REF] Lipparini | Quantum, classical, and hybrid QM/MM calculations in solution: General implementation of the ddCOSMO linear scaling strategy[END_REF]. With respect to existing implementations of CSM [START_REF] Scalmani | Achieving linear-scaling computational cost for the polarizable continuum model of solvation[END_REF][START_REF] Lipparini | Quantum calculations in solution for large to very large molecules: A new linear scaling QM/continuum approach[END_REF], ddCOSMO is two to three orders of magnitude faster, effectively allowing the use of a polarizable CSM for large systems. Furthermore, the ddCOSMO discretization is systematically improvable and controlled by a very limited number of parameters, making it easy to control and overall rigorous and sound [START_REF] Cancès | Domain decomposition for implicit solvation models[END_REF]. It also provides for a smooth energy as a function of the nuclear coordinates, making it suitable for geometry optimization and molecular dynamics simulation [START_REF] Lipparini | Polarizable molecular dynamics in a polarizable continuum solvent[END_REF][START_REF] Lipparini | Quantum calculations in solution for large to very large molecules: A new linear scaling QM/continuum approach[END_REF]. However, its implementation is slightly more cumbersome than standard CSM implementations based on the boundary element method (BEM) [START_REF] Cammi | Remarks on the use of the apparent surface-charges (ASC) methods in solvation problems -Iterative versus matrix-inversion procedures and the renormalization of the apparent charges[END_REF] or the York-Karplus (YK) method [START_REF] York | A smooth solvation potential based on the conductor-like screening model[END_REF][START_REF] Scalmani | Continuous surface charge polarizable continuum models of solvation. I. General formalism[END_REF][START_REF] Lange | A smooth, nonsingular, and faithful discretization scheme for polarizable continuum models: The switching/Gaussian approach[END_REF]. A stand-alone, opensource, modular implementation of ddCOSMO is distributed for free under the terms of the version 3 of the GNU Lesser General Public Licence (LGPL) license. It can be download from https://www.ddpcm.org, or directly from GitHub [START_REF] Lipparini | A fast domain decomposition based implementation of the COSMO solvation model[END_REF].

The goal of this communication is to provide the essential knowledge to successfully use and couple the stand-alone library to the code of choice, or implement the method from scratch. We provide all necessary technical details and review the implementation of ddCOSMO and its coupling with an existing QM or classical code in this communication.

This paper is organized as follows. In section 2, the main aspects of the ddCOSMO paradigm are succinctly presented. In section 3, we will discuss the coupling of ddCOSMO with an existing code and review what are the quantities that the existing code needs to assemble in order to compute the self-consistent field energy for a QM code and the energy and forces for a classical code. In section 4, the implementation of ddCOSMO is described with explicit references to the open source ddCOSMO module. Finally, we will provide some conclusion and perspectives in section 5.

| A BRIEF INTRODUCTION TO THE ddCOSMO PARADIGM

The idea of this section is to give a different derivation of the ddCOSMO equations than what was presented in the original articles [START_REF] Cancès | Domain decomposition for implicit solvation models[END_REF][START_REF] Lipparini | Fast domain decomposition algorithm for continuum solvation models: Energy and first derivatives[END_REF][START_REF] Lipparini | Quantum, classical, and hybrid QM/MM calculations in solution: General implementation of the ddCOSMO linear scaling strategy[END_REF][START_REF] Stamm | A new discretization for the polarizable continuum model within the domain decomposition paradigm[END_REF] that might be more intuitive to follow. In particular, we refer to the appendix B of [START_REF] Gatto | Computation of Forces arising from the Polarizable Continuum Model within the Domain-Decomposition Paradigm[END_REF] which presents an alternative derivation of the ddCOSMO equations. We start with recalling the Partial Differential Equation (PDE) associated to the COSMO (model). The reaction potential W satisfies

-∆W = 0, in Ω = Ω 1 ∪ • • • ∪ Ω M , W = -Φ, on ∂ Ω,
where ∆ denotes the Laplacian, Φ the free-space potential generated by the solute's charge distribution ρ and Ω is the solute's cavity consisting of a union of, possibly scaled, Van der Waals balls Ω i , each centered at x i and of radius r i .

The idea of ddCOSMO is to rewrite this problem into M coupled problems:

-∆W i = 0, in Ω i , (1)
W i = -Φ, on ∂ Ω i ∩ ∂ Ω, (2)
W i = j ∈N i ω i j W j , on ∂ Ω i ∩ Ω, (3)
where the right hand side of (3) at some point x ∈ ∂ Ω i denotes the average value of all neighboring potentials W j (x )

corresponding to all neighboring spheres Ω j intersecting Ω i at x , i.e. for all j i such that x ∈ Ω j . More precisely, we denote by N i the set of all spheres intersecting sphere i and the function ω i j denotes the weight function defined by

x ∈ ∂ Ω i : ω i j (x ) = χ j (x ) k ∈N i χ k (x )
, using the convention that 0/0 = 0 and where χ j is the characteristics function of Ω j . Therefore, ω i j (x ) can take the values 0 (exposed to the solvent at x ), 1 (one intersecting sphere at x ), 1/2 (two intersecting spheres at x ), 1/3, 1/4 etc.

Following this construction, we deduce that U i (x ) = 1j ∈N i ω i j (x ) denotes the characteristics function of the part of

∂ Ω i belonging to ∂ Ω taking only values 0 (inside Ω) or 1 (on ∂ Ω).
On the discrete level, we represent the local approximations W i to W in each Ω i by

x ∈ Ω i : W | Ω i (x ) ≈ W i (x ) = m 4π 2 + 1 [X i ] m r i (x )Y i ,m (x ), (4)
with radial scaling r i (x ,m up to degree max and integrating over ∂ Ω i . However, since exact integration cannot be carried out, we employed a Lebedev quadrature [START_REF] Lebedev | A quadrature formula for te sphere of the 131st algebraic order of accuracy[END_REF] given by weights ω n ∈ and integration points s n ∈ 2 on the unit sphere to provide a (scaled) numerical integration scheme on each ∂ Ω i :

) := |x -x i | r i and angular dependencyY i ,m (x ) := Y m x -x i |x -x i | relative
f , g n,i := Ng n=1 ω n f (x i n ) g (x i n ) ≈ 1 r 2 i ∫ ∂ Ω i f (s)g (s)ds = ∫ 2 f (x i + r i ŝ)g (x i + r i ŝ)d ŝ,
with x i n := x i + r i s n . Then, we require the following equations to hold for all i = 1, . . . , M , = 0, . . . , max and |m | ≤ :

W i ,Y i ,m n,i = j ∈N i ω i j W j -U i Φ,Y i ,m n,i . (5)
This results in a linear system

LX = g (6)
which is naturally block-sparse since only intersecting spheres provide non-zero block entries. While X is just a vector containing the coefficients [X i ] m , ordered block-wise for all spheres i , the matrix and right hand side coefficients are given by (for all i = 1, . . . , M , j ∈ N i ):

[L i i ] mm = 4π 2 +1 Y i ,m ,Y i ,m n,i = 4π 2 +1 δ δ mm , (7)
[L i j ] mm = -4π 2 +1 ω i j r j Y j ,m ,Y i ,m n,i = -4π 2 +1 Ng n=1 ω n ω i j (x i n ) r j (x i n )Y j ,m (x i n )Y m (s n ), (8)
[g i ] m = -U i Φ,Y i ,m n,i = - Ng n=1 ω n U i (x i n )Φ(x i n )Y m (s n ). (9)
If this linear system is solved by means of any numerical method for non-symmetric matrices such as the generalized minimal residual (GMRes) method or coupled Jacobi/direct inversion in the iterative subspace [START_REF] Pulay | Convergence acceleration of iterative sequences. the case of scf iteration[END_REF] (DIIS) iterations [START_REF] Lipparini | Scalable Evaluation of Polarization Energy and Associated Forces in Polarizable Molecular Dynamics: I. Toward Massively Parallel Direct Space Computations[END_REF],

the energy can then be computed following

E s = 1 2 f (ε s ) ∫ Ω ρ(x )W (x )dx or, its discrete variant E N s = 1 2 f (ε s ) i m [X i ] m [Ψ i ] m = 1 2 f (ε s ) X , Ψ ,
where we introduced the compact notation

A, B = i m [A i ] m [B i ] m
and where the definition of the vector Ψ with coefficients [Ψ i ] m depends on the nature of the solute's charge, i.e. whether the COSMO is coupled with a classical, quantum mechanical or hybrid model, which will be specified in the next section.

Please also note that in order to achieve a smooth energy while varying nuclear coordinates, the weight functions ω i j and U i are replaced by smooth approximate functions in practice. To do so, let us first introduce the polynomial of degree 5 given by

p η (t ) = η -5 (1 -t ) 3 (6t 2 + (15 -12)t + 1 -η 2 -15η + 6),
so that the regularized characteristic function defined by

χ η (t ) =              1 if t ≤ 1 -η, p η (t ) if 1 -η < t < 1, 0 if t ≥ 1,
is twice continuously differentiable, in order to ensure the continuity of the forces and second-derivatives of the energy.

Finally, the regularized versions of ω i j and U i are provided by the following expressions

ω i j η (x ) = d i (x i n ) χ η r j 1 (x ) resp. U i η (x ) = 1 - j ∈N i ω i j η (x ) with d i (x ) = min k ∈N i χ η r k 1 (x ) , 1 k ∈N i χ η (r k 1 (x ))
.

The construction of this regularization is a delicate task and its explanation goes beyond the scope of this communication [START_REF] Lipparini | Fast domain decomposition algorithm for continuum solvation models: Energy and first derivatives[END_REF].

In order to compute the ddCOSMO forces or other derived quantities, analytical derivatives of the ddCOSMO energy are required. Differentiating eq .2 with respect to a generic parameter λ we get

dE s dλ = E λ s = 1 2 f (ε s ) X λ , Ψ + X , Ψ λ , ( 10
)
where we denote the derivative with respect to λ with a superscript. Note that the function Ψ depends explicitly on λ, while the coefficients [X ] depend implicitly on λ. In order to compute the derivatives of the coefficients, we differentiate the ddCOSMO equation

(LX ) λ = L λ X + LX λ = g λ ⇒ X λ = L -1 (g λ -L λ X ). (11)
By substituting eq. 11 into eq. 10, we get

E λ s = 1 2 f (ε s ) L -1 (g λ -L λ X ), Ψ + X , Ψ λ = 1 2 f (ε s ) g λ -L λ X , (L -1 ) † Ψ + X , Ψ λ . (12)
If we now introduce the ddCOSMO adjoint equations

L † S = Ψ, (13)
the general expression for the analytical derivative of the ddCOSMO energy becomes

E λ s = 1 2 f (ε s ) g λ -L λ X , S + X , Ψ λ . (14)
In conclusion, in order to realize the ddCOSMO algorithm, one needs to be able to compute the right-hand side (g , eq. 9) to the ddCOSMO equation (eq. 6) and to solve the linear system. Then, to compute the energy, one has to assemble the function Ψ (eq. 2). Any further computation requires to solve the ddCOSMO adjoint equation (eq. 13) and to compute the derivative of the ddCOSMO energy (eq. 14), which in turn requires the computation of the derivatives of g and Ψ, and the product between the coefficients X and the differentiated ddCOSMO matrix. In section 3, all these terms will be detailed for QM and classical solutes.

| COUPLING ddCOSMO TO AN EXISTING CODE

In this section, we discuss how to implement ddCOSMO in an existing code (called the host code in the following) by coupling such a code to the ddCOSMO module 1 . We recall that the ddCOSMO code can be obtained for free [START_REF] Lipparini | A fast domain decomposition based implementation of the COSMO solvation model[END_REF] under the terms of the LGPL license (version 3), which is a weak copyleft license that allows developers and companies to use, modify, redistribute and integrate the code with the only constraint that the modified code has to be made available under the terms of the same license. If a proprietary code includes ddCOSMO, only the modified ddCOSMO library has to be made available with its source code, thus protecting the intellectual property of both the authors of the library and of the company that produces the proprietary code.

First, we give an overview of the main components of a modular implementation of ddCOSMO. Then, we analyze more in the specific the quantities that the host code needs to provide to the ddCOSMO module and how the host code has to digest the output of ddCOSMO in order to proceed further. We consider two different cases. First, we analyze the calculation of the (self-consistent) energy for a solute described at a QM level of theory such as Hartree-Fock or Kohn-Sham Density Functional Theory, namely, a method that requires the solution of the self-consistent field (SCF)

equations. Second, we consider the case of a classical force field energy and forces calculation, such as the one required to perform a step of molecular dynamics.

| The general picture

The most modular implementation possible would almost not require any new coding and be based on the input/output of the host code and of the ddCOSMO module, with the only coding required being limited to parsing and writing input and output files. While this is in principle possible, it would be terribly inefficient, especially when the solute is described with a QM method, or even a polarizable force field. We will therefore adopt the following philosophy. We divide the code in three sets of routines

• Host software routines. These routines are independent of ddCOSMO and are assumed to be already present in the host code . They perform general operations such as the evaluation of electrostatic quantities at provided grid points, the numerical quadrature for density functional theory or one and two electron integrals. No modification to such routines is needed in order to achieve the coupling.

• ddCOSMO routines. These routines are independent of the host code and deal exclusively with ddCOSMO-related quantities. They include setup routines, the routines needed to solve the ddCOSMO linear equations and routines that compute differentiated ddCOSMO quantities. No modification to such routines is needed in order to achieve the couplings.

• driver routines. These routines take care of the actual coupling by calling the appropriate host software or ddCOSMO routines. The implementation of driver routines is therefore the main task to be performed in order to achieve the coupling between ddCOSMO and the host software .

In this contribution, we only describe the functionalities that the host software routines need to provide. As these routines are strongly code-dependent we cannot provide specific details. Nevertheless, the operations to be performed are described in detail from an abstract point of view. In the following, we will discuss how to set up a driver. The ddCOSMO routines will be the subject of section 4.

| Initializing ddCOSMO

Independent of the method and level of theory to which ddCOSMO is coupled, any computation involving ddCOSMO requires some geometrical input data comprising: • the number of atoms forming the molecule (n);

TA B L E
• the position of the centers of the spheres that make the cavity (x, y, z, arrays of size n);

• the radii of the spheres that make the cavity (rvdw, array of size n) (for instance, the UFF radii scaled by 1.1);

as well as a few scalar parameters

• the order max of the spherical harmonics expansion used in the discretization (lmax) ;

• the number N g of Lebedev points to be used for the numerical integration (ngrid);

• the width η of the switching region relative to the radii of the spheres (eta);

• the solvent's dielectric constant ε (eps);

• the convergence threshold (iconv, the threshold is defined as 10 -iconv );

• a flag that, if set to a nonzero value, indicates that the computation of the forces is required (igrad).

The variable names in our implementation of these global variables is emphasized in the parentheses. Such variables need to be initialized by the user through the host code . With these input, an internal data-structure is built that comprises

• the list N i of intersecting spheres for each sphere (nl, see section 4.1 for further detail);

• Lebedev integration points s n and weights ω n (grid and w, arrays of dimension (3,ngrid); and (ngrid), respectively);

• the number and coordinates of the points that lie on a portion of the cavity which is exposed to the solvent (ncav and ccav, array of dimension (3,ncav));

• several quantities that are precomputed in order to speed up the computation (see section 4.1 ).

In a Fortran 90 implementation, these data as well as the above presented geometrical data and scalar parameters can be accessed when using our module by including use ddcosmo in the host code . The routine that handles the computations of this data-structure is called DDINIT and needs to be invoked before any ddCOSMO computation. It can be invoked as follows:

>> call ddinit(n,x,y,z,rvdw).

In table 1, a suggested set of default parameters is provided for both QM and classical applications. These parameters are, in our experience, safe combinations that result in a good compromise between the numerical accuracy and efficiency of ddCOSMO.

| Case 1: SCF energy with ddCOSMO

The coupling of ddCOSMO to the SCF procedure requires, at each SCF iteration, not only to compute the solvation energy and add it to the total energy, but also to add a contribution to the Fock (or Kohn-Sham) matrix. Both these contributions can be assembled right after the two-electron contribution to the Fock matrix has been computed. The main difficulty introduced by ddCOSMO in a QM computation [START_REF] Lipparini | Quantum, classical, and hybrid QM/MM calculations in solution: General implementation of the ddCOSMO linear scaling strategy[END_REF] lies in the way the ddCOSMO solvation energy is computed. Contrary to the usual discretizations of COSMO, ddCOSMO does not give access to an apparent surface charge, but relies entirely on the reaction potential W , which can be computed inside each sphere locally starting from the solution to the ddCOSMO equations X .

First, one has to assemble the right-hand side for the ddCOSMO linear system, defined in eq. 9. Note that the potential Φ is defined, for a QM solute, as

Φ(x i n ) = K Z K |x K -x i n | - µν P µν ∫ 3 χ µ (x )χ ν (x ) |x -x i n | dx, (15)
where we recall that x n i = x i + r i s n , Z K is the charge of the K -th nucleus and x K its position, and P µν is the density matrix in the basis of the atomic orbitals {χ µ }. Note that the value of the potential at the internal points is not needed, as the characteristic function of the external surface appears explicitly in the ddCOSMO right-hand side (eq. 9). The solvation energy is defined as

E s = 1 2 f (ε) ∫ Ω ρ(x )W (x )dx, (16)
where for a QM solute the density is given by

ρ(x ) = K Z K δ(x -x K ) - µν P µν χ µ (x )χ ν (x ) = ρ nuc + ρ ele . ( 17
)
As seen in eq. 4, the reaction potential can be expressed inside a given sphere Ω j as

W | Ω j (x ) ≈ W j (x ) = m 4π 2 + 1 [X j ] m r j (x )Y i ,m (x ).
In order to take advantage of the ddCOSMO local representation of W , eq. 16 can be rewritten as a sum of contributions over the spheres. However, in order to avoid double counting contributions arising from the intersection of spheres, a weight function θ j has to be introduced:

E s = 1 2 f (ε) j ∫ Ω j ρ(x )θ j (x )W j (x )dx . ( 18
)
The integral in eq. 18 can be evaluated numerically using the same quadrature that is used for the exchange-correlation (XC) potential in density functional theory (DFT). Such a quadrature [START_REF] Becke | A multicenter numerical integration scheme for polyatomic molecules[END_REF] is in fact performed with a set of atom-centered points and uses weights to avoid double counting. Let {z j g , θ j g } be the set of points and weights used for the DFT quadrature around the j -th atom. We can then write, for the contribution stemming from the electronic density [START_REF] Lipparini | Quantum, classical, and hybrid QM/MM calculations in solution: General implementation of the ddCOSMO linear scaling strategy[END_REF]:

E s = 1 2 f (ε) j g θ j g ρ ele (z j g )W j (z j g ) = 1 2 f (ε) j g θ j g ρ(z j g ) m 4π 2 + 1 r j (z j g )Y j ,m (z j g )[X j ] m . ( 19
)
If we compare eq. 19 with eq. 2, we get

[Ψ ele j ] m = g θ j g ρ ele (z j g ) 4π 2 + 1 r j (z j g )Y j ,m (z j g ). (20)
The vector Ψ ele can therefore be computed by using the same code that is used for XC quadrature in DFT. Note that DFT grids are usually built as the product of a radial and an angular grid. As the quantities evaluated at grid points in eq. 20 can also be clearly separated into a radial and an angular part, it is possible to precompute radial and angular intermediates and to assemble them while evaluating the integral. This makes the additional ddCOSMO cost negligible with respect to the cost of assembling the grids and evaluating the density at the grid points, two operations that are anyway needed in a DFT calculation. The nuclear contribution can be computed separately from the electronic one and added to form the total Ψ. For the nuclear charge, we get

[Ψ nuc j ] m = 2 √ π Z j δ 0 δ m0 (21)
A numerical quadrature is also necessary in order to assemble the ddCOSMO contribution to the Fock matrix, that is defined as the derivative of the ddCOSMO energy with respect to the density matrix

F s µν = dE s dP µν = 1 2 f (ε s ) ∂g ∂P µν , S + X , ∂ Ψ ∂P µν , (22)
where we used eq. 14 and the fact that the ddCOSMO matrix does not depend on the density. The first contribution can be rearranged as follows:

∂g ∂P µν , S = - j n m [S j ] m ω n U j η (x j n )Y m (s n ) ∂ Φ(x j n ) ∂P µν . ( 23
)
The sum over and m in eq. 23 can be computed before the evaluation of the integrals. Introducing the new quantities

ξ j n = -U j η (x j n )ω n m [S j ] m Y m (s n ), (24)
eq. 23 becomes

∂g ∂P µν , S = - j n ξ j n ∫ 3 χ µ (x )χ ν (x ) |x j n -x | dx . (25)
Note that ξ n i vanishes at points buried inside the cavity, so that it only needs to be evaluated and contracted with the integrals at the external points. The second contribution to the Fock matrix requires a second numerical integration. We get

X , ∂ Ψ ∂P µν = - j g m [X j ] m θ j g χ µ (z j g )χ ν (z j g ) 4π 2 + 1 r j (z j g )Y j ,m (z j g ). (26)
Again, the sum over and m can be computed before evaluating the integral. By introducing

η j g = m 4π 2 + 1 r j (z j g )Y j ,m (z j g )[X j ] m (27)
we get

X , ∂ Ψ ∂P µν = - j g
η j g θ j g χ µ (z j g )χ ν (z j g ).

We have now defined all the quantities required for a SCF/ddCOSMO calculation and can describe the structure of the driver routine. We focus on a single SCF iteration and underline the quantities that can be precomputed at the beginning. We also provide details on the specific calls to the ddCOSMO routines included in the ddCOSMO module, commenting on the interface.

1. Initial setup, as described in section 3.2.

2.

Computation of the solute's potential at the cavity points. The driver needs to call the appropriate, host code specific, integral routine to compute the potential (see eq. 15) at the external points (in ccav).

Φ i = K Z K |x K -x i n | - µν P µν ∫ 3 χ µ (x )χ ν (x ) |x -x i n | dx, (29)
where the potential needs to be computed only at the points x n i that lie on the external cavity surface (i.e., the ones contained in the ccav array). We assume that the potential is stored in an array phi of size ncav.

Computation of Ψ:

The driver needs to call the quadrature routines to compute the array psi, an array of size ((lmax+ 1) 2 ,n). The details of how this array is assembled are strictly dependent on the host code , but eq. 20 should be straightforward to implement given the DFT quadrature routines.

4.

Solution to the ddCOSMO equation and computation of the ddCOSMO energy. This operation is easily performed by calling the COSMO routine from the ddCOSMO module. The driver needs to allocate space for the X coefficients, which are stored in an array xcoef of dimension ((lmax+1) 2 , n). The COSMO routine can accept as an input either the ddCOSMO right-hand side (see eq. 9) or the potential at the external points, which is then processed internally.

Assuming the latter case, the routine should be called as follows:

>> call cosmo(.false.,.true.,phi, xx, psi, xcoef, esolv).

The first argument specifies whether the ddCOSMO equation (.false.) or its adjoint (.true.) should be solved.

xx is a dummy array that acts as a placeholder for the ddCOSMO right-hand side and is not referenced when the second argument is set to .true.. esolv is a real scalar that contains in output the ddCOSMO electrostatic solvation energy E s .

5.

Solution to the ddCOSMO adjoint equation. Again, this is achieved by calling the COSMO routine from the ddCOSMO module. The driver needs to allocate space for the S coefficients, which are stored in an array scoef of dimension ((lmax+1) 2 , n). With respect to the previous call, the only differences are that the first logical has to be set to .true.and that the output has to be saved in the array scoef:

>> call cosmo(.true.,.true.,phi, xx, psi, scoef, esolv).

6.

Computation of the ddCOSMO contributions to the Fock matrix. The driver needs here to call the proper integral and numerical quadrature routines from the host code in order to assemble the contributions in eqs. 23 and 26.

For improved efficiency, the intermediates ξ and η should be precomputed and used as in eqs. 25 and 28. The ξ intermediate can be computed by invoking the DDMKXI routine from the ddCOSMO module, which is invoked with the following arguments >> call ddmkxi(scoef, xi).

Note that an array xi of size ncav should be allocated by the user. The xi array can then be passed to the integrals evaluation routine. There is no routine to assemble the η intermediate in the ddCOSMO module, as it is assumed that the DFT quadrature points are not easily available in the driver. η can be evaluated on-the-fly inside the integration routine whenever a batch of integration points is processed. This can be done before the loop over the basis functions for improved efficiency.

| Case 2: Energy and forces with a MM force field

Interfacing ddCOSMO to a classical code is overall much simpler than interfacing it to a quantum mechanical code.

In this section, we discuss how to compute the ddCOSMO contributions to the energy and forces when the solute is described with a classical force field. We assume that the electrostatic interactions are modeled by endowing each atom with a point charge, as is commonly done in popular force fields such as AMBER [START_REF] Cornell | A 2nd Generation Force-Field for the Simulation of Proteins, Nucleic-Acids, and Organic-Molecules[END_REF], CHARMM [START_REF] Brooks | CHARMM: A program for macromolecular energy, minimization, and dynamics calculations[END_REF], OPLS [START_REF] Jorgensen | Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids[END_REF] and others. The ddCOSMO contribution to the forces encompasses both ddCOSMO specific quantities and quantities that depend on the solute and are to be computed by a routine from the host code . Differentiating the ddCOSMO energy with respect to the nuclear coordinates of the i -th atom, we get

i E s = S, i g -( i L)X , (30)
as Ψ does not depends on the position of the nuclei when the solute's density of charge is a collection of point charges, as it is apparent from eq. 21. Here, i denotes the derivative with respect to the nuclear coordinate x i . The terms stemming from the derivatives of the L matrix are ddCOSMO exclusive quantities. For the derivatives of the right-hand side, we get

i [g j ] m = - Ng n=1 ω n Y m (s n ) i U j η (x j n )Φ(x j n ) + U j η (x j n ) i Φ(x j n ) . (31)
The first contribution is again a ddCOSMO exclusive quantity and originates from the fact that the characteristic function U j is regularized by U j η , which is required to have a continuous energy as a function of the nuclear coordinates. The second term requires to compute the derivatives of the solute's potential and requires a host code routine that computes electrostatic quantities.

Let us now describe the structure of a MM/ddCOSMO driver.

1. Initial setup. This step is identical to the one described in section 3.2 as already mentioned in section 3.3 and is readily achieved by setting the ddCOSMO scalars and invoking the DDINIT routine from the ddCOSMO module.

2.

Compute the solute's potential at the external cavity points. For a set of point charges q K , the electrostatic potential is given by

Φ(x i n ) = Na K =1 q K |x K -x i n | . ( 32
)
For a simple-minded implementation, a naive double-loop routine can be easily written in a few lines of code. An example of such a naive implementation is provided in the ddCOSMO module (subroutine MKRHS). However, for molecular systems as large as the ones that are usually treated with classical molecular mechanics (MM), a fast summation technique is required in order to achieve linear scaling and speed up the computation. The most obvious choice is probably to resort to the Fast Multiple Method [START_REF] Greengard | A fast algorithm for particle simulations[END_REF] (FMM), for which many open source and highly optimized implementations exist. Nevertheless, the result of this step should be, as in section 3.3, an array phi of size ncav containing the potential at the external points. 3. Compute Ψ. This is trivially done for a classical solute, as only eq. 21 is involved. We get

[Ψ j ] m = 2 √ πq j δ 0 δ m0 . (33)
The MKRHS subroutine in the ddCOSMO module also computes Ψ.

4.

Solve the ddCOSMO equations. This step is independent of the level of theory used for the solute and can be performed as described in section 3.3, by invoking >> call cosmo(.false.,.true.,phi, xx, psi, xcoef, esolv).

Again, the ddCOSMO contribution to the solvation energy (in atomic units) is given as an output in esolv.

5.

Solution to the ddCOSMO adjoint equation. This step is needed if the computation of the forces is required and can be achieved by calling the COSMO routine >> call cosmo(.true.,.true.,phi, xx, psi, scoef, esolv).

6.

Computation of the forces -ddCOSMO exclusive terms. These contributions stem from the derivatives of the ddCOSMO L matrix (see eq. 14) and from the derivatives of the characteristic function (see eq. 31.

). The FORCES_DD subroutine from the ddCOSMO package takes care of assembling these terms. The user should allocate a temporary array for the ddCOSMO forces fx of dimension (3,n) and invoke the FORCES_DD subroutine as follows >> call forces_dd(n,phi,xcoef,scoef,fx).

7.

Computation of the forces -code-specific terms. The remaining contribution to the forces involves the derivatives of the potential. For the i -th atom, this terms contributes as follows (note that the forces are minus the derivative of the energy):

F non dd i = j m [S j ] m n w n U j η (x j n )Y m (s n ) i Φ(x j n ). (34)
Let

ζ j n = w n U j η (x j n ) m [S j ] m Y m (s n ). ( 35
)
Through simple algebra, we get

F non dd i = - n ζ i n k q k (x i n -R k ) |x i n -R k | 3 -q i j n ζ j n (R i -x j n ) |R i -x j n | 3 . ( 36
)
The first term can be interpreted as the product of (minus) the electric field produced by the solute's charges at the external cavity points interacting with the fictitious charge distribution ζ. The second, on the contrary, can be seen as the field produced by the distribution ζ at the atoms interacting with the MM charges. While the physical interpretation of these two terms is, of course, inessential, it makes the use of standard electrostatics machinery possible, including taking advantage of linear scaling techniques such as the FMM. A simple-minded routine (subroutine EFLD) that computes the electric field of a source at a set of target points with a double loop is included in the ddCOSMO package. Furthermore, a subroutine that, given the S coefficients produces the array ζ (subroutine DDMKZETA) is provided. Note that, due to the presence of the external characteristic function in the definition of ζ, such a quantity is only nonzero at the external cavity points, that are given in the ccav array. As a consequence, the two electric fields computations required to evaluate eq. 36 can be restricted to targets (sources) lying just on such points. Therefore, in order to compute this contribution, the user needs to allocate an array zeta of dimension ncav and fill it via >> call ddmkzeta(scoef,zeta).

The appropriate, host code specific routine to compute electric fields are then needed in order to evaluate the "fields" in eq. 36. Finally, eq. 36 can be evaluated.

This concludes the presentation of the MM/ddCOSMO driver. The ddCOSMO package includes a prototypical implementation of a MM driver (program main, in main.f90), which reads the MM charges, coordinates and Van der Waals radii from a text file, initializes the ddCOSMO environment and proceeds with the steps described in this section.

The implementation is extensively commented and can be used as a starting point for interfacing ddCOSMO with a classical MM code.

| THE ddCOSMO MODULE

In this section, a more detailed description of the implementation of ddCOSMO is provided. Besides providing a complete description of the existing open source code, this section provides guidelines for a new implementation, which can be used in alternative to the existing one for better integration in an existing code, in particular when such a code is written in a language different from Fortran 90.

There are three main operations required in order to perform a ddCOSMO computation. Such operations are mentioned as the ddCOSMO specific steps in section 3 and can be summarized as follows

initial setup;

2. solution to the ddCOSMO linear systems;

3. computation of the ddCOSMO specific contributions to the forces.

In the following, each of these steps is detailed.

| Initial setup of a ddCOSMO calculation

The setup of a ddCOSMO calculation requires the computation of several quantities that depend on the geometry of the solute and on the cavity (i.e., on the van der Waals radii used). Some of the quantities are precomputed for the sake of efficiency but can be also computed on-the-fly whenever needed. Given the limited amount of memory required to store such quantities, we believe that it is always worth precomputing them, as they allow for non-negligible savings in computation time.

As mentioned in section 3, the first operation that needs to be performed is to set the ddCOSMO scalars. We also set the number of spheres (nsph) to the number of atoms of the solute and compute nylm=(lmax+1)^2. The first quantity that needs to be assembled is the integration grid and weights. This requires to allocate two arrays, one for the grid (grid) of size (3,ngrid), and one for the weights (w) of size (ngrid). In our implementation, we use the publicly available routines by Laikov and van Wuellen (see http://www.ccl.net/cca/software/SOURCES/FORTRAN/Lebedev-Laikov-Grids).

We then create two arrays, one for the centers of the sphere (csph) of size (3,nsph), that coincide with the positions of the atoms, and one for their radii (rsph) of size (nsph).

A neighbor list, defined as the list of spheres that overlap each sphere, is then computed. The neighbor list is of crucial importance, as only overlapping spheres interact or, in more specific terms, only the off-diagonal blocks of the ddCOSMO matrix that correspond to two overlapping spheres are nonzero. Therefore, the neighbor list is also used to describe the sparsity pattern of the ddCOSMO matrix. Neighbor lists can be computed using several algorithms. In our implementation, we use a simple-minded double loop algorithm. This is of course not an efficient choice, which was made for the sake of simplicity. Routines that assemble neighbor lists with a computational cost that scales linearly with respect to the number of atoms are commonly available in any MD code and open source implementations can be easily found. We store the neighbor list using two integer arrays. The first, inl, of dimension nsph+1, contains a list of integers that indicate where the list of neighbors of the i -th sphere start. The second array, nl, contains the actual compressed neighbor list. According to this convention, the neighbors of the i -th sphere are nl(inl(isph)), . . . , nl(inl(isph+1))-1.

For the sake of convenience, we also precompute several quantities that are used repeatedly during a ddCOSMO calculation. We allocate two real arrays, facs and facl, of dimension nylm, that contain the normalization factors for the spherical harmonics functions (N M , in facs)

N m = (-1) m √ 2 2 + 1 4π ( -|m |)! ( + |m |)!
and the inverse diagonal of the diagonal blocks of the ddCOSMO L matrix (in all blocks are equal)

[L i i ] mm = 2 + 1 4π δ δ mm .
We then compute the set of spherical harmonics {Y m (s n )}, = 0, . . . , max , m = -, . . . , at the grid points s n for n = 1, . . . , N g , and save it in an array basis of dimension (nylm,ngrid). The computation of spherical harmonics at a given point on the unit sphere (performed by the subroutine ylmbas) is a key component of the ddCOSMO algorithm and a crucial factor in the overall efficiency of the code. The algorithm that we use is described in ref. [START_REF] Lipparini | Fast domain decomposition algorithm for continuum solvation models: Energy and first derivatives[END_REF] and is based on a recursion formula for the generalized Legendre polynomials (subroutine polleg) and on Chebyshev polynomials to evaluate the trigonometric functions (subroutine trgev).

Finally, we precompute two geometric quantities, that are defined at each grid point for each sphere and saved into two arrays, ui and fi, of dimension (ngrid,nsph). These arrays are used in order to evaluate the external and internal characteristic functions and are defined as follows

f i n = j ∈N i χ η (r j 1 (x i n )) and u i n = max(1 -f i n , 0) = U i η (x i n ). (37)
If the forces are to be computed, it is also convenient to assemble

z i n = j ∈N i 1 r j χ η (r j 1 (x i n ))e j (x i n ), (38)
where we introduced the unit vector e i (x ) =

x -x i |x -x i | and χ η denotes the derivative of the one-dimensional switching function χ η . The quantity z i n is stored in an array zi of dimension (3,ngrid,nsph).

| ddCOSMO linear system

Once the setup has been done, and the solute's potential and ψ have been computed, one has to solve the ddCOSMO linear system. In our implementation, this is done by the cosmo routine, that has been described in section 3. In this section, we describe such a routine in detail. The first task performed by the cosmo routine is to prepare the right-hand side, as in eq. 9. The solute's potential is only computed at the external cavity points, defined as the points for which u i n is nonzero. Therefore, the potential is first decompressed from ncav to (ngrid,nsph) and multiplied for the external characteristic function. Then, the numerical integration is carried out and the right-hand side assembled.

The ddCOSMO sparse linear system can be solved using any iterative solver capable of handling non symmetric linear systems, such as GMRES or the biconjugate gradient method (BCG). The use of a Jacobi solver [START_REF] Lipparini | Scalable Evaluation of Polarization Energy and Associated Forces in Polarizable Molecular Dynamics: I. Toward Massively Parallel Direct Space Computations[END_REF] coupled with Pulay's direct inversion in the iterative subspace [START_REF] Pulay | Convergence acceleration of iterative sequences. the case of scf iteration[END_REF] (DIIS) is a viable alternative to these standard solvers. Such a strategy has proven to be optimal and is the one that we follow in the ddCOSMO module. A Fortran 90 implementation of the Jacobi/DIIS solver is provided in the ddCOSMO module. Independent of the specific solver, we assume that its implementation requires as an input the name of a subroutine that performs a matrix-vector product y = Lx and possibly a routine that applies a preconditioner. Both routines are provided with a standard interface that requires as argument the length of the x and y vectors and the vectors themselves, i.e., the subroutine is defined as subroutine matvec(n,x,y). For the Jacobi solver, only the off-diagonal part of the L matrix is relevant for the matrix-vector product routine, while the preconditioner applies the inverse diagonal of the L matrix to a vector and has the same interface than the matrix-vector multiplication routine.

The matrix-vector multiplication routine performs the following operation

[Y i ] m = j ∈N i m [L i j ] mm [X j ] m = - j ∈N i m 4π 2 + 1 Ng n=1 ω n Y m (s n ) ω i j η (x i n ) r j (x i n )Y j ,m (x i n ) [X j ] m . ( 39
)
We implement this operation as follows. First, we compute the , m sum, by defining the intermediate quantity, for a given pair of spheres i and j ∈ N i ,

v i j n = m 4π 2 + 1 r j (x i n )Y j ,m (x i n ) [X j ] m . (40)
We accumulate in v i n the contributions of all neighbors j :

v i n = - j ∈N i v i j n ω i j η (x i n ) (41)
Finally, we take care of the sum over n:

[Y i ] m = Ng n=1 ω n Y m (s n )v i n . (42)
The first operation (eq. 40) is the evaluation of the potential produced by an internal multipolar distribution in the center of the j -th sphere at the points on the surface of the i -th sphere that are buried inside sphere j . Note that only these points need to be considered, as in eq. 41, v i j n is multiplied by the internal characteristic function ω i j η (x i n ). The latter is computed as follows

ω i j η (x i n ) = d i (x i n )χ η (r j 1 (x i n )), d i n = min(f i n , 1) f i n , (43)
where f i n is defined in eq. 37. We note that for each sphere, the matrix vector multiplication requires O(( max + 1) 2 N g ) operations, where the prefactor depends on the number of neighbors of the sphere. Therefore, the global cost of the matrix-vector multiplication is O(N ). Furthermore, the memory requirements for the computation are limited to the v i n intermediate, plus some O(( max + 1) 2 ) scratch space for the spherical harmonics. Thus, memory requirements are also linear in the size of the system. In our implementation, a loop over the spheres is the most external structure. For each sphere, a routine (calcv) is called to perform the operations in eqs. 40 and 41, then a second routine (intrhs) is called to perform the operations in eq. 42. The external loop can be trivially parallelized, as each sphere is treated independently.

| ddCOSMO adjoint linear system

The solution to the adjoint linear system can be obtained with the same solver used for the ddCOSMO equation. Again, the user must supply a routine that performs the matrix-vector multiplication with the transposed matrix, i.e., that computes

[Y i ] m = j ∈N i m [L j i ] m m [X j ] m = - j ∈N i m 4π 2 + 1 Ng n=1 ω n Y m (s n ) ω j i η (x j n ) r i (x j n )Y i ,m (x j n ) [X j ] m . ( 44
)
We proceed as follows. We first carry out the , m sum, by assembling

a j n = m Y m (s n ) [X j ] m . (45)
Then, we take care of the sum over n:

[u i j ] m = - n ω n ω j i η (x j n ) r i (x j n )Y i ,m (x j n ) a j n . (46)
Finally, we accumulate over the neighbors j ∈ N i :

[Y i ] m = j ∈N i 4π 2 + 1 [u i j ] m . (47)
Again, we note that the cost of the operations required to assemble the matrix vector product is, for each sphere, O(( max + 1) 2 N g ) operations, where the prefactor depends on the number of neighbors. Therefore, the overall matrixvector multiplication can be performed with O(N ) operations. As for memory, only the intermediates a i n and some scratch space for computing the spherical harmonics are required. The computation is performed via an external loop over the spheres. For each sphere, the intermediate a is computed and then a second routine (adjrhs) is called to take care of the operations in eq. 37 and 47. As for the ddCOSMO matrix-vector multiplication, the external loop is trivially parallelizable, as each sphere is treated independently

| ddCOSMO specific contributions to the forces

When computing the ddCOSMO forces, there are two different contributions that are specific to ddCOSMO. The first stems from the derivatives of the ddCOSMO matrix, the second from the derivatives of the U i external characteristic function that appears in the right-hand side of the ddCOSMO equation. The non ddCOSMO-specific contributions have already been discussed in section 3.4. The contributions to the force on the i -th atom stemming from the gradients of the ddCOSMO matrix are

K i = j m [S j ] m k ∈N j m ( i [L j k ] mm )[X k ] m . (48)
We require that j = k or k ∈ N j as [L j k ] mm = 0 otherwise. First, if j = k , the diagonal blocks of the L matrix do not depend on the positions of the nuclei, so no contribution to the forces arises. Second, for k ∈ N j , it is convenient to distinguish three different cases:

• i = j : We refer to this contribution as K A i ; • i = k : We refer to this contribution as K B i ;

• i ∈ N j but i k : This last contribution is referred to as K C i .
The derivation is very cumbersome and can be found in ref. [START_REF] Lipparini | Fast domain decomposition algorithm for continuum solvation models: Energy and first derivatives[END_REF]. We report here the final expressions. For the first contribution, we get

K A i = - m [S i ] m k ∈N i m 4π 2 + 1 Ng n=1 ω n Y m (s n ) 1 r k ω i k η (x i n )δ l ≥1 r k -1 (x i n ) Y k ,m (x i n )e k (x i n ) + i Y k ,m (x i n ) + 1 r k r k (x i n )Y k ,m (x i n )d i (x i n )χ η (r k 1 (x i n ))e k (x i n ) -δ f i n >1 r k (x i n )Y k ,m (x i n )d i (x i n )ω i k η (x i n )Z i n [X k ] m , (49)
where we introduced the notation

δ x >y = 1 x > y , 0 x ≤ y .
1η and 1. Finally, we accumulate over n, skipping the points for which we know a priori that there is no contribution:

K A i = - Ng n=1 ω n a i n [v a ] i n . ( 55
)
This contribution is handled by the fdoka subroutine in the ddCOSMO module. The cost of this contribution is, for each atom, O(( max + 1) 2 N g ), where the prefactor depends on the number of neighbors.

The second and third contributions are handled as follows. First, we assemble the intermediates for all neighbors j of the i -th sphere,

α j i n = m δ l ≥1 4π 2 + 1 r i -1 (x j n ) Y i ,m (x j n )e i (x j n ) + i Y i ,m (x j n ) [X i ] m , (56)
β j i n = m 4π 2 + 1 r i (x j n )Y i ,m (x j n )[X i ] m . (57)
While we are looping over the grid points and the neighbors j to i , we compute f j n . If such a quantity is larger than one, we loop over the neighbors k to j that are different from i and compute r k 1 (x

j n ). If 1 -η < r k 1 (x j n ) < 1,
there is a contribution from eq. 51. In order to take care of the latter, we assemble the intermediate

β j k n = m 4π 2 + 1 r k (x j n )Y k ,m (x j n )[X k ] m (58)
and we accumulate

[v c ] i j n = -δ f j n >1 k ∈N j \{i } 1 r i d j (x j n )ω j k η (x j n )χ η (r k 1 (x j n ))e k (x j n )β j k n . (59)
Then, we accumulate over the neighbors j

[v B ] i n = j ∈N i 1 r i ω j i η (x j n )[α j i n ] + 1 -δ f j n >1 ω j i η (x j n ) d j (x j n )χ η (r i 1 (x j n ))e i (x j n )β j i n + [v c ] j i n a j n . (60)
Finally, we accumulate over n and get

K B i + K C i = n ω n [v B ] i n . (61)
Again, it is important to note that many terms vanish and can be skipped using the same arguments used for the first one. These contributions are computed by the fdokb subroutine in the ddCOSMO module. The cost of this contribution is, for each atom, O(( max + 1) 2 N g ) operations, where the prefactor depends on the number of neighbors. This second computation is overall slightly more expensive than the first one, although it maintains the same scaling. This concludes the computation of the contribution arising from the derivatives of the ddCOSMO matrix.

The last contribution to the ddCOSMO forces is

G i = - m j Ng n=1 ω n Y m (s n ) ( i U j η (x j n ))Φ(x j n )[S j ] m . (62)
F I G U R E 1 Timings (in seconds) for the iterative solution of the ddCOSMO equation (direct), adjoint equation (adjoint) and for the computation of the ddCOSMO contribution to the forces (forces). Computations performed on a desktop computer equipped with an Intel i7-Q7700K quad-core processor running at 4.20GHz and 16GB of RAM. The following discretization parameters have been used: max = 6, N grid = 110, η = 0.2. We used a convergence threshold of 10 -7 for the iterative solver. number of atoms direct adjoint forces point of view. Second, thanks to the regularization adopted, the continuity of the energy and forces with respect to the position of the nuclei are guaranteed, which allows for energy-conserving molecular dynamics simulations and smooth geometry optimizations. Third, it is an intrinsically linear-scaling algorithm and does not require the use of fast summation techniques, that complicate the implementation and introduce a potential source of numerical errors.

Finally, and most importantly, it is overall much faster -up to three orders of magnitude -with respect to any other existing implementation. While an extensive benchmark of ddCOSMO is beyond the aims of this software review and can be found elsewhere, we report in figure 1 the timings to solve the ddCOSMO equation and adjoint equation and to compute the ddCOSMO contribution to the forces for larger and larger systems. For this example, we choose the first eight chains of the GFP-like fluorescent protein DRONPA (PDB reference 2GX2), which is composed of three tetramers, each in turn made of four chains. The calculations have been performed on a desktop computer equipped with an Intel Core i7-Q7700K quad-core processor running at 4.20GHz and equipped with 16GB of RAM. We would like to stress that, to the best of our knowledge, COSMO computations on systems as large as the ones here reported have never been attempted with any other implementation, as they would be far too expensive and memory demanding. With ddCOSMO, such computations can be routinely performed on a standard desktop computer.

In conclusion, ddCOSMO is a new, powerful algorithm that extends the applicability of polarizable continuum solva-tion methods to large and very large systems. Currently, it has been implemented for molecular dynamics simulation with standard and polarizable force fields [START_REF] Lipparini | Polarizable molecular dynamics in a polarizable continuum solvent[END_REF] and for QM and QM/MM calculations [START_REF] Lipparini | Quantum calculations in solution for large to very large molecules: A new linear scaling QM/continuum approach[END_REF][START_REF] Lipparini | Quantum, classical, and hybrid QM/MM calculations in solution: General implementation of the ddCOSMO linear scaling strategy[END_REF] of the energy, structure, electric response properties and excitation energies. A distributed memory parallel implementation of ddCOSMO is available in the highly scalable Tinker HP molecular dynamics software [START_REF] Lagardere | Tinker-HP: a massively parallel molecular dynamics package for multiscale simulations of large complex systems with advanced point dipole polarizable force fields[END_REF].

The domain decomposition paradigm has recently been extended [START_REF] Stamm | A new discretization for the polarizable continuum model within the domain decomposition paradigm[END_REF][START_REF] Gatto | Computation of Forces arising from the Polarizable Continuum Model within the Domain-Decomposition Paradigm[END_REF] to the Integral Equation Formalism -Polarizable Continuum Model (IEF-PCM) [START_REF] Cancès | A new integral equation formalism for the polarizable continuum model: Theoretical background and applications to isotropic and anisotropic dielectrics[END_REF] and to the solution of the linearized Poisson-Boltzmann equation [START_REF] Quan | A Domain Decomposition Method for the Poisson-Boltzmann Solvation Models[END_REF].

Extensions of the domain decomposition paradigm to the Solvent Excluded Surface (SES) type of molecular surface [START_REF] Quan | Mathematical analysis and calculation of molecular surfaces[END_REF][START_REF] Quan | Meshing molecular surfaces based on analytical implicit representation[END_REF] has also been developed [START_REF] Quan | Polarizable Continuum Model based on the Solvent Excluded Surface[END_REF].

The extension of ddCOSMO and its generalizations to other molecular properties is being actively investigated in our groups. Better and more efficient parallel implementations are also being investigated, including a CPU implementation based on MPI and a GPU implementation.

A C K N O W L E D G E M E N T S

R E F E R E N C E S

  to the i -th atom. The symbol m denotes max =0 m=-representing a truncated series of real spherical harmonics Y m . Since W i is by construction harmonic in Ω i , we only need to match the boundary conditions on each ∂ Ω i , which is done by multiplying (2)-(3) on any ∂ Ω i by any spherical harmonics Y i

1

 1 Suggested default values for the ddCOSMO parameters for QM and classical (

  Stamm acknowledges the funding from the German Academic Exchange Service (DAAD) from funds of the "Bundesministeriums für Bildung und Forschung" (BMBF) for the project Aa-Par-T (Project-ID 57317909). Yvon Maday acknowledges the funding from the PICS-CNRS and the PHC PROCOPE 2017 (Project N • 37855ZK). This work was made possible thanks to the French state funds managed by the CalSimLab LABEX and the ANR within the Investissements d'Avenir program (reference ANR-11-LABX-0037-01).

Available at https://github.com/filippolipparini/ddPCM

The second and third contributions read

We implement these contributions as follows. All three terms share the fact that the sum over , m can be precomputed.

As a consequence, we assemble before the computation of the forces a vector

We assemble separately the contribution in eq. 49 and the ones from 50 and eq. 51, the latter being treated together.

Let us start with eq. 49. We take care of the , m sums first, by assembling, for each grid point and each neighbor k of the i -th sphere,

All contributions in eq. 52 and 53 contain quantities that have already been used for the computation of a matrix-vector product for the ddCOSMO linear system, with the exception of the spherical harmonics gradients. We compute all the ( max + 1) 2 spherical harmonics and their gradients at the same time, in order to reuse intermediates as much as possible, with an algorithm described in ref. [START_REF] Lipparini | Fast domain decomposition algorithm for continuum solvation models: Energy and first derivatives[END_REF] and implemented in the dbasis routine. We emphasize again the key importance of this step for the overall efficiency of the computation, as it has a tremendous impact on the performance of the implementation. Note that the computation of β corresponds to computing the potential of the multipolar distribution [X k ] at the point x i n . Such an operation is also required for the ddCOSMO matrix-vector products. Once the intermediate quantities α i k n (a vector of dimension 3) and β i k n (a scalar) have been computed, we accumulate over the neighbors k , introducing the vector (dimension 3) [v A ] i n defined as follows:

It is particularly important to stress that the contributions in eq. 54 are very often zero. In particular, no contribution arises for points on sphere i that are not buried in sphere k . Furthermore, the derivative of the χ η (t ) switching functions are nonzero only if 1η < t < 1. Therefore, many points in the n loop can be skipped, either by checking whether

which is equivalent to requiring that the point is buried, and by checking whether

Using the same a j n intermediate used for the other contributions to the forces, we get

The gradient of U j n is nonzero only if j = i or j ∈ N i . In the former case, we get

In the latter case, we get

We can now accumulate

We finally accumulate over n and get

The computation of this contribution is carried out by the fdoga subroutine. Its computational cost is O(N g ), with a proportionality factor that depends on the number of neighbors.

We conclude this section by pointing out that the computation is performed independently for each atom (sphere), so that it can be trivially parallelized. Furthermore, the overall computation scales linearly with the number of atoms.

| CONCLUSIONS AND PERSPECTIVES

In this Software news report, we have briefly recapitulated the derivation and described in detail the implementation of ddCOSMO, a new paradigm for the solution of the Conductor-like Screening Model equations. We have provided a complete description of our ddCOSMO implementation and commented extensively on efficiency and memory issues.

Furthermore, by making explicit reference to our implementation we have provided an extensive documentation and guidelines to interface it with an existing code, either for quantum chemistry or classical molecular dynamics simulation.

The modular ddCOSMO code is written in Fortran 90. Interfacing it with existing Fortran codes is thus straightforward and can be done by following the guidelines given in section 3. The same guidelines can be used, with some additional technical work, to interface the ddCOSMO module with codes written in a different language such as C, C++, python and others. The detailed description of the implementation provided in section 4, together with the existing Fortran code that can be used as a reference, are meant to provide a thorough guide in support of a new implementation, possibly in a different language.

The implementation of ddCOSMO is not a straightforward task. However, the new domain-decomposition paradigm offers several advantages with respect to other implementations. First, the discretization is systematically improvable and controlled by a very limited number of parameters, which makes it easy to have a detailed control from a numerical


First investigations on the feasibility of integration of a smart sensor in harsh environment 1 S. Ben-Krit, 1 W. Rahajandraibe, 1 This very first study takes place in the framework of the I_SMART European project. The final system will consist of one or more detectors in silicon carbide (spectroscopic detector, neutron, gamma …), a readout electronic, a conditioning chain and microprocessor for the signal recognition. The work presented here is dedicated to a first approach of the feasibility of integration of the full system. This system will have to to work in harsh environment in terms of temperature and radiations what makes necessary the development of specifications for operation and reliability of the components and the investigation of margins for the interplay of the components.

When an incident particle, emitted by a source of radiation tracks a semi-conductor (p-n junction for example), direct generated ions or secondary products from proton or neutron/nucleus interactions generate electron-holes pairs [START_REF] Castellani-Coulié | Development of a CMOS Oscillator Concept for Particle Detection and Tracking[END_REF]. Then, electrons drift within the semi-conductor and diffuse towards a collecting anode inducing a parasitic current. Such kind of behavior is called single effect which differs from cumulative effects (dose effects) that are linked to the presence of oxide in the component [START_REF] Gatti | The Quest for Low-Loise Processing of Random Signals, Integrated Circuits, and Nuclear Science[END_REF]. Thus, radiation sensors require pulse processing electronics so that energy or time information involved with radiation interactions (parasitic currents) can be properly extracted.

Several types of detectors can be found in literature. Depending on the needs, they can be divided on two broad categories. A lot of detectors are based on SRAM detection which is only able to ensure information about the number of incident particles detected by the sensor.

Other detectors are based on collection and amplification by the use of the silicon sensor (p-n junction) as volume in order to provide information about the generated currents characteristics. Generally, the collection is assured by a matrix of pixels to ensure high sensitivity [START_REF] Castellani-Coulié | Development of a CMOS Oscillator Concept for Particle Detection and Tracking[END_REF]. Although this second type requires a whole analog preprocessing chain comprising several components, it is still the more required one due to the importance of information that can be provided on its outputs. The readout chain is mainly composed by i) blocks ensuring amplification, based on a charge sensitive amplifier (CSA) or a current sensitive amplifier (ISA), ii) pulse-shaping that is assured by a pulse shaping amplifier (PSA) providing a semi-Gaussian pulse and iii) an analog to digital conversion (ADC) which provides a digital data ready to pass on to iv) digital signal processing (DSP) [3] [4].

Since the readout electronics has to be placed as close as possible to the sensor to ensure high sensitivity, the full system must be able to work under several constraints (very high temperature and irradiations). However, providing electronic components able to run normally under a temperature higher than 400° and under radiative environment is still a real challenge in the area of the electronic instrumentation. In this context, this first work aims to provide first investigations on the feasibility of integration of the full system in these conditions. It focuses essentially on the development of the specifications for operation (signal to noise ratio (SNR), bandwidth, gain…) of the components. An example of the proposed chain can be shown in the following figure (Fig. 1). Moreover, several solutions will be proposed to decrease pulse pile-up occurrence. Pulse pile-up happens when pulses arrives closer in time than the pulse resolution time for the system. Thus, in Fig. 3, the pile-up phenomenon is visible. The peak pulse of the output signal is 1.1V at the beginning of the process (before pile-up) and increase up to 1.33V as soon as the second pulse is detected. In the final paper, the input currents will be noise source currents modeled using the systematic tool SITARe [START_REF] Micolau | SITARe: a SImulation Tool for Analysis and diagnosis of Radiation Effects[END_REF]. In SITARe, the current model takes as input the ion physical characteristics (initial energy, LET, particle species) and location (point of generation, strike angle). The current modeling is based on a simplified drift-diffusion model.

The model outputs are the currents generated at the sensitive zones of the pixel of detection, what are the outputs of the detection chain.

Fig. 1 :!Fig. 2 .Fig. 2 :

 122 Fig. 1: Electrical circuit diagram of the analog signal conditioning (1 channel)

!

  "#$%&'$(()"%* +,-&.$%#$%* +/01"&.$%#$%* 2.3%04)*526* !*5$-6* 718)*5$96* Finally, the association of the systematic tool SITARe with Cadence will provide a complete simulation platform, necessary for the characterization and the design of the full conditioning chain.

Fig. 3 :

 3 Fig. 3: Transient response of the conditioning chain showing pulse pile-up (Frequency=572 KHz)

K. Coulié- Castellani, 2 G. Micolau

	1 IM2NP-UMR CNRS 7334 / Aix-Marseille University
	IM2NP, IMT Technopôle de Château -Gombert, 13451 Marseille Cedex 20, France
	2 UMR 1114 EMMAH (INRA-UAPV)