1
1
! > High-level program to perform SFC cycles in Extended Huckel Hamiltonian.
2
- ! ! This program takes coordinates in xyz or pdb format and extracts information
3
- ! ! about the system.
2
+ ! ! This program takes coordinates in xyz or pdb format and extracts information
3
+ ! ! about the system.
4
4
program gpscf
5
5
6
6
! BML lib.
7
- use bml
7
+ use bml
8
8
! Progress and LATTE lib modules.
9
9
use prg_progress_mod
10
10
use prg_system_mod
11
11
use prg_ptable_mod
12
12
use latteparser_latte_mod
13
13
use huckel_latte_mod
14
- use tbparams_latte_mod
14
+ use tbparams_latte_mod
15
15
use ham_latte_mod
16
16
use coulomb_latte_mod
17
17
use prg_charges_mod
@@ -30,7 +30,7 @@ program gpscf
30
30
use prg_subgraphLoop_mod
31
31
use prg_homolumo_mod
32
32
33
- implicit none
33
+ implicit none
34
34
35
35
integer , parameter :: dp = kind (1.0d0 )
36
36
integer :: i, j, nel, norb
@@ -56,56 +56,56 @@ program gpscf
56
56
call prg_progress_init()
57
57
58
58
if (printRank() .eq. 1 ) then
59
- write (* ,* ) " gptest start ..."
59
+ write (* ,* ) " gptest start ..."
60
60
endif
61
61
62
- ! > Parsing input file. This file contains all the variables needed to
63
- ! run the scf including the sp2 (solver) variables. lt is "latte_type" structure
64
- ! containing all the variables.
62
+ ! > Parsing input file. This file contains all the variables needed to
63
+ ! run the scf including the sp2 (solver) variables. lt is "latte_type" structure
64
+ ! containing all the variables.
65
65
! file://~/progress/build/doc/html/structlatteparser__latte__mod_1_1latte__type.html
66
- call parse_latte(lt," input.in" )
66
+ call parse_latte(lt," input.in" )
67
67
68
- ! > Parsing system coordinates. This reads the coords.pdb file to get the position of every
68
+ ! > Parsing system coordinates. This reads the coords.pdb file to get the position of every
69
69
! atom in the system. sy is the "system_type" structure containing all the variables.
70
70
! file://~/progress/build/doc/html/classsystem__latte__mod.html
71
- call prg_parse_system(sy," coords" ," pdb" )
71
+ call prg_parse_system(sy," coords" ," pdb" )
72
72
73
73
! > Allocate bounds vactor.
74
74
allocate (gbnd(2 ))
75
75
76
- ! > Get Huckel hamiltonian. Computes the Extended Huckel Hamiltonian from the
76
+ ! > Get Huckel hamiltonian. Computes the Extended Huckel Hamiltonian from the
77
77
! atom coordinates. The main inputs are the huckelTBparams and the system coordinate (sy%coordinate)
78
78
! The main outputs are Hamiltonian (ham_bml) and Overlap (over_bml) matrices.
79
79
call get_hshuckel(ham_bml,over_bml,sy% coordinate,sy% spindex,sy% spatnum,&
80
- " ../../huckelTBparams" ,lt% bml_type,lt% mdim,lt% threshold&
81
- ,tb% nsp,tb% splist,tb% basis,tb% numel,tb% onsite_energ,&
82
- tb% norbi,tb% hubbardu)
80
+ " ../../huckelTBparams" ,lt% bml_type,lt% mdim,lt% threshold&
81
+ ,tb% nsp,tb% splist,tb% basis,tb% numel,tb% onsite_energ,&
82
+ tb% norbi,tb% hubbardu)
83
83
84
- ! > Get the mapping of the Hamiltonian index with the atom index
84
+ ! > Get the mapping of the Hamiltonian index with the atom index
85
85
! hindex(1,i)=starting Hindex for atom i.
86
86
! hindex(2,i)=final Hindex for atom i.
87
87
! file://~/progress/build/doc/html/ham__latte__mod_8F90_source.html
88
- call get_hindex(sy% spindex,tb% norbi,hindex,norb)
89
-
88
+ call get_hindex(sy% spindex,tb% norbi,hindex,norb)
89
+
90
90
if (printRank() .eq. 1 ) then
91
91
write (* ,* ) " Number of orbitals = " , norb
92
92
write (* ,* )
93
93
call bml_print_matrix(" ham0_bml" ,ham_bml,0 ,6 ,0 ,6 )
94
94
endif
95
95
96
- ! > Get occupation based on last shell population.
96
+ ! > Get occupation based on last shell population.
97
97
! WARNING: This could change depending on the TB method being used.
98
98
nel = sum (element_numel(sy% atomic_number(:)),&
99
- size (sy% atomic_number,dim= 1 ))
99
+ size (sy% atomic_number,dim= 1 ))
100
100
bndfil = nel/ (2.0_dp * norb)
101
101
if (printRank() .eq. 1 ) then
102
102
write (* ,* ) " bndfil = " , bndfil
103
103
write (* ,* ) " nel = " , nel
104
104
endif
105
105
106
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
107
- ! > First Charge computation
108
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
106
+ ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
107
+ ! > First Charge computation
108
+ ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
109
109
110
110
! > Initialize the density matrix (rho_bml) and inverse overlap factor (zmat_bml).
111
111
call prg_init_pzmat(rho_bml,zmat_bml,lt% bml_type,lt% mdim,norb)
@@ -119,15 +119,15 @@ program gpscf
119
119
! > Orthogonalize ham.
120
120
call prg_timer_start(ortho_timer)
121
121
call prg_orthogonalize(ham_bml,zmat_bml,orthoh_bml,&
122
- lt% threshold,lt% bml_type,lt% verbose)
122
+ lt% threshold,lt% bml_type,lt% verbose)
123
123
call prg_timer_stop(ortho_timer)
124
124
125
125
#ifdef DO_MPI_BLOCK
126
126
call prg_allGatherParallel(orthoh_bml)
127
127
#endif
128
128
129
- ! > The SP2 algorithm is used to get the first orthogonal Density matrix (orthop).
130
- call prg_parse_gsp2(gsp2," input.in" )
129
+ ! > The SP2 algorithm is used to get the first orthogonal Density matrix (orthop).
130
+ call prg_parse_gsp2(gsp2," input.in" )
131
131
132
132
! > Calculate gershgorin bounds
133
133
call bml_gershgorin(orthoh_bml, gbnd)
@@ -139,7 +139,7 @@ program gpscf
139
139
! > SP2 algorithm.
140
140
call prg_timer_start(sp2_timer)
141
141
call prg_sp2_alg2_genseq(orthoh_bml,orthop_bml,lt% threshold,bndfil,gsp2% minsp2iter,&
142
- gsp2% maxsp2iter,gsp2% sp2conv,gsp2% sp2tol, pp, icount, vv)
142
+ gsp2% maxsp2iter,gsp2% sp2conv,gsp2% sp2tol, pp, icount, vv)
143
143
call prg_timer_stop(sp2_timer)
144
144
#ifdef DO_MPI_BLOCK
145
145
call prg_allGatherParallel(orthop_bml)
@@ -166,20 +166,20 @@ program gpscf
166
166
call bml_zero_matrix(lt% bml_type,bml_element_real,dp,norb,norb,g_bml)
167
167
call bml_copy(orthop_bml, g_bml)
168
168
169
- ! > Deprg_orthogonalize rho.
169
+ ! > Deprg_orthogonalize rho.
170
170
call prg_timer_start(deortho_timer)
171
171
call prg_deorthogonalize(orthop_bml,zmat_bml,rho_bml,&
172
- lt% threshold,lt% bml_type,lt% verbose)
172
+ lt% threshold,lt% bml_type,lt% verbose)
173
173
call prg_timer_stop(deortho_timer)
174
174
#ifdef DO_MPI_BLOCK
175
175
call prg_allGatherParallel(rho_bml)
176
176
#endif
177
177
178
178
if (printRank() .eq. 1 ) then
179
- call bml_print_matrix(" rho_bml" ,rho_bml,0 ,6 ,0 ,6 )
179
+ call bml_print_matrix(" rho_bml" ,rho_bml,0 ,6 ,0 ,6 )
180
180
endif
181
181
182
- ! > Get charges based on rho. rho_bml is the input and sy%net_charge is the outputs vector containing
182
+ ! > Get charges based on rho. rho_bml is the input and sy%net_charge is the outputs vector containing
183
183
! the charges.
184
184
call prg_get_charges(rho_bml, over_bml, hindex, sy% net_charge, tb% numel, sy% spindex, lt% mdim, lt% threshold)
185
185
charges_old = sy% net_charge
@@ -188,9 +188,9 @@ program gpscf
188
188
write (* ,* )" Total charges =" , sum (sy% net_charge(:),size (sy% net_charge,dim= 1 ))
189
189
endif
190
190
191
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
191
+ ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
192
192
! > SCF loop
193
- ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
193
+ ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
194
194
195
195
! > Save actual hamiltonian as the non-scf Hamiltonian (H0)
196
196
call bml_copy_new(ham_bml,ham0_bml)
@@ -199,7 +199,7 @@ program gpscf
199
199
call prg_get_recip_vects(sy% lattice_vector,sy% recip_vector,sy% volr,sy% volk)
200
200
201
201
! > Beginning of the SCF loop.
202
- do i= 1 ,lt% maxscf
202
+ do i= 1 ,lt% maxscf
203
203
204
204
if (printRank() .eq. 1 ) then
205
205
write (* ,* )" SCF iter" , i
@@ -210,34 +210,34 @@ program gpscf
210
210
write (* ,* )" In real Coul ..."
211
211
endif
212
212
call get_ewald_real(sy% spindex,sy% splist,sy% coordinate&
213
- ,sy% net_charge,tb% hubbardu,sy% lattice_vector,&
214
- sy% volr,lt% coul_acc,coul_forces_r,coul_pot_r);
213
+ ,sy% net_charge,tb% hubbardu,sy% lattice_vector,&
214
+ sy% volr,lt% coul_acc,coul_forces_r,coul_pot_r);
215
215
216
216
! > Reciprocal contribution to the Coul energy. The outputs are coul_forces_k,coul_pot_k.
217
217
if (printRank() .eq. 1 ) then
218
- write (* ,* )" In recip Coul ..."
218
+ write (* ,* )" In recip Coul ..."
219
219
endif
220
220
call get_ewald_recip(sy% spindex,sy% splist,sy% coordinate&
221
- ,sy% net_charge,tb% hubbardu,sy% lattice_vector,&
222
- sy% recip_vector,sy% volr,lt% coul_acc,coul_forces_k,coul_pot_k);
221
+ ,sy% net_charge,tb% hubbardu,sy% lattice_vector,&
222
+ sy% recip_vector,sy% volr,lt% coul_acc,coul_forces_k,coul_pot_k);
223
223
224
224
! > Get the scf hamiltonian. The outputs is ham_bml.
225
225
if (printRank() .eq. 1 ) then
226
226
write (* ,* )" in prg_get_hscf ..."
227
227
endif
228
228
call prg_get_hscf(ham0_bml,over_bml,ham_bml,sy% spindex,hindex,tb% hubbardu,sy% net_charge,&
229
- coul_pot_r,coul_pot_k,lt% mdim,lt% threshold)
229
+ coul_pot_r,coul_pot_k,lt% mdim,lt% threshold)
230
230
231
231
! > Initialize the orthogonal versions of ham and rho.
232
232
call prg_init_ortho(orthoh_bml,orthop_bml,lt% bml_type,lt% mdim,norb)
233
233
234
234
! > Orthogonalize the Hamiltonian
235
- ! if (printRank() .eq. 1) then
236
- write (* ,* )" in prg_orthogonalize H ..."
237
- ! endif
235
+ ! if (printRank() .eq. 1) then
236
+ write (* ,* )" in prg_orthogonalize H ..."
237
+ ! endif
238
238
call prg_timer_start(ortho_timer)
239
239
call prg_orthogonalize(ham_bml,zmat_bml,orthoh_bml,&
240
- lt% threshold,lt% bml_type,lt% verbose)
240
+ lt% threshold,lt% bml_type,lt% verbose)
241
241
call prg_timer_stop(ortho_timer)
242
242
#ifdef DO_MPI_BLOCK
243
243
call prg_allGatherParallel(orthoh_bml)
@@ -246,7 +246,7 @@ program gpscf
246
246
if (printRank() .eq. 1 ) then
247
247
call bml_print_matrix(" orthoh_bml" ,orthoh_bml,0 ,6 ,0 ,6 )
248
248
call bml_print_matrix(" ham_bml" ,ham_bml,0 ,6 ,0 ,6 )
249
- call bml_print_matrix(" zmat_bml" ,zmat_bml,0 ,6 ,0 ,6 )
249
+ call bml_print_matrix(" zmat_bml" ,zmat_bml,0 ,6 ,0 ,6 )
250
250
endif
251
251
252
252
! > Threshold the graph
@@ -260,27 +260,27 @@ program gpscf
260
260
call prg_equalPartition(gp, gsp2% nodesPerPart, norb)
261
261
call prg_timer_stop(part_timer)
262
262
263
- ! call prg_timer_start(dyn_timer,"metis+SA")
264
- ! call metis
265
- ! call SA
266
- ! call prg_timer_start(dyn_timer,1)
263
+ ! call prg_timer_start(dyn_timer,"metis+SA")
264
+ ! call metis
265
+ ! call SA
266
+ ! call prg_timer_start(dyn_timer,1)
267
+
268
+ ! call prg_timer_start(dyn_timer,"metisSA")
269
+ ! call system("echo 'hello'")
270
+ ! call prg_timer_stop(dyn_timer,1)
267
271
268
- ! call prg_timer_start(dyn_timer,"metisSA")
269
- ! call system("echo 'hello'")
270
- ! call prg_timer_stop(dyn_timer,1)
271
-
272
272
! > Calculate gershgorin bounds
273
273
call bml_gershgorin(orthoh_bml, gbnd)
274
274
gp% mineval = gbnd(1 )
275
- gp% maxeval = gbnd(2 )
275
+ gp% maxeval = gbnd(2 )
276
276
if (printRank() .eq. 1 ) then
277
277
write (* ,* ) " Gershgorin: mineval = " , gbnd(1 ), " maxeval = " , gbnd(2 )
278
278
write (* ,* )
279
279
endif
280
280
281
281
! ! Calculate SP2 sequence
282
282
call prg_sp2sequence(gp% pp, gp% maxIter, gp% mineval, gp% maxeval, ehomo, elumo, &
283
- gsp2% errlimit)
283
+ gsp2% errlimit)
284
284
if (printRank() .eq. 1 ) then
285
285
write (* ,* )
286
286
write (* ,* ) " SP2Sequence: Max iterations = " , gp% maxIter
@@ -290,9 +290,9 @@ program gpscf
290
290
! > Now use the graph-based SP2 algorithm to get the orthogonal Density
291
291
! matrix.
292
292
call prg_timer_start(graphsp2_timer)
293
-
294
- call prg_subgraphSP2Loop(orthoh_bml, g_bml, orthop_bml, gp, lt% threshold)
295
- ! call prg_sp2_alg1_seq(orthoh_bml,orthop_bml,lt%threshold, gp%pp, gp%maxIter, gp%vv)
293
+
294
+ call prg_subgraphSP2Loop(orthoh_bml, g_bml, orthop_bml, gp, lt% threshold)
295
+ ! call prg_sp2_alg1_seq(orthoh_bml,orthop_bml,lt%threshold, gp%pp, gp%maxIter, gp%vv)
296
296
297
297
call prg_timer_stop(graphsp2_timer)
298
298
#ifdef DO_MPI_BLOCK
@@ -328,7 +328,7 @@ program gpscf
328
328
! > Deprg_orthogonalize orthop_bml to get the density matrix rho_bml.
329
329
call prg_timer_start(deortho_timer)
330
330
call prg_deorthogonalize(orthop_bml,zmat_bml,rho_bml,&
331
- lt% threshold,lt% bml_type,lt% verbose)
331
+ lt% threshold,lt% bml_type,lt% verbose)
332
332
call prg_timer_stop(deortho_timer)
333
333
#ifdef DO_MPI_BLOCK
334
334
call prg_allGatherParallel(rho_bml)
@@ -343,38 +343,38 @@ program gpscf
343
343
344
344
! > Get the system charges.
345
345
call prg_get_charges(rho_bml,over_bml,hindex,sy% net_charge,tb% numel,&
346
- sy% spindex,lt% mdim,lt% threshold)
346
+ sy% spindex,lt% mdim,lt% threshold)
347
347
348
348
! if (printRank() .eq. 1) then
349
- write (* ,* )" Total charge" , sum (sy% net_charge(:)),size (sy% net_charge,dim= 1 )
349
+ write (* ,* )" Total charge" , sum (sy% net_charge(:)),size (sy% net_charge,dim= 1 )
350
350
! endif
351
351
352
352
call prg_qmixer(sy% net_charge,charges_old,dqin,&
353
- dqout,scferror,i,lt% pulaycoeff,lt% mpulay,lt% verbose)
353
+ dqout,scferror,i,lt% pulaycoeff,lt% mpulay,lt% verbose)
354
354
355
355
charges_old = sy% net_charge
356
-
356
+
357
357
if (printRank() .eq. 1 ) then
358
358
write (* ,* )" System charges:"
359
359
do j= 1 ,4
360
360
write (* ,* )sy% symbol(j),sy% net_charge(j)
361
361
enddo
362
362
endif
363
- if (scferror.lt. lt% scftol.and. i.gt. 5 ) then
363
+ if (scferror.lt. lt% scftol.and. i.gt. 5 ) then
364
364
if (printRank() .eq. 1 ) then
365
365
write (* ,* )" SCF converged within" ,i," steps ..."
366
366
write (* ,* )" SCF error =" ,scferror
367
367
endif
368
368
exit
369
- endif
369
+ endif
370
370
371
371
enddo
372
372
! > End of SCF loop.
373
373
374
374
allocate (eigenvals(norb))
375
375
call prg_get_eigenvalues(ham_bml,eigenvals,lt% verbose)
376
376
call prg_write_tdos(eigenvals, 0.01_dp , 1000 , - 25.0_dp , 20.0_dp , " dos.dos" )
377
-
377
+
378
378
deallocate (gbnd)
379
379
380
380
call bml_deallocate(ham_bml)
0 commit comments