~jose-soler/siesta/newdynamics : revision 745

1

!*****************************************************************************

2

module md_state_module

3

!*****************************************************************************

4

!

5

! Project: MolecularDynamics

6

!

7

! Created on: Mon 4 Jun 10:42:01 2018 by Eduardo R. Hernandez

8

!

9

! Last Modified: Mon 5 Nov 18:15:31 2018

10

!

11

! ---

12

13

! This file is distributed under the terms of the

14

! GNU General Public License: see COPYING in the top directory

15

! or http://www.gnu.org/copyleft/gpl.txt .

16

! ---

17

!

18

! This module defines the State derived type. This type contains information

19

! about the atomic positions, velocities (momenta), cell, extended variables

20

! (thermostat and/or barostat) etc. In short, it contains all the information

21

! needed to characterise the state of a simulated system, bundled into one

22

! object. All the data in State is private, and must be accessed through

23

! the corresponding get() and put() routines, also provided here.

24

!

25

! Internally, positions, momenta and forces are usually stored most commonly

26

! in Cartesian representation; however, some MD integration algorithms,

27

! in particular those that involve constant-pressure, require them in a lattice

28

! representation. It is therefore possible to transform from one representation

29

! to the other; the transformation is different and independent for each type

30

! of variable. A logical flag is used to keep track of the current representation

31

! of each variable.

32

!

33

!*****************************************************************************

34

! modules used

35

36

use md_precision_module

37

use md_constants_module

38

39

!*****************************************************************************

40

! No implicits please!

41

42

implicit none

43

44

!*****************************************************************************

45

! derived type definitions

46

47

type barostat_type

48

49

real (dp) :: mass

50

! there is no position; that role is played by the cell metricTensor

51

real (dp), dimension (3,3) :: momentum

52

53

end type barostat_type

54

55

type energy_type

56

57

real (dp) :: kinetic ! thermostat is used

58

real (dp) :: potential

59

real (dp) :: thermostat

60

real (dp) :: barostat

61

real (dp) :: constantOfMotion ! may be different from total energy

62

63

end type energy_type

64

65

type energyOrigin_type

66

67

logical :: set ! true if origin is set

68

real (dp) :: E

69

70

end type energyOrigin_type

71

72

type lattice_type

73

74

real (dp), dimension(3) :: a ! lattice vectors

75

real (dp), dimension(3) :: b

76

real (dp), dimension(3) :: c

77

real (dp) :: a_modulus ! their modulus

78

real (dp) :: b_modulus

79

real (dp) :: c_modulus

80

real (dp) :: alpha ! their angles

81

real (dp) :: beta

82

real (dp) :: gamma

83

real (dp) :: volume ! cell volume

84

real (dp), dimension (3,3) :: metricTensor

85

86

end type lattice_type

87

88

type representation_type

89

90

logical :: position ! .true. if Cartesian rep; .false. otherwise

91

logical :: momentum

92

logical :: force

93

94

end type representation_type

95

96

type stress_type

97

98

logical :: calculated

99

real (dp), dimension(3,3) :: cartesian

100

real (dp), dimension(3,3) :: lattice

101

102

end type stress_type

103

104

type thermostat_type

105

106

real (dp) :: mass

107

real (dp) :: position

108

real (dp) :: momentum

109

110

! the following are parameters used only

111

! in the particular case of Recursive-Nosé-Poincaré chains

112

113

real (dp) :: a

114

real (dp) :: C

115

116

end type thermostat_type

117

118

type state_type

119

120

type ( representation_type ), private :: representation

121

122

! atom-related variables

123

124

integer, private :: nAtoms

125

integer, private :: nDegreesOfFreedom

126

real (dp), private, pointer, dimension (:,:) :: position

127

real (dp), private, pointer, dimension (:,:) :: force

128

real (dp), private, pointer, dimension (:,:) :: momentum

129

real (dp), private, pointer, dimension (:) :: mass

130

131

! cell-related variables

132

type ( lattice_type ), private :: Cell

133

type ( lattice_type ), private :: reciprocalCell

134

type ( stress_type ), private :: stress

135

136

! extended-system variables, if any

137

integer :: nThermostats

138

type ( thermostat_type ), private, pointer, dimension (:) :: thermostat

139

type ( barostat_type ), private :: barostat

140

141

! energy of the system

142

143

type ( energy_type ), private :: energy

144

145

type ( energyOrigin_type ), private :: H0 ! constant, only used if Nose thermostat(s) used

146

147

end type state_type

148

149

!*****************************************************************************

150

! private module variables

151

152

logical, private, parameter :: debug = .false.

153

logical, private :: PBC ! periodic boundary conditions (true if)

154

155

real (dp) :: energyConversionFactor

156

real (dp), private :: forceConversionFactor

157

real (dp), private :: lengthConversionFactor

158

real (dp), private :: massConversionFactor

159

real (dp) :: pressureConversionFactor

160

real (dp) :: timeConversionFactor

161

real (dp), private :: velConversionFactor

162

163

type ( state_type ), private :: state ! the current state of the system

164

165

!*****************************************************************************

166

! private module subroutines

167

168

private generateMomenta

169

private calculateKineticEnergy

170

private transformCoordinates2Cartesian

171

private transformCoordinates2Lattice

172

private transformForces2Cartesian

173

private transformForces2Lattice

174

private transformMomenta2Cartesian

175

private transformMomenta2Lattice

176

177

contains

178

179

!*****************************************************************************

180

subroutine calculateKineticEnergy( state )

181

!*****************************************************************************

182

!

183

! Project: MolecularDynamics

184

!

185

! Module: state

186

!

187

! Created on: Wed 27 Jun 17:37:12 2018 by Eduardo R. Hernandez

188

!

189

! Last Modified: Wed 18 Jul 17:16:46 2018

190

!

191

! ---

192

193

! This file is distributed under the terms of the

194

! GNU General Public License: see COPYING in the top directory

195

! or http://www.gnu.org/copyleft/gpl.txt .

196

! ---

197

!

198

! Given the momentum contained in state (and if appropriate, any thermostats)

199

! this routine calculates the kinetic energy of the system.

200

!

201

!*****************************************************************************

202

! modules used

203

204

!*****************************************************************************

205

! No implicits please!

206

207

implicit none

208

209

!*****************************************************************************

210

! variables

211

212

type ( state_type ), intent (INOUT) :: state

213

214

!*****************************************************************************

215

! variables

216

217

integer :: i, n

218

219

real (dp) E_kin, factor, prefactor

220

221

!*****************************************************************************

222

! Start of subroutine

223

224

if ( debug ) write(*,*) 'Entering calculateKineticEnergy()'

225

226

! consider the different possible cases regarding thermostats

227

228

E_kin = zero

229

230

! make sure momenta are in Cartesian coordinates

231

232

call transformMomenta2Cartesian( state )

233

234

do i=1, state % nAtoms

235

236

factor = half / state % mass(i)

237

238

E_kin = E_kin + factor * &

239

dot_product( state % momentum(1:3,i), state % momentum(1:3,i) )

240

241

end do

242

243

prefactor = one

244

245

if ( state % nThermostats > 0 ) then ! this is the standard NVE case

246

247

do n=1, state % nThermostats

248

249

prefactor = prefactor / state % thermostat(n) % position

250

251

end do

252

253

prefactor = prefactor * prefactor

254

255

end if

256

257

state % energy % kinetic = prefactor * E_kin

258

259

if ( debug ) write(*,*) 'Exiting calculateKineticEnergy()'

260

261

end subroutine calculateKineticEnergy

262

263

!*****************************************************************************

264

subroutine createState( state, nAtoms, &

265

nThermostats, thermostatMass, barostatMass )

266

!*****************************************************************************

267

!

268

! Project: MolecularDynamics

269

!

270

! Module: State

271

!

272

! Created on: Tue 12 Jun 11:42:03 2018 by Eduardo R. Hernandez

273

!

274

! Last Modified: Mon 5 Nov 18:17:18 2018

275

!

276

! ---

277

278

! This file is distributed under the terms of the

279

! GNU General Public License: see COPYING in the top directory

280

! or http://www.gnu.org/copyleft/gpl.txt .

281

! ---

282

!

283

! This subroutine is called to create a new instance of State; the state instance

284

! is allocated, but not "filled"; this must be done by either initialiseState (in

285

! case of a brand new simulation) or by readRestart (in case of a continuation

286

! simulation).

287

!

288

! variables:

289

!

290

! type ( state_type ), intent (INOUT) :: state

291

! this is the state that is going to be allocated and filled with the data

292

! passed next.

293

!

294

! integer, intent (IN) :: nAtoms

295

! the number of atoms in the system; used to dimension several arrays below

296

!

297

! integer, intent (IN), optional :: nThermostats

298

! the number of thermostats in the thermostat chain. For NVE or NPH simulations

299

! that have no thermostat, this variable is either set to zero or absent (default = 0);

300

! for simulations with a thermostat or a thermostat chain (NVT, NPT), nThermostats

301

! is a positive integer indicating the number of Nosé-Poincaré thermostats in the

302

! chain.

303

!

304

! real (dp), dimension (:), intent (IN), optional :: thermostatMass

305

! if nThermostats > 0, thermostatMass is an array of length nThermostats providing the

306

! masses of each thermostat in the chain, in Daltons (1/12 of the mass of C12)

307

!

308

! real (dp), intent (IN), optional :: barostatMass

309

! if the simulation is a constant-pressure one, this specifies the mass of the barostat

310

! in Daltons (1/12 of the mass of C12)

311

!

312

!*****************************************************************************

313

! modules used

314

315

!*****************************************************************************

316

! No implicits please!

317

318

implicit none

319

320

!*****************************************************************************

321

! arguments

322

323

type ( state_type ), intent (INOUT) :: state

324

325

integer, intent (IN) :: nAtoms

326

327

integer, intent (IN), optional :: nThermostats

328

329

real (dp), dimension (:), intent (IN), optional :: thermostatMass

330

real (dp), intent (IN), optional :: barostatMass

331

332

!*****************************************************************************

333

! local variables

334

335

integer :: i

336

337

real (dp), dimension (3,3) :: LatticeVectors

338

339

real (dp), allocatable, dimension (:,:) :: momentum

340

341

!*****************************************************************************

342

! start subroutine

343

344

if ( debug ) write(*,*) 'Entering createState()'

345

346

! assign nAtoms and allocate variables that depend on it

347

348

state % nAtoms = nAtoms

349

state % nDegreesOfFreedom = 3 * state % nAtoms - 3

350

351

allocate( state % position( 3, state % nAtoms ) )

352

allocate( state % force( 3, state % nAtoms ) )

353

allocate( state % momentum( 3, state % nAtoms ) )

354

allocate( state % mass( state % nAtoms ) )

355

356

if ( present( nThermostats ) ) then

357

358

state % nThermostats = nThermostats

359

360

allocate( state % thermostat( nThermostats ) )

361

362

do i = 1, state % nThermostats

363

state % thermostat(i) % mass = massConversionFactor * thermostatMass(i)

364

state % thermostat(i) % position = one

365

state % thermostat(i) % momentum = zero

366

state % thermostat(i) % a = one

367

state % thermostat(i) % C = one

368

write(*,*) 'thermostat mass ', i, state % thermostat(i) % mass

369

end do

370

371

else

372

373

state % nThermostats = 0 ! no thermostat is defined

374

375

end if

376

377

if ( present( barostatMass ) ) then

378

379

state % barostat % mass = massConversionFactor * barostatMass

380

381

state % barostat % momentum = zero

382

383

end if

384

385

! these are initialised to zero; they will be set prorperly if needed later

386

387

state % energy % thermostat = zero

388

state % energy % barostat = zero

389

390

state % H0 % set = .false.

391

state % H0 % E = zero

392

393

if ( debug ) write(*,*) 'Exiting createState'

394

395

end subroutine createState

396

397

!*****************************************************************************

398

subroutine deleteState( state )

399

!*****************************************************************************

400

!

401

! Project: MolecularDynamics

402

!

403

! Module: State

404

!

405

! Created on: Thu 19 Jul 11:28:36 2018 by Eduardo R. Hernandez

406

!

407

! Last Modified: Fri 2 Nov 17:27:20 2018

408

!

409

! ---

410

411

! This file is distributed under the terms of the

412

! GNU General Public License: see COPYING in the top directory

413

! or http://www.gnu.org/copyleft/gpl.txt .

414

! ---

415

!

416

! type ( state_type ), intent (INOUT) :: state

417

! this is the state that is going to be deallocated.

418

!

419

!*****************************************************************************

420

! modules used

421

422

!*****************************************************************************

423

! No implicits please!

424

425

implicit none

426

427

!*****************************************************************************

428

! arguments

429

430

type ( state_type ), intent (INOUT) :: state

431

432

!*****************************************************************************

433

! local variables

434

435

!*****************************************************************************

436

! start subroutine

437

438

if ( debug ) write(*,*) 'Entering deleteState()'

439

440

deallocate( state % position )

441

deallocate( state % force )

442

deallocate( state % momentum )

443

deallocate( state % mass )

444

if ( state % nThermostats .gt. 0 ) deallocate( state % thermostat )

445

446

if ( debug ) write(*,*) 'Exiting deleteState'

447

448

end subroutine deleteState

449

450

!*****************************************************************************

451

subroutine generateMomenta( state, target_temperature )

452

!*****************************************************************************

453

!

454

! Project: MolecularDynamics

455

!

456

! Created on: Wed 11 Apr 12:27:49 2018 by Eduardo R. Hernandez

457

!

458

! Last Modified: Thu 26 Jul 16:37:05 2018

459

!

460

! ---

461

462

! This file is distributed under the terms of the

463

! GNU General Public License: see COPYING in the top directory

464

! or http://www.gnu.org/copyleft/gpl.txt .

465

! ---

466

!

467

!*****************************************************************************

468

! modules used

469

470

!*****************************************************************************

471

! No implicits please!

472

473

implicit none

474

475

!*****************************************************************************

476

! variables

477

478

type ( state_type ), intent (INOUT) :: state

479

480

real (dp), intent (IN) :: target_temperature

481

482

!**************************************************************************

483

! Local variables

484

485

integer :: i, n

486

487

real (dp) :: r, r2, rnum, total_mass, variance, v1, v2

488

real (dp) :: determinant, determinant_1, E_kin, factor, temperature

489

490

real (dp) :: angular_momentum( 3 ), angular_velocity( 3 )

491

real (dp) :: centre_of_mass( 3 )

492

real (dp) :: cm_velocity( 3 )

493

real (dp) :: inertia_tensor(3,3)

494

real (dp) :: matrix(3,3)

495

496

!**************************************************************************

497

! Start of subroutine

498

499

if ( debug ) write(*,*) 'Entering generateMomenta()'

500

501

total_mass = zero

502

503

! loop over atoms and assign each atom a random velocity

504

! according to its mass and the target temperature

505

506

! let's try our luck

507

508

do i=1, state % nAtoms

509

510

total_mass = total_mass + state % mass(i)

511

variance = dsqrt( boltzmann_k * target_temperature / state % mass(i) )

512

513

do n= 1,3

514

515

r = two

516

do while ( r .gt. one )

517

call random_number( rnum )

518

v1 = two * rnum - one

519

call random_number( rnum )

520

v2 = two * rnum - one

521

r = v1 * v1 + v2 * v2

522

end do

523

r = variance * v1 * sqrt( -two * log( r ) / r )

524

state % momentum(n,i) = state % mass(i) * r

525

end do

526

end do

527

528

! now calculate the amount of velocity to take out for each atom

529

530

do n=1, 3

531

cm_velocity( n ) = zero

532

do i=1, state % nAtoms

533

cm_velocity(n) = cm_velocity(n) + state % momentum(n,i)

534

end do

535

cm_velocity(n) = cm_velocity(n) / total_mass

536

end do

537

538

! now substract from each atom the velocity necessary so as

539

! to ensure that the total momentum is zero

540

541

do n=1, 3

542

do i=1, state % nAtoms

543

state % momentum(n,i) = state % momentum(n,i) - &

544

state % mass(i) * cm_velocity(n)

545

end do

546

end do

547

548

! calculate the instantaneous temperature

549

550

E_kin = zero

551

552

do i=1, state % nAtoms

553

554

factor = half / state % mass(i)

555

556

E_kin = E_kin + factor * ( &

557

dot_product( state % momentum(1:3,i), state % momentum(1:3,i) ) )

558

559

end do

560

561

temperature = two * E_kin / ( ( three * dfloat(state % nAtoms) - three ) * boltzmann_k )

562

563

! calculate the scaling factor and scale the velocities accordingly

564

565

write(*,*) 'temperature = ', temperature, target_temperature

566

567

factor = dsqrt( target_temperature / temperature )

568

569

state % momentum = factor * state % momentum

570

571

E_kin = zero

572

573

do i=1, state % nAtoms

574

575

factor = half / state % mass(i)

576

577

E_kin = E_kin + factor * ( &

578

dot_product( state % momentum(1:3,i), state % momentum(1:3,i) ) )

579

580

end do

581

582

temperature = two * E_kin / ( ( three * dfloat(state % nAtoms) - three ) * boltzmann_k )

583

584

! just checking

585

586

do n=1, 3

587

do i=1, state % nAtoms

588

cm_velocity(n) = cm_velocity(n) + state % momentum(n,i)

589

end do

590

cm_velocity(n) = cm_velocity(n) / total_mass

591

end do

592

593

write(*,*) 'centre of mass velocity = ', cm_velocity(1:3), temperature

594

write(*,*) 'kinetic energy = ', E_kin

595

write(*,*) 'momentum 1 ', state % momentum(1:3,1)

596

597

! if the system is not periodic in any direction (i.e. cluster or molecule)

598

! then we must make sure that it does have zero angular momentum (no rotation)

599

600

if ( .not. PBC ) then

601

602

! we will refer momentarily the positions with respect to the centre of mass

603

604

centre_of_mass = zero

605

606

do i=1, state % nAtoms

607

608

centre_of_mass(1) = centre_of_mass(1) + &

609

state % mass(i) * state % position(1,i)

610

centre_of_mass(2) = centre_of_mass(2) + &

611

state % mass(i) * state % position(2,i)

612

centre_of_mass(3) = centre_of_mass(3) + &

613

state % mass(i) * state % position(3,i)

614

615

end do

616

617

centre_of_mass = centre_of_mass / total_mass

618

619

state % position(1,1:state % nAtoms) = &

620

state % position(1,1:state % nAtoms) - centre_of_mass(1)

621

state % position(2,1:state % nAtoms) = &

622

state % position(2,1:state % nAtoms) - centre_of_mass(2)

623

state % position(3,1:state % nAtoms) = &

624

state % position(3,1:state % nAtoms) - centre_of_mass(3)

625

626

angular_momentum = zero

627

inertia_tensor = zero

628

629

do i=1, state % nAtoms

630

631

angular_momentum(1) = angular_momentum(1) + &

632

( state % position(2,i) * state % momentum(3,i) - &

633

state % position(3,i) * state % momentum(2,i) )

634

angular_momentum(2) = angular_momentum(2) + &

635

( state % position(3,i) * state % momentum(1,i) - &

636

state % position(1,i) * state % momentum(3,i) )

637

angular_momentum(3) = angular_momentum(3) + &

638

( state % position(1,i) * state % momentum(2,i) - &

639

state % position(2,i) * state % momentum(1,i) )

640

641

r2 = dot_product( state % position(1:3,i), state % position(1:3,i))

642

643

inertia_tensor(1,1) = inertia_tensor(1,1) + state % mass(i) * &

644

( r2 - state % position(1,i) * state % position(1,i) )

645

inertia_tensor(1,2) = inertia_tensor(1,2) - state % mass(i) * &

646

state % position(1,i) * state % position(2,i)

647

inertia_tensor(1,3) = inertia_tensor(1,3) - state % mass(i) * &

648

state % position(1,i) * state % position(3,i)

649

650

inertia_tensor(2,1) = inertia_tensor(2,1) - state % mass(i) * &

651

state % position(2,i) * state % position(1,i)

652

inertia_tensor(2,2) = inertia_tensor(2,2) + state % mass(i) * &

653

( r2 - state % position(2,i) * state % position(2,i) )

654

inertia_tensor(2,3) = inertia_tensor(2,3) - state % mass(i) * &

655

state % position(2,i) * state % position(3,i)

656

657

inertia_tensor(3,1) = inertia_tensor(3,1) - state % mass(i) * &

658

state % position(3,i) * state % position(1,i)

659

inertia_tensor(3,2) = inertia_tensor(3,2) - state % mass(i) * &

660

state % position(3,i) * state % position(2,i)

661

inertia_tensor(3,3) = inertia_tensor(3,3) + state % mass(i) * &

662

( r2 - state % position(3,i) * state % position(3,i) )

663

664

end do

665

666

determinant = inertia_tensor(1,1) * ( &

667

inertia_tensor(2,2) * inertia_tensor(3,3) - &

668

inertia_tensor(2,3) * inertia_tensor(3,2) ) - &

669

inertia_tensor(1,2) * ( &

670

inertia_tensor(2,1) * inertia_tensor(3,3) - &

671

inertia_tensor(2,3) * inertia_tensor(3,1) ) + &

672

inertia_tensor(1,3) * ( &

673

inertia_tensor(2,1) * inertia_tensor(3,2) - &

674

inertia_tensor(2,2) * inertia_tensor(3,1) )

675

676

! calculate the components of the angular velocity

677

678

matrix = inertia_tensor

679

matrix(1:3,1) = angular_momentum

680

681

determinant_1 = matrix(1,1) * ( &

682

matrix(2,2) * matrix(3,3) - matrix(2,3) * matrix(3,2) ) - &

683

matrix(1,2) * ( &

684

matrix(2,1) * matrix(3,3) - matrix(2,3) * matrix(3,1) ) + &

685

matrix(1,3) * ( &

686

matrix(2,1) * matrix(3,2) - matrix(2,2) * matrix(3,1) )

687

688

angular_velocity(1) = determinant_1 / determinant

689

690

matrix = inertia_tensor

691

matrix(1:3,2) = angular_momentum

692

693

determinant_1 = matrix(1,1) * ( &

694

matrix(2,2) * matrix(3,3) - matrix(2,3) * matrix(3,2) ) - &

695

matrix(1,2) * ( &

696

matrix(2,1) * matrix(3,3) - matrix(2,3) * matrix(3,1) ) + &

697

matrix(1,3) * ( &

698

matrix(2,1) * matrix(3,2) - matrix(2,2) * matrix(3,1) )

699

700

angular_velocity(2) = determinant_1 / determinant

701

702

matrix = inertia_tensor

703

matrix(1:3,3) = angular_momentum

704

705

determinant_1 = matrix(1,1) * ( &

706

matrix(2,2) * matrix(3,3) - matrix(2,3) * matrix(3,2) ) - &

707

matrix(1,2) * ( &

708

matrix(2,1) * matrix(3,3) - matrix(2,3) * matrix(3,1) ) + &

709

matrix(1,3) * ( &

710

matrix(2,1) * matrix(3,2) - matrix(2,2) * matrix(3,1) )

711

712

angular_velocity(3) = determinant_1 / determinant

713

714

angular_velocity = -angular_velocity

715

716

do i=1, state % nAtoms

717

718

state % momentum(1,i) = state % momentum(1,i) + &

719

state % mass(i) * ( &

720

state % position(3,i) * angular_velocity(2) - &

721

state % position(2,i) * angular_velocity(3) )

722

723

state % momentum(2,i) = state % momentum(2,i) + &

724

state % mass(i) * ( &

725

state % position(1,i) * angular_velocity(3) - &

726

state % position(3,i) * angular_velocity(1) )

727

728

state % momentum(3,i) = state % momentum(3,i) + &

729

state % mass(i) * ( &

730

state % position(2,i) * angular_velocity(1) - &

731

state % position(1,i) * angular_velocity(2) )

732

733

end do

734

735

! just checking

736

737

angular_momentum = zero

738

739

do i=1, state % nAtoms

740

741

angular_momentum(1) = angular_momentum(1) + &

742

( state % position(2,i) * state % momentum(3,i) - &

743

state % position(3,i) * state % momentum(2,i) )

744

angular_momentum(2) = angular_momentum(2) + &

745

( state % position(3,i) * state % momentum(1,i) - &

746

state % position(1,i) * state % momentum(3,i) )

747

angular_momentum(3) = angular_momentum(3) + &

748

( state % position(1,i) * state % momentum(2,i) - &

749

state % position(2,i) * state % momentum(1,i) )

750

751

end do

752

753

write(*,*) 'angular momentum components ', angular_momentum(1:3)

754

755

state % position(1,1:state % nAtoms) = &

756

state % position(1,1:state % nAtoms) + centre_of_mass(1)

757

state % position(2,1:state % nAtoms) = &

758

state % position(2,1:state % nAtoms) + centre_of_mass(2)

759

state % position(3,1:state % nAtoms) = &

760

state % position(3,1:state % nAtoms) + centre_of_mass(3)

761

762

end if

763

764

! finally, calculate the kinetic energy and store it in type energy

765

766

E_kin = zero

767

768

do i=1, state % nAtoms

769

770

factor = half / state % mass(i)

771

772

E_kin = E_kin + factor * ( &

773

dot_product( state % momentum(1:3,i), state % momentum(1:3,i) ) )

774

775

end do

776

777

state % energy % kinetic = E_kin

778

779

state % representation % momentum = .true. ! momenta are in cartesian rep.

780

781

if ( debug ) write(*,*) 'Exiting generateMomenta()'

782

783

end subroutine generateMomenta

784

785

!*****************************************************************************

786

subroutine get( state, nAtoms, nDegreesOfFreedom, nThermostats, &

787

position, force, momentum, mass, &

788

thermostat, barostat, energy, H0, &

789

Cell, reciprocalCell, stress, representation )

790

!*****************************************************************************

791

!

792

! Project: MolecularDynamics

793

!

794

! Module: State

795

!

796

! Created on: Fri 8 Jun 14:53:35 2018 by Eduardo R. Hernandez

797

!

798

! Last Modified: Mon 5 Nov 18:16:11 2018

799

!

800

! ---

801

802

! This file is distributed under the terms of the

803

! GNU General Public License: see COPYING in the top directory

804

! or http://www.gnu.org/copyleft/gpl.txt .

805

! ---

806

!

807

! This subroutine extracts from an object of type State_type specific values

808

! as requested by the caller. The first variable, which is mandatory, is the

809

! state object being interrogated. Then follows a keyword-list of desired

810

! values (temperature, pressure, etc), all of which are optional, that are

811

! returned upon request.

812

!

813

!*****************************************************************************

814

! modules used

815

816

!*****************************************************************************

817

! No implicits please!

818

819

implicit none

820

821

!*****************************************************************************

822

! arguments

823

824

type ( state_type ), intent (INOUT) :: state ! the state being interrogated

825

826

integer, intent (OUT), optional :: nAtoms

827

integer, intent (OUT), optional :: nDegreesOfFreedom

828

integer, intent (OUT), optional :: nThermostats

829

830

real (dp), intent (OUT), optional :: H0

831

832

real (dp), dimension (:,:), intent (OUT), optional :: position

833

real (dp), dimension (:,:), intent (OUT), optional :: force

834

real (dp), dimension (:,:), intent (OUT), optional :: momentum

835

real (dp), dimension (:), intent (OUT), optional :: mass

836

837

type ( thermostat_type ), dimension (:), intent (OUT), optional :: thermostat

838

type ( barostat_type ), intent (OUT), optional :: barostat

839

840

type ( energy_type ), intent (OUT), optional :: energy

841

842

type ( lattice_type ), intent (OUT), optional :: Cell

843

type ( lattice_type ), intent (OUT), optional :: reciprocalCell

844

845

type ( stress_type ), intent (OUT), optional :: stress

846

847

type ( representation_type ), intent (IN), optional :: representation

848

849

!*****************************************************************************

850

! variables

851

852

integer :: i

853

854

!*****************************************************************************

855

! start of subroutine

856

857

if ( debug ) write(*,*) 'Entering get()'

858

859

if ( present( nAtoms ) ) then

860

861

nAtoms = state % nAtoms

862

863

end if

864

865

if ( present( nDegreesOfFreedom ) ) then

866

867

nDegreesOfFreedom = state % nDegreesOfFreedom

868

869

end if

870

871

if ( present( nThermostats ) ) then

872

873

nThermostats = state % nThermostats

874

875

end if

876

877

if ( present( position ) ) then

878

879

if ( present( representation ) ) then

880

881

if ( representation % position ) then ! coordinates are needed in cartesian rep.

882

883

call transformCoordinates2Cartesian( state )

884

885

else ! coordinates are needed in lattice rep.

886

887

call transformCoordinates2Lattice( state )

888

889

end if

890

891

else ! default; make sure they are in cartesian

892

893

call transformCoordinates2Cartesian( state )

894

895

end if

896

897

position = state % position

898

899

end if

900

901

if ( present( force ) ) then

902

903

if ( present( representation ) ) then

904

905

if ( representation % force ) then ! forces are needed in cartesian rep.

906

907

call transformForces2Cartesian( state )

908

909

else ! forces are needed in lattice rep.

910

911

call transformForces2Lattice( state )

912

913

end if

914

915

else ! default, we return them in cartesian, so make sure they are in that rep.

916

917

call transformForces2Cartesian( state )

918

919

end if

920

921

force = state % force ! copy forces

922

923

end if

924

925

if ( present( momentum ) ) then

926

927

if ( present( representation ) ) then

928

929

if ( representation % momentum ) then ! momenta are needed in cartesian rep.

930

931

call transformMomenta2Cartesian( state )

932

933

else ! momenta are needed in lattice rep.

934

935

call transformMomenta2Lattice( state )

936

937

end if

938

939

else ! default case, we return them in cartesian

940

941

call transformMomenta2Cartesian( state )

942

943

end if

944

945

momentum = state % momentum ! copy momenta

946

947

end if

948

949

if ( present( mass ) ) then

950

951

mass = state % mass

952

953

end if

954

955

if ( present( thermostat ) ) then

956

957

do i = 1, state % nThermostats

958

959

thermostat(i) = state % thermostat(i)

960

961

end do

962

963

end if

964

965

if ( present( barostat ) ) then

966

967

barostat = state % barostat

968

969

end if

970

971

if ( present( energy ) ) then

972

973

energy = state % energy

974

975

end if

976

977

if ( present( H0 ) ) then

978

979

H0 = state % H0 % E

980

981

end if

982

983

if ( present( Cell ) ) then

984

985

Cell = state % Cell

986

987

end if

988

989

if ( present( reciprocalCell ) ) then

990

991

reciprocalCell = state % reciprocalCell

992

993

end if

994

995

if ( present( stress ) ) then

996

997

stress = state % stress

998

999

end if

1000

1001

if ( debug ) write(*,*) 'Exiting get()'

1002

1003

end subroutine get

1004

1005

!*****************************************************************************

1006

subroutine getConfiguration( state, nAtoms, a, b, c, position, lattice )

1007

!*****************************************************************************

1008

!

1009

! Project: MolecularDynamics

1010

!

1011

! Module: State

1012

!

1013

! Created on: Wed 13 Jun 18:04:27 2018 by Eduardo R. Hernandez

1014

!

1015

! Last Modified: Wed 18 Jul 17:12:54 2018

1016

!

1017

! ---

1018

1019

! This file is distributed under the terms of the

1020

! GNU General Public License: see COPYING in the top directory

1021

! or http://www.gnu.org/copyleft/gpl.txt .

1022

! ---

1023

!

1024

! This subroutine is designed as a means of communicating with the client

1025

! program; the latter will need to be told the lattice vectors and atom positions

1026

! in units that it understands. Lattice vectors and ionic positions are

1027

! privately held by state, and are thus directly not accessible to the

1028

! client program. This subroutine provides, given an instance of state, the

1029

! lattice vectors a, b, c and the ionic positions, in the chosen length

1030

! units of the client program. Optionally, atomic positions can be returned

1031

! in lattice coordinates, by appropriately setting the optional flag Lattice

1032

! as .true.

1033

!

1034

!*****************************************************************************

1035

! modules used

1036

1037

!*****************************************************************************

1038

! No implicits please!

1039

1040

implicit none

1041

1042

!*****************************************************************************

1043

! variables

1044

1045

type ( state_type ), intent (INOUT) :: state

1046

1047

integer , intent (OUT) :: nAtoms

1048

1049

real (dp), dimension (3), intent (OUT) :: a

1050

real (dp), dimension (3), intent (OUT) :: b ! the three lattice vectors

1051

real (dp), dimension (3), intent (OUT) :: c

1052

1053

real (dp), dimension ( 3, state % nAtoms ), intent (OUT) :: position

1054

1055

logical, optional, intent (IN) :: lattice

1056

1057

!*****************************************************************************

1058

! variables

1059

1060

!*****************************************************************************

1061

! start of subroutine

1062

1063

if ( debug ) write(*,*) 'Entering getConfiguration()'

1064

1065

! convert lattice vectors to client program length units

1066

1067

a = lengthConversionFactor * state % Cell % a

1068

b = lengthConversionFactor * state % Cell % b

1069

c = lengthConversionFactor * state % Cell % c

1070

1071

! the number of atoms

1072

1073

nAtoms = state % nAtoms

1074

1075

! now atomic coordinates

1076

1077

if ( present( lattice ) ) then ! check if client wants coordinates in lattice rep.

1078

1079

if ( lattice ) then

1080

1081

call transformCoordinates2Lattice( state ) ! make sure they are in lattice rep.

1082

1083

position = state % position

1084

1085

else

1086

1087

call transformCoordinates2Cartesian( state ) ! make sure they are in cartesian rep.

1088

1089

position = lengthConversionFactor * state % position

1090

1091

end if

1092

1093

else ! the default case is that coordinates are returned in Cartesian rep.

1094

1095

call transformCoordinates2Cartesian( state ) ! make sure they are in cartesian rep.

1096

1097

position = lengthConversionFactor * state % position

1098

1099

end if

1100

1101

if ( debug ) write(*,*) 'Exiting getConfiguration()'

1102

1103

end subroutine getConfiguration

1104

1105

!*****************************************************************************

1106

subroutine initialiseState( state, &

1107

a, b, c, &

1108

mass, position, temperature, velocity, lattice, &

1109

periodic )

1110

!*****************************************************************************

1111

!

1112

! Project: MolecularDynamics

1113

!

1114

! Module: State

1115

!

1116

! Created on: Fri 2 Nov 15:55:03 2018 by Eduardo R. Hernandez

1117

!

1118

! Last Modified: Fri 2 Nov 16:17:40 2018

1119

!

1120

! ---

1121

1122

! This file is distributed under the terms of the

1123

! GNU General Public License: see COPYING in the top directory

1124

! or http://www.gnu.org/copyleft/gpl.txt .

1125

! ---

1126

!

1127

! This subroutine is called to "fill" an empty but allocated instance of State

1128

! with state information (atom coordinates, cell parameters, etc)

1129

!

1130

! variables:

1131

!

1132

! type ( state_type ), intent (INOUT) :: state

1133

! this is the state that is going to be allocated and filled with the data

1134

! passed next.

1135

!

1136

! integer, intent (IN) :: nAtoms

1137

! the number of atoms in the system; used to dimension several arrays below

1138

!

1139

! real (dp), dimension (3), intent (IN) :: a, b, c

1140

! the lattice vectors of the simulation cell

1141

!

1142

! real (dp), dimension (:), intent (IN) :: mass

1143

! an array of dimension (nAtoms) containing the masses of all atoms in

1144

! Daltons (1/12 of the mass of C12)

1145

!

1146

! real (dp), dimension (:,:), intent (IN) :: position

1147

! an array of dimension (3,nAtoms) containing the coordinates of atoms;

1148

! if the optional flag 'lattice' (see below) is passed as .true., these

1149

! are assumed to be passed in lattice representation (x,y,z in [0,10])

1150

!

1151

! real (dp), intent(IN), optional :: temperature

1152

! the temperature of the simulation; the user can provide either the simulation

1153

! temperature (in which case random initial momenta are generated according to

1154

! the Maxwel-Boltzmann distribution at that temperature), or else a set of

1155

! initial velocities (see next entry)

1156

!

1157

! real (dp), dimension (:,:), intent (IN), optional :: velocity

1158

! optionally, the client program can pass initial velocities to be used

1159

! in the construction of the momenta that are used internally; if absent,

1160

! random momental will be sampled from the Maxwell-Boltzmann distribution

1161

! corresponding to the temperature of the simulation. If provided, velocities

1162

! must be in units compatible with the position; also note that if the flag

1163

! 'lattice' is set to .true., the velocities must also be in this representation,

1164

! to be compatible with the atomic coordinates (see position, above)

1165

!

1166

! logical, intent (IN), optional :: lattice

1167

! if .true., this variable indicates that coordinates (position) and velocity are

1168

! given in lattice representation. If set to .false. (or absent), it is assumed

1169

! that coordinates (and velocities) are in Cartesian representation in units

1170

! of the client program.

1171

!

1172

! logical, intent (IN), optional :: periodic

1173

! if .true. (default) this indicates that periodic boundary conditions will be

1174

! imposed; if .false. then a cluster geometry is assumed. The only relevance of this

1175

! here is that in the latter case the momenta are generated such that not only the

1176

! total linear momentum is zero, but also the total angular momentum is zero, to

1177

! avoid whole-cluster rotation, which makes analysis more cumbersome.

1178

!

1179

!*****************************************************************************

1180

! modules used

1181

1182

!*****************************************************************************

1183

! No implicits please!

1184

1185

implicit none

1186

1187

!*****************************************************************************

1188

! arguments

1189

1190

type ( state_type ), intent (INOUT) :: state

1191

1192

real (dp), dimension (3), intent (IN) :: a ! lattice vectors in

1193

real (dp), dimension (3), intent (IN) :: b ! length units of the

1194

real (dp), dimension (3), intent (IN) :: c ! client program

1195

1196

real (dp), dimension (:), intent (IN) :: mass

1197

real (dp), dimension (:,:), intent (IN) :: position

1198

1199

real (dp), intent (IN), optional :: temperature

1200

real (dp), dimension (:,:), intent (IN), optional :: velocity

1201

1202

logical, intent (IN), optional :: lattice

1203

logical, intent (IN), optional :: periodic

1204

1205

!*****************************************************************************

1206

! local variables

1207

1208

integer :: i

1209

1210

real (dp), dimension (3,3) :: LatticeVectors

1211

1212

real (dp), allocatable, dimension (:,:) :: momentum

1213

1214

!*****************************************************************************

1215

! start subroutine

1216

1217

if ( debug ) write(*,*) 'Entering initialiseState()'

1218

1219

! set the cell

1220

1221

LatticeVectors(1:3,1) = lengthConversionFactor * a

1222

LatticeVectors(1:3,2) = lengthConversionFactor * b

1223

LatticeVectors(1:3,3) = lengthConversionFactor * c

1224

1225

call updateCell( state % Cell, state % reciprocalCell, LatticeVectors = LatticeVectors )

1226

1227

! store masses internally

1228

1229

state % mass = massConversionFactor * mass

1230

1231

! now the atomic positions

1232

1233

if ( present( lattice ) ) then

1234

1235

if ( lattice ) then ! coordinates are in lattice representation

1236

1237

state % position = position

1238

state % representation % position = .false.

1239

1240

else

1241

1242

state % position = lengthConversionFactor * position

1243

state % representation % position = .true.

1244

1245

end if

1246

1247

else

1248

1249

state % position = lengthConversionFactor * position

1250

state % representation % position = .true.

1251

1252

end if

1253

1254

if ( present( periodic ) ) then

1255

1256

PBC = periodic

1257

1258

else

1259

1260

PBC = .true.

1261

1262

end if

1263

1264

if ( present( temperature ) .and. ( .not. present( velocity ) ) ) then

1265

1266

call generateMomenta( state, temperature )

1267

1268

else if ( present( velocity ) ) then

1269

1270

if ( present( lattice ) ) then

1271

1272

if ( lattice ) then ! velocities are in lattice units

1273

1274

allocate( momentum( 3, state % nAtoms ) )

1275

1276

do i = 1, state % nAtoms

1277

momentum(1:3,i) = state % mass(i) * velocity(1:3,i)

1278

end do

1279

1280

state % momentum = momentum

1281

state % representation % momentum = .false.

1282

1283

deallocate( momentum )

1284

1285

else

1286

1287

do i = 1, state % nAtoms

1288

state % momentum(1:3,i) = &

1289

velConversionFactor * state % mass(i) * velocity(1:3,i)

1290

end do

1291

1292

state % representation % momentum = .true. ! cartesian rep. for momenta

1293

1294

end if

1295

1296

else ! lattice not present, so Cartesian

1297

1298

do i = 1, state % nAtoms

1299

state % momentum(1:3,i) = &

1300

state % mass(i) * velocity(1:3,i)

1301

!velConversionFactor * state % mass(i) * velocity(1:3,i)

1302

end do

1303

1304

state % representation % momentum = .true. ! cartesian rep. for momenta

1305

1306

end if

1307

1308

! write(*,*) 'velConversionFactor = ', velConversionFactor

1309

! write(*,*) 'vel in createState ', velocity(1:3,1)

1310

! write(*,*) 'createState ', state % momentum(1:3,1)

1311

1312

end if

1313

1314

if ( debug ) write(*,*) 'Exiting initialiseState'

1315

1316

end subroutine initialiseState

1317

1318

!*****************************************************************************

1319

subroutine inquire( state, temperature, pressure, stress, &

1320

kineticEnergy, potentialEnergy, thermostatEnergy, &

1321

barostatEnergy, constantOfMotion )

1322

!*****************************************************************************

1323

!

1324

! Project: MolecularDynamics

1325

!

1326

! Module: State

1327

!

1328

! Created on: Fri 27 Apr 12:15:38 2018 by Eduardo R. Hernandez

1329

!

1330

! Last Modified: Fri 29 Jun 16:43:24 2018

1331

!

1332

! ---

1333

1334

! This file is distributed under the terms of the

1335

! GNU General Public License: see COPYING in the top directory

1336

! or http://www.gnu.org/copyleft/gpl.txt .

1337

! ---

1338

!

1339

! This subroutine inquires state to obtain properties of the current system

1340

! configuration, such as the temperature, pressure, or different energy

1341

! components. Properties are output in the chosen units of the client program.

1342

!

1343

!*****************************************************************************

1344

! modules used

1345

1346

!*****************************************************************************

1347

! No implicits please!

1348

1349

implicit none

1350

1351

!*****************************************************************************

1352

! argument

1353

1354

type ( state_type ), intent (INOUT) :: state

1355

1356

real (dp), intent (OUT), optional :: temperature

1357

real (dp), intent (OUT), optional :: pressure

1358

real (dp), intent (OUT), dimension (3,3), optional :: stress

1359

real (dp), intent (OUT), optional :: kineticEnergy

1360

real (dp), intent (OUT), optional :: potentialEnergy

1361

real (dp), intent (OUT), optional :: thermostatEnergy

1362

real (dp), intent (OUT), optional :: barostatEnergy

1363

real (dp), intent (OUT), optional :: constantOfMotion

1364

1365

!*****************************************************************************

1366

! variables

1367

1368

integer :: i

1369

1370

real(dp) :: const

1371

real(dp) :: factor

1372

1373

!*****************************************************************************

1374

! start of subroutine

1375

1376

if ( debug ) write(*,*) 'Entering inquire()'

1377

1378

! make sure that the kinetic energy is consistent with current momenta

1379

1380

call calculateKineticEnergy( state )

1381

1382

if ( present( temperature ) ) then

1383

1384

temperature = two * state % energy % kinetic / ( state % nDegreesOfFreedom * boltzmann_k )

1385

1386

end if

1387

1388

if ( present ( pressure ) ) then

1389

1390

end if

1391

1392

if ( present( stress ) ) then

1393

1394

stress = pressureConversionFactor * state % stress % cartesian

1395

1396

end if

1397

1398

if ( present( kineticEnergy ) ) then

1399

1400

kineticEnergy = energyConversionFactor * state % energy % kinetic

1401

1402

end if

1403

1404

if ( present( potentialEnergy ) ) then

1405

1406

potentialEnergy = energyConversionFactor * state % energy % potential

1407

1408

end if

1409

1410

if ( present( thermostatEnergy ) ) then

1411

1412

thermostatEnergy = energyConversionFactor * state % energy % thermostat

1413

1414

end if

1415

1416

if ( present( barostatEnergy ) ) then

1417

1418

barostatEnergy = energyConversionFactor * state % energy % barostat

1419

1420

end if

1421

1422

if ( present( constantOfMotion ) ) then

1423

1424

constantOfMotion = energyConversionFactor * state % energy % constantOfMotion

1425

1426

end if

1427

1428

if ( debug ) write(*,*) 'Exiting inquire()'

1429

1430

end subroutine inquire

1431

1432

!*****************************************************************************

1433

subroutine set( state, H, G, &

1434

position, force, momentum, mass, &

1435

thermostat, barostat, energy, &

1436

Cell, reciprocalCell, representation )

1437

!*****************************************************************************

1438

!

1439

! Project: MolecularDynamics

1440

!

1441

! Module: State

1442

!

1443

! Created on: Fri 8 Jun 14:32:53 2018 by Eduardo R. Hernandez

1444

!

1445

! Last Modified: Tue 6 Nov 09:38:52 2018

1446

!

1447

! ---

1448

1449

! This file is distributed under the terms of the

1450

! GNU General Public License: see COPYING in the top directory

1451

! or http://www.gnu.org/copyleft/gpl.txt .

1452

! ---

1453

!

1454

! This subroutine has the purpose of setting the values of variables contained

1455

! in an instance of State to new numerical values. An object of type State

1456

! must be passed as first argument; the remaining arguments are optional and

1457

! identified by keyword, so that a single variable, a subset or the whole set

1458

! can be specified in one go. The updated State object is then returned.

1459

!

1460

!*****************************************************************************

1461

! modules used

1462

1463

!*****************************************************************************

1464

! No implicits please!

1465

1466

implicit none

1467

1468

!*****************************************************************************

1469

! arguments

1470

1471

type ( state_type ), intent (INOUT) :: state ! the state being modified

1472

1473

real (dp), dimension (3,3), intent (IN), optional :: H ! the lattice vectors (one per row)

1474

real (dp), dimension (3,3), intent (IN), optional :: G ! the metric tensor

1475

1476

real (dp), dimension ( 3, state % nAtoms ), intent (IN), optional :: position

1477

real (dp), dimension ( 3, state % nAtoms ), intent (IN), optional :: force

1478

real (dp), dimension ( 3, state % nAtoms ), intent (IN), optional :: momentum

1479

real (dp), dimension ( state % nAtoms ), intent (IN), optional :: mass

1480

1481

type (thermostat_type), dimension ( state % nThermostats ), intent (IN), optional :: thermostat

1482

type (barostat_type), intent (IN), optional :: barostat

1483

1484

type (energy_type), intent (IN), optional :: energy

1485

1486

type (lattice_type), intent (IN), optional :: Cell

1487

type (lattice_type), intent (IN), optional :: reciprocalCell

1488

1489

type (representation_type), intent (IN), optional :: representation

1490

1491

!*****************************************************************************

1492

! local variables

1493

1494

integer :: n

1495

1496

!*****************************************************************************

1497

! start subroutine

1498

1499

if ( debug ) write(*,*) 'Entering set()'

1500

1501

! the user must provide either cell vectors or metric tensor, but not both,

1502

! hence the if-then-else

1503

1504

if ( present( H ) ) then

1505

1506

call updateCell( state % Cell, state % reciprocalCell, LatticeVectors = H )

1507

1508

else if ( present( G ) ) then

1509

1510

call updateCell( state % Cell, state % reciprocalCell, MetricTensor = G )

1511

1512

end if

1513

1514

if ( present( position ) ) then

1515

1516

if ( present( representation ) ) then

1517

1518

if ( representation % position ) then

1519

1520

state % representation % position = .true. ! coordinates are input in cartesian rep.

1521

1522

else

1523

1524

state % representation % position = .false. ! coordinates are input in lattice rep.

1525

1526

end if

1527

1528

else ! default case; it is assumed coordinates provided in cartesian

1529

1530

state % representation % position = .true.

1531

1532

end if

1533

1534

state % position = position

1535

1536

end if

1537

1538

if ( present( force ) ) then

1539

1540

if ( present( representation ) ) then

1541

1542

if ( representation % force ) then

1543

1544

state % representation % force = .true. ! forces are input in cartesian rep.

1545

1546

else

1547

1548

state % representation % force = .false. ! they are in lattice rep.

1549

1550

end if

1551

1552

else

1553

1554

state % representation % force = .true. ! default, they are in cartesian

1555

1556

end if

1557

1558

state % force = force

1559

1560

end if

1561

1562

if ( present( momentum ) ) then

1563

1564

if ( present( representation ) ) then

1565

1566

if ( representation % momentum ) then

1567

1568

state % representation % momentum = .true. ! momenta are input in cartesian rep.

1569

1570

else

1571

1572

state % representation % momentum = .false. ! nope, they are in lattice rep.

1573

1574

end if

1575

1576

else

1577

1578

state % representation % momentum = .true. ! default case, they are in cartesian

1579

1580

end if

1581

1582

state % momentum = momentum

1583

1584

end if

1585

1586

if ( present( mass ) ) state % mass = mass

1587

1588

if ( present( thermostat ) ) then

1589

1590

do n = 1, state % nThermostats

1591

state % thermostat(n) = thermostat(n)

1592

end do

1593

1594

end if

1595

1596

if ( present( barostat ) ) then

1597

1598

state % barostat = barostat

1599

1600

end if

1601

1602

if ( present( energy ) ) then

1603

1604

state % energy = energy

1605

1606

if ( .not. state % H0 % set ) then

1607

1608

state % H0 % E = state % energy % kinetic + state % energy % potential + &

1609

state % energy % barostat + state % energy % thermostat

1610

1611

state % H0 % set = .true.

1612

1613

end if

1614

1615

end if

1616

1617

if ( present( Cell ) ) then

1618

1619

state % Cell = Cell

1620

1621

end if

1622

1623

if ( present( reciprocalCell ) ) then

1624

1625

state % reciprocalCell = reciprocalCell

1626

1627

end if

1628

1629

if ( debug ) write(*,*) 'Exiting set()'

1630

1631

end subroutine set

1632

1633

!*****************************************************************************

1634

subroutine setForce( state, potential_energy, force, stress, lattice )

1635

!*****************************************************************************

1636

!

1637

! Project: MolecularDynamics

1638

!

1639

! Module: State

1640

!

1641

! Created on: Thu 14 Jun 09:33:11 2018 by Eduardo R. Hernandez

1642

!

1643

! Last Modified: Thu 26 Jul 17:16:33 2018

1644

!

1645

! ---

1646

1647

! This file is distributed under the terms of the

1648

! GNU General Public License: see COPYING in the top directory

1649

! or http://www.gnu.org/copyleft/gpl.txt .

1650

! ---

1651

!

1652

! This subroutine receives the forces, and optionally the stress tensor, as

1653

! calculated by the client program, and stores this information in State.

1654

! It is assumed that force and stress information are input in Cartesian

1655

! representation. The client program will normally pass the forces in

1656

! Cartesian representation, but if the optional flag lattice is set tot .true.,

1657

! then it is understood that forces are passed in lattice representation.

1658

!

1659

!*****************************************************************************

1660

! modules used

1661

1662

!*****************************************************************************

1663

! No implicits please!

1664

1665

implicit none

1666

1667

!*****************************************************************************

1668

! variables

1669

1670

type ( state_type ), intent (INOUT) :: state

1671

1672

real (dp), intent (IN) :: potential_energy

1673

real (dp), dimension ( 3, state % nAtoms ), intent (IN) :: force

1674

real (dp), dimension ( 3, 3 ), intent (IN), optional :: stress

1675

1676

logical, intent (IN), optional :: lattice

1677

1678

!*****************************************************************************

1679

! local variables

1680

1681

real (dp), dimension ( 3, 3 ) :: matrix_a

1682

real (dp), dimension ( 3, 3 ) :: matrix_b

1683

1684

!*****************************************************************************

1685

! start of subroutine

1686

1687

if ( debug ) write(*,*) 'Entering setForce()'

1688

1689

state % energy % potential = potential_energy / energyConversionFactor

1690

1691

if ( present( lattice ) ) then

1692

1693

if ( lattice ) then ! forces are given in lattice rep;

1694

! TODO: make sure conversion factor

1695

! is right in this case

1696

state % representation % force = .false.

1697

1698

else

1699

1700

state % representation % force = .true. ! cartesian

1701

1702

end if

1703

1704

else

1705

1706

state % representation % force = .true.

1707

1708

end if

1709

1710

state % force = force / forceConversionFactor

1711

1712

! if stress has been passed, convert it too

1713

1714

if ( present( stress ) ) then

1715

1716

! state % stress % cartesian = stress / pressureConversionFactor

1717

state % stress % cartesian = stress

1718

1719

matrix_a(1:3,1) = state % reciprocalCell % a

1720

matrix_a(1:3,2) = state % reciprocalCell % b

1721

matrix_a(1:3,3) = state % reciprocalCell % c

1722

1723

matrix_b = matmul( state % stress % cartesian, matrix_a )

1724

1725

state % stress % lattice = half * matmul( transpose( matrix_a ), matrix_b )

1726

1727

else

1728

1729

state % stress % cartesian = zero

1730

state % stress % lattice = zero

1731

1732

end if

1733

1734

if ( debug ) write(*,*) 'Exiting setForce()'

1735

1736

end subroutine setForce

1737

1738

!*****************************************************************************

1739

subroutine setUnits( length, energy, pressure )

1740

!*****************************************************************************

1741

!

1742

! Project: MolecularDynamics

1743

!

1744

! Module: Structure

1745

!

1746

! Created on: Thu 14 Jun 11:22:27 2018 by Eduardo R. Hernandez

1747

!

1748

! Last Modified: Mon 25 Jun 18:00:04 2018

1749

!

1750

! ---

1751

1752

! This file is distributed under the terms of the

1753

! GNU General Public License: see COPYING in the top directory

1754

! or http://www.gnu.org/copyleft/gpl.txt .

1755

! ---

1756

!

1757

! This subroutine is called once at the begining, and its purpose is to

1758

! set a number of unit conversion factors in order to manage communication

1759

! to/from the client program. The idea is that both the State type and

1760

! the routines in the MolecularDynamics module work with the same unit system,

1761

! namely atomic units (Bohr, Hartree, etc); the client program will probably

1762

! use other units (typical choices are Angstrom, eV, etc), so we must

1763

! cater for the appropriate conversion factors.

1764

!

1765

!!!!!!!!!!!!!!! MANDATORY ARGUMENTS !!!!!!!!!!!!!!!

1766

!

1767

! length (INT): 0 specifies client length unit is bohr (atomic unit)

1768

! 1 client length unit is Angstrom

1769

!

1770

! energy (INT): 0 specifies client energy unit is hartree (atomic unit)

1771

! 1 specifies client energy unit is rydberg

1772

! 2 specifies client energy unit is eV

1773

!

1774

!!!!!!!!!!!!!!! OPTIONAL ARGUMENTS !!!!!!!!!!!!!!!

1775

!

1776

! pressure (INT): if this is a constant-pressure simulation, or if the

1777

! user is going to ask for values of the instantaneous

1778

! internal pressure, specify here the units pressure is to be

1779

! input/output in

1780

! 0: specifies client pressure/stress is in kB

1781

! 1: specifies client pressure/stress is in GPa

1782

!

1783

!*****************************************************************************

1784

! modules used

1785

1786

!*****************************************************************************

1787

! No implicits please!

1788

1789

implicit none

1790

1791

!*****************************************************************************

1792

! variables

1793

1794

integer, intent (IN) :: length

1795

integer, intent (IN) :: energy

1796

integer, intent (IN), optional :: pressure

1797

1798

!*****************************************************************************

1799

! begin subroutine

1800

1801

if ( debug ) write(*,*) 'Entering setUnits()'

1802

1803

! set the length unit conversion factor

1804

1805

select case ( length )

1806

1807

case ( 0 ) ! client length units are in bohr, good choice!

1808

1809

lengthConversionFactor = one

1810

1811

case ( 1 ) ! uh uh, user chose Angstrom

1812

1813

lengthConversionFactor = Ang

1814

1815

end select

1816

1817

! set the energy unit conversion factor

1818

1819

! set the energy unit conversion factor

1820

1821

select case ( energy )

1822

1823

case ( 0 ) ! client is clever and is using hartree units

1824

1825

energyConversionFactor = one

1826

1827

case ( 1 ) ! client is using rydbergs; well, could be worse

1828

1829

energyConversionFactor = half

1830

1831

case ( 2 ) ! client is one of those renegades using eV

1832

1833

energyConversionFactor = 1.0_dp / 27.2116_dp

1834

1835

end select

1836

1837

! with length and energy units, we can set the conversion factors

1838

! for other magnitudes

1839

1840

forceConversionFactor = energyConversionFactor / lengthConversionFactor

1841

massConversionFactor = 1822.8885302342505_dp ! convert Dalton to electron masses

1842

timeConversionFactor = sqrt( au_to_joules / amu ) / adu

1843

velConversionFactor = lengthConversionFactor / timeConversionFactor

1844

1845

! finally, pressure units, if required

1846

1847

if ( present( pressure ) ) then

1848

1849

select case( pressure )

1850

1851

case ( 0 ) ! client units are kB

1852

1853

pressureConversionFactor = 3.398923e-4_dp

1854

1855

case ( 1 ) ! client units are GPa

1856

1857

pressureConversionFactor = 3.398923e-5_dp

1858

1859

end select

1860

1861

else ! use default as GPa

1862

1863

pressureConversionFactor = 3.398923e-5_dp

1864

1865

end if

1866

1867

if ( debug ) write(*,*) 'Exiting setUnits()'

1868

1869

end subroutine setUnits

1870

1871

!*****************************************************************************

1872

subroutine transformCoordinates2Lattice( state )

1873

!*****************************************************************************

1874

!

1875

! Project: MolecularDynamics

1876

!

1877

! Module: State

1878

!

1879

! Created on: Thu 8 Mar 11:39:29 2018 by Eduardo R. Hernandez

1880

!

1881

! Last Modified: Wed 18 Jul 14:29:25 2018

1882

!

1883

! ---

1884

1885

! This file is distributed under the terms of the

1886

! GNU General Public License: see COPYING in the top directory

1887

! or http://www.gnu.org/copyleft/gpl.txt .

1888

! ---

1889

!

1890

! This subroutine transforms the coordinates from Cartesian to Lattice

1891

! representation. Internally, type State contains the coordinates in Cartesian

1892

! representation, but they may be required by MD routines or client program

1893

! in Lattice representation. This soubroutine takes care of the transformation.

1894

! The internal representation is not changed, only the copied array is.

1895

!

1896

!*****************************************************************************

1897

! modules used

1898

1899

!*****************************************************************************

1900

! No implicits please!

1901

1902

implicit none

1903

1904

!*****************************************************************************

1905

! interface variables

1906

1907

type ( state_type ), intent (INOUT) :: state

1908

1909

!*****************************************************************************

1910

! internal variables

1911

1912

integer :: i

1913

1914

real (dp), dimension (3) :: coordinates

1915

1916

!*****************************************************************************

1917

! begin subroutine

1918

1919

if ( debug ) write(*,*) 'Entering transformCoordinates2Lattice()'

1920

1921

if ( state % representation % position ) then ! coordinates are in cartesian rep., so work to do

1922

1923

do i=1, state % nAtoms

1924

1925

coordinates = state % position(1:3,i)

1926

1927

state % position(1,i) = dot_product( state % reciprocalCell % a, &

1928

coordinates )

1929

state % position(2,i) = dot_product( state % reciprocalCell % b, &

1930

coordinates )

1931

state % position(3,i) = dot_product( state % reciprocalCell % c, &

1932

coordinates )

1933

1934

end do

1935

1936

state % representation % position = .false. ! now they are in lattice rep.

1937

1938

else ! coordinates are in lattice already, so there is nothing to do

1939

1940

return

1941

1942

end if

1943

1944

if ( debug ) write(*,*) 'Exiting transformCoordinates2Lattice()'

1945

1946

end subroutine transformCoordinates2Lattice

1947

1948

!*****************************************************************************

1949

subroutine transformCoordinates2Cartesian( state )

1950

!*****************************************************************************

1951

!

1952

! Project:

1953

!

1954

! Created on: Tue 19 Jun 13:44:20 2018 by Eduardo R. Hernandez

1955

!

1956

! Last Modified: Wed 18 Jul 14:30:26 2018

1957

!

1958

! ---

1959

1960

! This file is distributed under the terms of the

1961

! GNU General Public License: see COPYING in the top directory

1962

! or http://www.gnu.org/copyleft/gpl.txt .

1963

! ---

1964

!

1965

! This subroutine transforms the coordinates from Lattice to Cartesian representation,

1966

! when the client program provides them in that representation.

1967

!

1968

!*****************************************************************************

1969

! modules used

1970

1971

!*****************************************************************************

1972

! No implicits please!

1973

1974

implicit none

1975

1976

!*****************************************************************************

1977

! interface variables

1978

1979

type ( state_type ), intent (INOUT) :: state

1980

1981

!*****************************************************************************

1982

! internal variables

1983

1984

integer :: i

1985

1986

real (dp), dimension (3) :: coordinates, r, s, t

1987

1988

!*****************************************************************************

1989

! begin subroutine

1990

1991

if ( debug ) write(*,*) 'Entering transformCoordinates2Cartesian()'

1992

1993

if ( state % representation % position ) then ! if coordinates are already cartesian, do nothing

1994

1995

return

1996

1997

else ! coordinates are in lattice rep., so transform them

1998

1999

r = (/ state % cell % a(1), state % cell % b(1), state % cell % c(1) /)

2000

s = (/ state % cell % a(2), state % cell % b(2), state % cell % c(2) /)

2001

t = (/ state % cell % a(3), state % cell % b(3), state % cell % c(3) /)

2002

2003

do i=1, state % nAtoms

2004

2005

coordinates = state % position(1:3,i)

2006

2007

state % position(1,i) = dot_product( r, coordinates )

2008

state % position(2,i) = dot_product( s, coordinates )

2009

state % position(3,i) = dot_product( t, coordinates )

2010

2011

end do

2012

2013

state % representation % position = .true.

2014

2015

end if

2016

2017

if ( debug ) write(*,*) 'Exiting transformCoordinates2Cartesian()'

2018

2019

end subroutine transformCoordinates2Cartesian

2020

2021

!*****************************************************************************

2022

subroutine transformForces2Lattice( state )

2023

!*****************************************************************************

2024

!

2025

! Project: MolecularDynamics

2026

!

2027

! Module: State

2028

!

2029

! Created on: Thu 8 Mar 12:10:34 2018 by Eduardo R. Hernandez

2030

!

2031

! Last Modified: Wed 18 Jul 16:43:32 2018

2032

!

2033

! ---

2034

2035

! This file is distributed under the terms of the

2036

! GNU General Public License: see COPYING in the top directory

2037

! or http://www.gnu.org/copyleft/gpl.txt .

2038

! ---

2039

!

2040

! This subroutine transforms the forces from Cartesian representation (the

2041

! way they are stored internally) to Lattice, when they are requested in this

2042

! representation by the program calling get() (usually this call will be from

2043

! a subroutine in MolecularDynamics module).

2044

!

2045

!*****************************************************************************

2046

! modules used

2047

2048

!*****************************************************************************

2049

! No implicits please!

2050

2051

implicit none

2052

2053

!*****************************************************************************

2054

! interface variables

2055

2056

type ( state_type ), intent (INOUT) :: state

2057

2058

!*****************************************************************************

2059

! internal variables

2060

2061

integer :: i

2062

2063

real(dp), dimension (3) :: vector

2064

2065

!*****************************************************************************

2066

! begin subroutine

2067

2068

if ( debug ) write(*,*) 'Entering transformForces2Lattice()'

2069

2070

if ( state % representation % force ) then ! forces are in cartesian rep., so transform

2071

2072

do i=1, state % nAtoms

2073

2074

vector = state % force(1:3,i)

2075

2076

state % force(1,i) = dot_product( state % cell % a, vector )

2077

state % force(2,i) = dot_product( state % cell % b, vector )

2078

state % force(3,i) = dot_product( state % cell % c, vector )

2079

2080

end do

2081

2082

state % representation % force = .false. ! now they are in lattice rep.

2083

2084

else ! forces already in lattice rep., so do nothing

2085

2086

return

2087

2088

end if

2089

2090

if ( debug ) write(*,*) 'Exiting transformForces2Lattice()'

2091

2092

end subroutine transformForces2Lattice

2093

2094

!*****************************************************************************

2095

subroutine transformForces2Cartesian( state )

2096

!*****************************************************************************

2097

!

2098

! Project: MolecularDynamics

2099

!

2100

! Created on: Thu 8 Mar 12:10:34 2018 by Eduardo R. Hernandez

2101

!

2102

! Last Modified: Wed 18 Jul 16:42:30 2018

2103

!

2104

! ---

2105

2106

! This file is distributed under the terms of the

2107

! GNU General Public License: see COPYING in the top directory

2108

! or http://www.gnu.org/copyleft/gpl.txt .

2109

! ---

2110

!

2111

! This is a subroutine of module configuration_module; it transforms the

2112

! forces from Lattice representation, if the client provides them in that way

2113

! to Cartesian, the way in which they are stored by State type.

2114

!

2115

!*****************************************************************************

2116

! modules used

2117

2118

!*****************************************************************************

2119

! No implicits please!

2120

2121

implicit none

2122

2123

!*****************************************************************************

2124

! interface variables

2125

2126

type ( state_type ), intent (INOUT) :: state

2127

2128

!*****************************************************************************

2129

! internal variables

2130

2131

integer :: i

2132

2133

real(dp), dimension (3) :: force, r, s, t

2134

2135

!*****************************************************************************

2136

! begin subroutine

2137

2138

if ( debug ) write(*,*) 'Entering transformForces2Cartesian()'

2139

2140

if ( state % representation % force ) then ! forces already Cartesian; do nothing

2141

2142

return

2143

2144

else ! forces in lattice rep., so transform

2145

2146

r = (/ state % cell % a(1), state % cell % b(1), state % cell % c(1) /)

2147

s = (/ state % cell % a(2), state % cell % b(2), state % cell % c(2) /)

2148

t = (/ state % cell % a(3), state % cell % b(3), state % cell % c(3) /)

2149

2150

do i=1, state % nAtoms

2151

2152

force = state % force(1:3,i)

2153

2154

state % force(1,i) = dot_product( r, force(1:3) )

2155

state % force(2,i) = dot_product( s, force(1:3) )

2156

state % force(3,i) = dot_product( t, force(1:3) )

2157

2158

end do

2159

2160

state % representation % force = .true.

2161

2162

end if

2163

2164

if ( debug ) write(*,*) 'Exiting transformForces2Cartesian()'

2165

2166

end subroutine transformForces2Cartesian

2167

2168

!*****************************************************************************

2169

subroutine transformMomenta2Lattice( state )

2170

!*****************************************************************************

2171

!

2172

! Project: MolecularDynamics

2173

!

2174

! Module: State

2175

!

2176

! Created on: Thu 8 Mar 12:23:32 2018 by Eduardo R. Hernandez

2177

!

2178

! Last Modified: Wed 18 Jul 14:49:27 2018

2179

!

2180

! ---

2181

2182

! This file is distributed under the terms of the

2183

! GNU General Public License: see COPYING in the top directory

2184

! or http://www.gnu.org/copyleft/gpl.txt .

2185

! ---

2186

!

2187

! This subroutine transforms the atomic momenta from Cartesian to Lattice

2188

! representation. Internally, State stores all atomic variables in Cartesian

2189

! representation, but they may be requested by the calling program in

2190

! Lattice representation; this routine performs the appropriate transformation

2191

! for momenta.

2192

!

2193

!*****************************************************************************

2194

! modules used

2195

2196

!*****************************************************************************

2197

! No implicits please!

2198

2199

implicit none

2200

2201

!*****************************************************************************

2202

! interface variables

2203

2204

type ( state_type ), intent (INOUT) :: state

2205

2206

!*****************************************************************************

2207

! local variables

2208

2209

integer :: i

2210

2211

real (dp), dimension (3) :: vector

2212

2213

!*****************************************************************************

2214

! begin subroutine

2215

2216

if ( debug ) write(*,*) 'Entering transformMomenta2Lattice()'

2217

2218

if ( state % representation % momentum ) then ! momenta are in cartesian, so transform

2219

2220

do i=1, state % nAtoms

2221

2222

vector(1) = dot_product( state % reciprocalCell % a, &

2223

state % momentum(1:3,i) )

2224

vector(2) = dot_product( state % reciprocalCell % b, &

2225

state % momentum(1:3,i) )

2226

vector(3) = dot_product( state % reciprocalCell % c, &

2227

state % momentum(1:3,i) )

2228

2229

state % momentum(1:3,i) = matmul( state % cell % metricTensor, vector )

2230

2231

end do

2232

2233

state % representation % momentum = .false. ! now they are in lattice

2234

2235

else ! momenta already in lattice rep, so nothing to do

2236

2237

return

2238

2239

end if

2240

2241

if ( debug ) write(*,*) 'Exiting transformMomenta2Lattice()'

2242

2243

end subroutine transformMomenta2Lattice

2244

2245

!*****************************************************************************

2246

subroutine transformMomenta2Cartesian( state )

2247

!*****************************************************************************

2248

!

2249

! Project: MolecularDynamics

2250

!

2251

! Created on: Thu 8 Mar 12:23:32 2018 by Eduardo R. Hernandez

2252

!

2253

! Last Modified: Wed 18 Jul 14:34:42 2018

2254

!

2255

! ---

2256

2257

! This file is distributed under the terms of the

2258

! GNU General Public License: see COPYING in the top directory

2259

! or http://www.gnu.org/copyleft/gpl.txt .

2260

! ---

2261

!

2262

!*****************************************************************************

2263

! modules used

2264

2265

!*****************************************************************************

2266

! No implicits please!

2267

2268

implicit none

2269

2270

!*****************************************************************************

2271

! interface variables

2272

2273

type ( state_type ), intent (INOUT) :: state

2274

2275

!*****************************************************************************

2276

! local variables

2277

2278

integer :: i

2279

2280

real(dp), dimension (3) :: p, r, s, t

2281

2282

!*****************************************************************************

2283

! begin subroutine

2284

2285

if ( debug ) write(*,*) 'Entering transformMomenta2Cartesian()'

2286

2287

if ( state % representation % momentum ) then ! momenta already cartesian; do nothing

2288

2289

return

2290

2291

else ! no, we need to transform to cartesian, so do it here

2292

2293

r = (/ state % cell % a(1), state % cell % b(1), state % cell % c(1) /)

2294

s = (/ state % cell % a(2), state % cell % b(2), state % cell % c(2) /)

2295

t = (/ state % cell % a(3), state % cell % b(3), state % cell % c(3) /)

2296

2297

do i=1, state % nAtoms

2298

2299

p = matmul( state % reciprocalCell % metricTensor, state % momentum(1:3,i) )

2300

2301

state % momentum(1,i) = dot_product( r, p )

2302

state % momentum(2,i) = dot_product( s, p )

2303

state % momentum(3,i) = dot_product( t, p )

2304

2305

end do

2306

2307

state % representation % momentum = .true.

2308

2309

end if

2310

2311

if ( debug ) write(*,*) 'Exiting transformMomenta2Cartesian()'

2312

2313

end subroutine transformMomenta2Cartesian

2314

2315

!*****************************************************************************

2316

subroutine updateCell( Cell, reciprocalCell , LatticeVectors, MetricTensor )

2317

!*****************************************************************************

2318

!

2319

! Project: MolecularDynamics

2320

!

2321

! Module: State

2322

!

2323

! Created on: Mon 19 Mar 10:53:13 2018 by Eduardo R. Hernandez

2324

!

2325

! Last Modified: Tue 3 Jul 17:30:42 2018

2326

!

2327

! ---

2328

2329

! This file is distributed under the terms of the

2330

! GNU General Public License: see COPYING in the top directory

2331

! or http://www.gnu.org/copyleft/gpl.txt .

2332

! ---

2333

!

2334

! From the current lattice vectors, or alternatively, the metric tensor,

2335

! this subroutine updates (recalculates) the

2336

! current cell (cell vectors, angles, moduli, etc) as well as its reciprocal

2337

! cell; this is only needed in case of flexible-cell dynamics.

2338

!

2339

!*****************************************************************************

2340

! modules used

2341

2342

!*****************************************************************************

2343

! No implicits please!

2344

2345

implicit none

2346

2347

!*****************************************************************************

2348

! arguments

2349

2350

type ( lattice_type ), intent (INOUT) :: Cell

2351

type ( lattice_type ), intent (INOUT) :: reciprocalCell

2352

2353

real (dp), dimension (3,3), intent (IN), optional :: LatticeVectors

2354

! a lattice vector per column, a(1:3) = LatticeVectors(1:3,1), etc

2355

real (dp), dimension (3,3), intent (IN), optional :: MetricTensor

2356

2357

!*****************************************************************************

2358

! local variables

2359

2360

integer :: i

2361

2362

double precision :: cos_angle, det_metric, factor

2363

2364

double precision :: vector(3)

2365

2366

double precision :: matrix_a(3,3)

2367

double precision :: matrix_b(3,3)

2368

2369

!*****************************************************************************

2370

! start of subroutine

2371

2372

if ( debug ) write(*,*) 'Entering updateCell()'

2373

2374

if ( present( LatticeVectors ) ) then

2375

2376

Cell % a(1:3) = LatticeVectors(1:3,1)

2377

Cell % b(1:3) = LatticeVectors(1:3,2)

2378

Cell % c(1:3) = LatticeVectors(1:3,3)

2379

2380

Cell % metricTensor(1,1) = dot_product( Cell % a, Cell % a )

2381

Cell % metricTensor(1,2) = dot_product( Cell % a, Cell % b )

2382

Cell % metricTensor(1,3) = dot_product( Cell % a, Cell % c )

2383

2384

Cell % metricTensor(2,1) = dot_product( Cell % b, Cell % a )

2385

Cell % metricTensor(2,2) = dot_product( Cell % b, Cell % b )

2386

Cell % metricTensor(2,3) = dot_product( Cell % b, Cell % c )

2387

2388

Cell % metricTensor(3,1) = dot_product( Cell % c, Cell % a )

2389

Cell % metricTensor(3,2) = dot_product( Cell % c, Cell % b )

2390

Cell % metricTensor(3,3) = dot_product( Cell % c, Cell % c )

2391

2392

else if ( present( MetricTensor ) ) then

2393

2394

Cell % metricTensor = MetricTensor

2395

2396

else if ( ( .not. present( LatticeVectors ) ) .and. ( .not. present( MetricTensor ) ) ) then

2397

2398

write(*,*) 'Error in updateCell: either lattice vectors or metric tensor must be provided!'

2399

stop

2400

2401

end if

2402

2403

Cell % a(1) = sqrt( Cell % metricTensor(1,1) )

2404

Cell % a(2) = zero

2405

Cell % a(3) = zero

2406

2407

Cell % b(1) = Cell % metricTensor(1,2) / sqrt( Cell % metricTensor(1,1) )

2408

Cell % b(2) = sqrt( Cell % metricTensor(2,2) - &

2409

Cell % metricTensor(1,2) * Cell % metricTensor(1,2) / &

2410

Cell % metricTensor(1,1) )

2411

Cell % b(3) = zero

2412

2413

Cell % c(1) = Cell % metricTensor(1,3) / sqrt( Cell % metricTensor(1,1) )

2414

Cell % c(2) = ( Cell % metricTensor(1,1) * Cell % metricTensor(2,3) - &

2415

Cell % metricTensor(1,2) * Cell % metricTensor(1,3) ) / &

2416

sqrt( Cell % metricTensor(1,1) * Cell % metricTensor(1,1) * &

2417

Cell % metricTensor(2,2) - Cell % metricTensor(1,1) * &

2418

Cell % metricTensor(1,2) * Cell % metricTensor(1,2) )

2419

Cell % c(3) = sqrt( Cell % metricTensor(3,3) - &

2420

Cell % c(1) * Cell % c(1) - Cell % c(2) * Cell % c(2) )

2421

2422

Cell % a_modulus = sqrt( dot_product( Cell % a, Cell % a ) )

2423

Cell % b_modulus = sqrt( dot_product( Cell % b, Cell % b ) )

2424

Cell % c_modulus = sqrt( dot_product( Cell % c, Cell % c ) )

2425

2426

cos_angle = dot_product( Cell % b, Cell % c ) / &

2427

( Cell % b_modulus * Cell % c_modulus )

2428

2429

Cell % alpha = acos ( cos_angle ) * rad_to_deg

2430

2431

cos_angle = dot_product( Cell % a, Cell % c ) / &

2432

( Cell % a_modulus * Cell % c_modulus )

2433

2434

Cell % beta = acos ( cos_angle ) * rad_to_deg

2435

2436

cos_angle = dot_product( Cell % a, Cell % b ) / &

2437

( Cell % a_modulus * Cell % b_modulus )

2438

2439

Cell % gamma = acos ( cos_angle ) * rad_to_deg

2440

2441

vector(1) = Cell % b(2) * Cell % c(3) - Cell % c(2) * Cell % b(3)

2442

vector(2) = Cell % b(3) * Cell % c(1) - Cell % c(3) * Cell % b(1)

2443

vector(3) = Cell % b(1) * Cell % c(2) - Cell % c(1) * Cell % b(2)

2444

2445

Cell % volume = dot_product( vector, Cell % a )

2446

2447

! calculate the reciprocal lattice vectors

2448

! these are the reciprocal lattice vectors but for a factor of 2pi

2449

2450

reciprocalCell % a = vector / Cell % volume

2451

2452

vector(1) = Cell % c(2) * Cell % a(3) - Cell % a(2) * Cell % c(3)

2453

vector(2) = Cell % c(3) * Cell % a(1) - Cell % a(3) * Cell % c(1)

2454

vector(3) = Cell % c(1) * Cell % a(2) - Cell % a(1) * Cell % c(2)

2455

2456

reciprocalCell % b = vector / Cell % volume

2457

2458

vector(1) = Cell % a(2) * Cell % b(3) - Cell % b(2) * Cell % a(3)

2459

vector(2) = Cell % a(3) * Cell % b(1) - Cell % b(3) * Cell % a(1)

2460

vector(3) = Cell % a(1) * Cell % b(2) - Cell % b(1) * Cell % a(2)

2461

2462

reciprocalCell % c = vector / Cell % volume

2463

2464

reciprocalCell % a_modulus = &

2465

sqrt( dot_product( reciprocalCell % a, reciprocalCell % a ) )

2466

reciprocalCell % b_modulus = &

2467

sqrt( dot_product( reciprocalCell % b, reciprocalCell % b ) )

2468

reciprocalCell % c_modulus = &

2469

sqrt( dot_product( reciprocalCell % c, reciprocalCell % c ) )

2470

2471

cos_angle = dot_product( reciprocalCell % b, reciprocalCell % c ) / &

2472

( reciprocalCell % b_modulus * reciprocalCell % c_modulus )

2473

2474

reciprocalCell % alpha = acos ( cos_angle ) * rad_to_deg

2475

2476

cos_angle = dot_product( reciprocalCell % a, reciprocalCell % c ) / &

2477

( reciprocalCell % a_modulus * reciprocalCell % c_modulus )

2478

2479

reciprocalCell % beta = acos ( cos_angle ) * rad_to_deg

2480

2481

cos_angle = dot_product( reciprocalCell % a, reciprocalCell % b ) / &

2482

( reciprocalCell % a_modulus * reciprocalCell % b_modulus )

2483

2484

reciprocalCell % gamma = acos ( cos_angle ) * rad_to_deg

2485

2486

reciprocalCell % volume = one / Cell % volume

2487

2488

! calculate the reciprocal metric tensor

2489

2490

reciprocalCell % metricTensor(1,1) = &

2491

dot_product( reciprocalCell % a, reciprocalCell % a )

2492

reciprocalCell % metricTensor(1,2) = &

2493

dot_product( reciprocalCell % a, reciprocalCell % b )

2494

reciprocalCell % metricTensor(1,3) = &

2495

dot_product( reciprocalCell % a, reciprocalCell % c )

2496

2497

reciprocalCell % metricTensor(2,1) = &

2498

dot_product( reciprocalCell % b, reciprocalCell % a )

2499

reciprocalCell % metricTensor(2,2) = &

2500

dot_product( reciprocalCell % b, reciprocalCell % b )

2501

reciprocalCell % metricTensor(2,3) = &

2502

dot_product( reciprocalCell % b, reciprocalCell % c )

2503

2504

reciprocalCell % metricTensor(3,1) = &

2505

dot_product( reciprocalCell % c, reciprocalCell % a )

2506

reciprocalCell % metricTensor(3,2) = &

2507

dot_product( reciprocalCell % c, reciprocalCell % b )

2508

reciprocalCell % metricTensor(3,3) = &

2509

dot_product( reciprocalCell % c, reciprocalCell % c )

2510

2511

if ( debug ) write(*,*) 'Exiting updateCell()'

2512

2513

end subroutine updateCell

2514

2515

!*****************************************************************************

2516

subroutine readRestart( state, restartFile, nSteps )

2517

!*****************************************************************************

2518

!

2519

! Project: MolecularDynamics

2520

!

2521

! Created on: Fri 2 Nov 12:33:45 2018 by Eduardo R. Hernandez

2522

!

2523

! Last Modified: Mon 5 Nov 18:17:59 2018

2524

!

2525

! ---

2526

2527

! This file is distributed under the terms of the

2528

! GNU General Public License: see COPYING in the top directory

2529

! or http://www.gnu.org/copyleft/gpl.txt .

2530

! ---

2531

!

2532

!*****************************************************************************

2533

! modules used

2534

2535

!*****************************************************************************

2536

! No implicits please!

2537

2538

implicit none

2539

2540

!*****************************************************************************

2541

! arguments

2542

2543

type ( state_type ), intent (INOUT) :: state ! the state to be read

2544

2545

character ( len = 30 ), intent (IN) :: restartFile

2546

2547

integer, intent (INOUT), optional :: nSteps

2548

2549

!*****************************************************************************

2550

! local variables

2551

2552

integer :: i

2553

integer :: n

2554

2555

!*****************************************************************************

2556

! begin subroutine

2557

2558

if ( debug ) write(*,*) 'Entering readRestart()'

2559

2560

open( unit = 9, file = restartFile )

2561

2562

! if the number of steps has been passed, write it; else write a zero

2563

! this is so that the continuation run can start counting steps from a

2564

! the counter of a previous run, if desired

2565

2566

read(9,*) nSteps

2567

2568

! write the H0 constant, needed for calculating the constant of motion

2569

2570

read(9,*) state % H0 % E

2571

state % H0 % set = .true.

2572

2573

! write the representation flags

2574

2575

read(9,*) state % representation % position, &

2576

state % representation % momentum, &

2577

state % representation % force

2578

2579

! now write the number of atoms, positions, forces, momenta and mass

2580

2581

read(9,*) state % nAtoms

2582

read(9,*) state % nDegreesOfFreedom

2583

2584

do i=1, state % nAtoms

2585

read(9,*) state % position(1:3,i)

2586

end do

2587

2588

do i=1, state % nAtoms

2589

read(9,*) state % momentum(1:3,i)

2590

end do

2591

2592

do i=1, state % nAtoms

2593

read(9,*) state % mass(i)

2594

end do

2595

2596

! write the cell metric tensor; with this, the rest of the cell info

2597

! can be easily reconstructed

2598

2599

read(9,*) state % cell % metricTensor(1:3,1)

2600

read(9,*) state % cell % metricTensor(1:3,2)

2601

read(9,*) state % cell % metricTensor(1:3,3)

2602

2603

call updateCell( state % cell, state % reciprocalCell, &

2604

MetricTensor = state % cell % metricTensor )

2605

2606

! write associated barostat variables

2607

2608

read(9,*) state % barostat % mass

2609

2610

read(9,*) state % barostat % momentum(1:3,1)

2611

read(9,*) state % barostat % momentum(1:3,2)

2612

read(9,*) state % barostat % momentum(1:3,3)

2613

2614

! now write the thermostat information

2615

2616

read(9,*) state % nThermostats

2617

2618

do i=1, state % nThermostats

2619

2620

read(9,*) state % thermostat(i) % mass, &

2621

state % thermostat(i) % position, &

2622

state % thermostat(i) % momentum, &

2623

state % thermostat(i) % a, &

2624

state % thermostat(i) % C

2625

2626

end do

2627

2628

close( unit = 9 )

2629

2630

if ( debug ) write(*,*) 'Exiting readRestart()'

2631

2632

end subroutine readRestart

2633

!*****************************************************************************

2634

subroutine writeRestart( state, restartFileName, nSteps )

2635

!*****************************************************************************

2636

!

2637

! Project: MolecularDynamics

2638

!

2639

! Created on: Fri 2 Nov 11:05:47 2018 by Eduardo R. Hernandez

2640

!

2641

! Last Modified: Tue 6 Nov 09:40:31 2018

2642

!

2643

! ---

2644

2645

! This file is distributed under the terms of the

2646

! GNU General Public License: see COPYING in the top directory

2647

! or http://www.gnu.org/copyleft/gpl.txt .

2648

! ---

2649

!

2650

! This subroutine writes to a file the current state of the system, with the

2651

! purpose of serving as a continuation file for subsequent runs.

2652

!

2653

!*****************************************************************************

2654

! modules used

2655

2656

!*****************************************************************************

2657

! No implicits please!

2658

2659

implicit none

2660

2661

!*****************************************************************************

2662

! arguments

2663

2664

type ( state_type ), intent (IN) :: state ! the state to be stored

2665

2666

character ( len = 30 ), intent (IN) :: restartFileName

2667

2668

integer, intent (IN), optional :: nSteps

2669

2670

!*****************************************************************************

2671

! local variables

2672

2673

integer :: i

2674

integer :: n

2675

2676

!*****************************************************************************

2677

! begin subroutine

2678

2679

if ( debug ) write(*,*) 'Entering writeRestart()'

2680

2681

open( unit = 9, file = restartFileName )

2682

2683

! if the number of steps has been passed, write it; else write a zero

2684

! this is so that the continuation run can start counting steps from a

2685

! the counter of a previous run, if desired

2686

2687

if ( present( nSteps ) ) then

2688

n = nSteps

2689

else

2690

n = 0

2691

end if

2692

2693

write(9,*) n

2694

2695

! write the H0 constant, needed for calculating the constant of motion

2696

2697

write(9,*) state % H0 % E

2698

2699

! write the representation flags

2700

2701

write(9,*) state % representation % position, &

2702

state % representation % momentum, &

2703

state % representation % force

2704

2705

! now write the number of atoms, positions, forces, momenta and mass

2706

2707

write(9,*) state % nAtoms

2708

write(9,*) state % nDegreesOfFreedom

2709

2710

do i=1, state % nAtoms

2711

write(9,*) state % position(1:3,i)

2712

end do

2713

2714

do i=1, state % nAtoms

2715

write(9,*) state % momentum(1:3,i)

2716

end do

2717

2718

do i=1, state % nAtoms

2719

write(9,*) state % mass(i)

2720

end do

2721

2722

! write the cell metric tensor; with this, the rest of the cell info

2723

! can be easily reconstructed

2724

2725

write(9,*) state % cell % metricTensor(1:3,1)

2726

write(9,*) state % cell % metricTensor(1:3,2)

2727

write(9,*) state % cell % metricTensor(1:3,3)

2728

2729

! write associated barostat variables

2730

2731

write(9,*) state % barostat % mass

2732

2733

write(9,*) state % barostat % momentum(1:3,1)

2734

write(9,*) state % barostat % momentum(1:3,2)

2735

write(9,*) state % barostat % momentum(1:3,3)

2736

2737

! now write the thermostat information

2738

2739

write(9,*) state % nThermostats

2740

2741

do i=1, state % nThermostats

2742

2743

write(9,*) state % thermostat(i) % mass, &

2744

state % thermostat(i) % position, &

2745

state % thermostat(i) % momentum, &

2746

state % thermostat(i) % a, &

2747

state % thermostat(i) % C

2748

2749

end do

2750

2751

close( unit = 9 )

2752

2753

if ( debug ) write(*,*) 'Exiting writeRestart()'

2754

2755

end subroutine writeRestart

2756

2757

end module md_state_module