Cálculo de fuerzas de largo alcance en Dinámica Molecular - Suma de Ewald

Estoy tratando de escribir un código para calcular el potencial y las fuerzas, para lo mismo usando la suma de Ewald. Para este propósito, la fórmula para el potencial y la fuerza que he usado es:

tu = tu ( r ) + tu ( k ) + tu ( b C ) + tu s mi yo F

donde la contribución del potencial en el espacio k viene dada por

tu ( k ) = 1 2 π L 3 k 0 4 π 2 k 2 mi k 2 4 k 2 | S ( k ) | 2 S ( k ) = i = 1 norte z i mi i k . r i
tu ( r ) = j < i norte = 0 z i z j erfc ( k | r i j + norte | ) | r i j + norte |
tu ( b C ) = 2 π 3 L 3 | i = 1 norte z i r i | 2 tu ( s mi yo F ) = k π i = 1 norte z i 2

y las ecuaciones de fuerza -

F i = F i ( k ) + F i ( r ) + F i ( b C )

F i ( k ) = 4 π z i 2 π L 3 k 0 k mi k 2 4 k 2 k 2 ( pecado ( k i . r i ) Re ( S ( k ) ) + porque ( k i . r i ) Soy ( S ( k ) ) )
y dos ecuaciones más, que estoy cansado de escribir pero estoy seguro de que son correctas.

Actualmente estoy usando este método para resolver la simulación de dinámica molecular de moléculas de agua. El problema al que me enfrento es que el cálculo del potencial depende en gran medida de k , que puede ser, pero hay un gran cambio en el potencial si cambio k levemente. En segundo lugar, ¿no debería el cambio con respecto disminuir más allá de cierto valor? No veo que eso también suceda.

Aquí deseo saber si hay una forma precisa de averiguar si mi código basado en las fórmulas anteriores funciona correctamente.

EDICIÓN 1: Como nota, he simulado el modelo de agua SPC/E usando esto ahora. El problema al que me enfrento en este momento es que la energía potencial que se calcula es de 5 a 6 veces mayor que la que se informó anteriormente en las revistas, pero sin embargo, la conservación de energía durante la ejecución está intacta. No entiendo la razón aquí !!

Respuestas (1)

Solo he usado la suma de Ewald, nunca la he implementado yo mismo.

Sin embargo, mencionas que no estás convergiendo como k aumenta ni está convergiendo al valor correcto. Parecería que, independientemente del problema, si su implementación es correcta, debería converger en algún punto.

Si alcanza la convergencia wrt k ; en cuanto al punto de que no está obteniendo los valores informados. Recuerda que el tamaño del cuadro de simulación también es importante. ¿Está seguro de que está utilizando un espacio de simulación lo suficientemente grande con suficientes moléculas de agua?

EDITAR:

Basado en sus comentarios a continuación:

Bueno, la energía debe conservarse, :P. Si su campo de fuerza, sea lo que sea que implementó, es conservador, entonces la energía debe conservarse (o su integrador de movimiento está mal).

¿Qué conjunto estás usando? ¿El papel? Si es NVE, ¿está seguro de que la V es exactamente correcta? recuerde que incluso un cambio del 0,1% en el tamaño de la celda de la caja provocará un gran cambio en la energía de un material condensado. (¡intenta comprimir agua al 0,1 %!)

¿Estás simulando agua fluida o agua cristalina? De cualquier manera, ¿estás seguro de que no estás atrapado en un mínimo local? Intenta pulsar el volumen de la caja del sim para salir.

¡Estoy usando un tamaño de caja de simulación de 18,63 angstroms con 216 moléculas! ¡He usado los detalles dados de un documento antiguo, truncando la suma de Ewald por la condición mínima de imagen! ¿Cómo puedo depurar mi código? ¿Por qué la conservación de la energía sigue intacta?
+1 Gracias por las ideas. Solo estoy usando el conjunto NVE, el artículo al que me refiero es uno de M. Prevost et al. . Una parte que no entendí sobre el documento fue que mencionó el uso de 342 vectores en el espacio recíproco con una magnitud de 3 en la suma sobre la suma de la red recíproca. Ahora no entendí a qué magnitud se refiere, así que no pude adherirme estrictamente a la magnitud de los k-vectores, ¡aunque lo mantuve por debajo de un límite superior!