-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathneural_network.F90
217 lines (185 loc) · 6.51 KB
/
neural_network.F90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
subroutine force_nnp_matlantis
!!! === Nproc much be 1 at current verion === !!!
use Parameters, &
only: Eenergy, r, fr, Natom, AUtoAng, eVtoAU, dp_inv, alabel, &
addresstmp, Ista, Iend, laddress, MyRank, Nbead, eVAng2AU, &
Lperiodic
use utility, only: program_abort
implicit none
integer :: Nfile, Imode, i, j
integer :: Uinp, Uout, ios
character(len=13), allocatable :: Fout(:)
character(len=12) :: char_num
character(len=:), allocatable :: command, Fforce, Fenergy
integer :: access, Iaccess
Fforce = 'forces.out'
Fenergy = 'energy.out'
Nfile = Iend - Ista + 1
allocate(Fout(Nfile))
if ( Lperiodic .eqv. .True. ) then
call input_periodic
else
call input_nonperio
end if
if (access(Fforce, ' ') == 0) call system('rm -f '//Fforce)
if (access(Fenergy, ' ') == 0) call system('rm -f '//Fenergy)
! +++ Execution Matlantis +++
write(char_num,'(" ",I0, " ",I0, L2)') Ista, Iend, Lperiodic
command = ' ./run_matlantis.py '//trim(addresstmp)//trim(char_num)
call system(command)
! +++ Execution Matlantis +++
!if ( MyRank == 0 ) then
open(newunit=Uinp,file=Fforce,status='old',iostat=ios)
if ( ios /= 0 ) then
call program_abort('ERROR!!! There is no "'//Fforce//'"')
end if
do j = 1, Nbead
do i = 1, Natom
read(Uinp,*) fr(:,i,j)
end do
end do
close(Uinp)
open(newunit=Uinp,file=Fenergy,status='old',iostat=ios)
if ( ios /= 0 ) then
call program_abort('There is no "'//Fenergy//'"')
end if
do j = 1, Nbead
read(Uinp,*) Eenergy(j)
end do
close(Uinp)
fr(:,:,:) = fr(:,:,:) * eVAng2AU * dp_inv
Eenergy(:) = Eenergy(:) * eVtoAU
!end if
contains
subroutine input_periodic
do Imode = Ista, Iend
write(Fout(Imode),'("str",I5.5,".vasp")') Imode
end do
do Imode = Ista, Iend
call system( 'cp LATTICE '//trim(addresstmp)//trim(Fout(Imode)) )
open(newunit=Uout,file=trim(addresstmp)//trim(Fout(Imode)),status='old',position='append')
do i = 1, Natom
write(Uout,*) r(:,i,Imode) * AUtoAng
end do
write(Uout,*)
close(Uout)
end do
end subroutine input_periodic
subroutine input_nonperio
do Imode = Ista, Iend
write(Fout(Imode),'("str",I5.5,".xyz")') Imode
end do
do Imode = Ista, Iend
open( newunit=Uout,file=trim(addresstmp)//trim(Fout(Imode)) )
write(Uout,*) Natom
write(Uout,*) 'Properties=species:S:1:pos:R:3 pbc="F F F"'
do i = 1, Natom
write(Uout,*) alabel(i), r(:,i,Imode)*AUtoAng
end do
close(Uout)
end do
end subroutine input_nonperio
end subroutine force_nnp_matlantis
subroutine set_nnp_matlantis
use Parameters, only: Lperiodic, Nproc
use utility, only: program_abort
implicit none
character(len=:), allocatable :: name_file
integer :: access
if ( Nproc /= 1 ) then
call program_abort('ERROR!!! current version only for Nproc == 1')
end if
name_file = 'run_matlantis.py'
if ( access(name_file,' ') /= 0 ) then
call program_abort('ERROR!!! There is no "'//name_file//'"')
end if
if ( Lperiodic .eqv. .True. ) then
if ( access("./LATTICE", " ") /= 0) call err_LATTICE
end if
contains
subroutine err_LATTICE
integer :: Nunit
open(newunit=Nunit,file='LATTICE',status='new')
write(Nunit,'(a)') trim("H2O molecule ")
write(Nunit,'(a)') trim("1.0 ")
write(Nunit,'(a)') trim(" 7.9376997948 0.0000000000 0.0000000000")
write(Nunit,'(a)') trim(" 0.0000000000 7.9376997948 0.0000000000")
write(Nunit,'(a)') trim(" 0.0000000000 0.0000000000 7.9376997948")
write(Nunit,'(a)') trim(" O H ")
write(Nunit,'(a)') trim(" 1 2 ")
write(Nunit,'(a)') trim("Cartesian ")
close(Nunit)
print *, 'ERROR!! "LATTICE" is NOT exist!'
print *, 'Please use "LATTICE" template'
call program_abort('')
end subroutine err_LATTICE
end subroutine set_nnp_matlantis
subroutine set_nnp_araidai
use Parameters
implicit none
Integer :: i,j,k
integer :: Imode
do Imode = Ista, Iend
write(addresstmp(laddress+1:laddress+6),'(i5.5,A1)') Imode,'/'
call system('mkdir -p '//trim(addresstmp))
call system('cp -r nnp_files '//trim(addresstmp))
call system('cp training.data_1 '//trim(addresstmp))
call system('cp n2training '//trim(addresstmp))
enddo
return
end subroutine set_nnp_araidai
subroutine force_nnp_araidai
use Parameters, &
only: Eenergy, r, fr, Natom, AUtoAng, eVtoAU, dp_inv, alabel, &
addresstmp, Ista, Iend, laddress, eVAng2AU
use utility, only: program_abort
implicit none
character(Len=130) :: line, inp_train(8)
integer :: Imode, i, Uout, Uinp, ios
character :: dummyC
!real(8), parameter :: eVAng_HartBohr = 0.5291772108d0 / 27.21138505d0
open(newunit=Uinp,file='training.data_1',status='old',iostat=ios)
if ( ios /= 0) then
call program_abort('ERROR!!: There is no "training.data_1"')
end if
do i = 1, 5
read(Uinp,'(a)') inp_train(i)
end do
do i = 1, Natom
read(Uinp,'()')
end do
do i = 6, 8
read(Uinp,'(a)') inp_train(i)
end do
close(Uinp)
do Imode = Ista, Iend
write(addresstmp(laddress+1:laddress+6),'(i5.5,A1)') Imode,'/'
open(newunit=Uout,file=trim(addresstmp)//'training.data_1',status='replace')
do i = 1, 5
write(Uout,'(a)') trim(inp_train(i))
end do
do i=1,Natom
write(Uout,9998) &
"atom", r(:,i,Imode)*AUtoAng, alabel(i), 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0
enddo
do i = 6, 8
write(Uout,'(a)') trim(inp_train(i))
end do
close(Uout)
call system('cd '//trim(addresstmp)//' ; ./n2training > /dev/null ; cd ../.. ')
!call system('cd '//trim(addresstmp)//' ; mpiexec.hydra -n 1 ./n2training ; cd ../.. ')
open(newunit=Uinp,file=trim(addresstmp)//'foce_npp.out',status='old')
read(Uinp,*) Eenergy(Imode)
do i = 1, Natom
read(Uinp,*) dummyC, fr(:,i,Imode)
end do
close(Uinp)
fr(:,:,Imode)=fr(:,:,Imode)*eVAng2AU*dp_inv
Eenergy(Imode) = Eenergy(Imode) * eVtoAU
end do
return
9999 format(3F24.16)
9998 format(A, 3E17.9, x, A, 6E12.4)
end subroutine force_nnp_araidai
subroutine force_nnp_aenet
end subroutine force_nnp_aenet