1 subroutine w2x(nvarr,x,*)
4 include "DIMENSIONS.ZSCOPT"
5 include "COMMON.CONTROL"
6 include "COMMON.INTERACT"
7 include "COMMON.TORSION"
8 include "COMMON.WEIGHTS"
9 include "COMMON.ENERGIES"
10 include "COMMON.IOUNITS"
11 include "COMMON.NAMES"
13 include "COMMON.VMCPAR"
14 include "COMMON.CLASSES"
15 include "COMMON.WEIGHTDER"
16 include "COMMON.FFIELD"
17 include "COMMON.SCCOR"
18 integer itor2typ(-2:2) /-20,-20,10,9,20/
20 double precision x(max_paropt),xchuj,epsi
21 integer nvarr,i,ii,j,k,l,ll,iprot,iblock
26 if (imask(i).eq.1) then
35 c Change 9/21/07 AL Variables are epsilons and not the increments on epsilons
36 c to enable control on their deviations from original values.
37 c Single-body epsilons are now defined as epsilons of pairs of same type side
40 write (iout,*) "W2X: mod_side",mod_side," nsingle_sc",nsingle_sc,
41 & " npair_sc",npair_sc
46 eps_orig(i,j)=eps(i,j)
54 c write (iout,*) "ii",ii," x",x(ii)
57 c if (epss_up(i).eq.0.0d0) then
58 c x_up(ii)=dabs(epsi)*epss_low(i)
59 c x_low(ii)=epsi-x_up(ii)
60 c x_up(ii)=epsi+x_up(ii)
62 c x_low(ii)=epss_low(i)
72 if (ityp_ssc(j).eq.iti) then
73 write (iout,*) "Error - pair ",restyp(iti),"-",
75 & " specified in variable epsilon list but type defined",
76 & " in single-type list."
85 c if (epsp_up(i).eq.0.0d0) then
86 c x_up(ii)=dabs(epsi)*epsp_low(i)
87 c x_low(ii)=epsi-x_up(ii)
88 c x_up(ii)=epsi+x_up(ii)
90 c x_low(ii)=epsp_low(i)
93 c write (iout,*) "ii",ii," x",x(ii)
97 if (mod_fourier(nloctyp).gt.0) then
102 if (mod_fourier(i).gt.0) then
106 if (mask_bnew1(k,j,i).gt.0) then
109 x_low(ii)=bnew1_low(k,j,i)
110 x_up(ii)=bnew1_up(k,j,i)
117 if (mask_bnew2(k,j,i).gt.0) then
120 x_low(ii)=bnew2_low(k,j,i)
121 x_up(ii)=bnew2_up(k,j,i)
128 if (mask_ccnew(k,j,i).gt.0) then
131 x_low(ii)=ccnew_low(k,j,i)
132 x_up(ii)=ccnew_up(k,j,i)
139 if (mask_ddnew(k,j,i).gt.0) then
142 x_low(ii)=ddnew_low(k,j,i)
143 x_up(ii)=ddnew_up(k,j,i)
149 if (mask_e0new(j,i).gt.0) then
152 x_low(ii)=e0new_low(j,i)
153 x_up(ii)=e0new_up(j,i)
160 if (mask_eenew(l,k,j,i).gt.0) then
163 x_low(ii)=eenew_low(l,k,j,i)
164 x_up(ii)=eenew_up(l,k,j,i)
176 if (mod_fourier(i).gt.0) then
178 if (mask_fourier(j,i).gt.0) then
196 write (iout,*) "Indices of optimizable torsionals"
197 do i=-ntortyp+1,ntortyp-1
198 do j=-ntortyp+1,ntortyp-1
199 c write (iout,*) "mask",i,j,mask_tor(j,i)
201 if (tor_mode.eq.0) then
204 if (mask_tor(0,j,i,iblock).eq.1) then
207 x(ii)=weitor(0,j,i,iblock)
208 x_low(ii)=weitor_low(0,j,i,iblock)
209 x_up(ii)=weitor_up(0,j,i,iblock)
210 write (iout,'(2a4,2i4)') restyp(itor2typ(j)),
211 & restyp(itor2typ(i)),iblock,ii
212 else if (mask_tor(0,j,i,iblock).eq.2) then
213 do k=1,nterm(j,i,iblock)
216 x(ii)=v1(k,j,i,iblock)
217 x_low(ii)=weitor_low(0,j,i,iblock)
218 x_up(ii)=weitor_up(0,j,i,iblock)
219 write (iout,'(2a4,3i4)') restyp(itor2typ(j)),
220 & restyp(itor2typ(i)),iblock,k,ii
222 else if (mask_tor(0,j,i,iblock).gt.2) then
223 do k=1,nterm(j,i,iblock)
226 x(ii)=v1(k,j,i,iblock)
227 x_low(ii)=weitor_low(0,j,i,iblock)
228 x_up(ii)=weitor_up(0,j,i,iblock)
229 write (iout,'(2a4,3i4)') restyp(itor2typ(j)),
230 & restyp(itor2typ(i)),iblock,k,ii
232 do k=1,nterm(j,i,iblock)
235 x(ii)=v2(k,j,i,iblock)
236 x_low(ii)=weitor_low(0,j,i,iblock)
237 x_up(ii)=weitor_up(0,j,i,iblock)
238 write (iout,'(2a4,3i4)') restyp(itor2typ(j)),
239 & restyp(itor2typ(i)),ii,k,ii
244 else if (tor_mode.eq.1) then
246 if (mask_tor(0,j,i,1).eq.1) then
249 x(ii)=weitor(0,j,i,1)
250 x_low(ii)=weitor_low(0,j,i,1)
251 x_up(ii)=weitor_up(0,j,i,1)
252 write (iout,'(2a4,i4)') restyp(itor2typ(j)),
253 & restyp(itor2typ(i)),ii
254 else if (mask_tor(0,j,i,1).eq.2) then
255 do k=1,nterm_kcc(j,i)
256 do l=1,nterm_kcc_Tb(j,i)
257 do ll=1,nterm_kcc_Tb(j,i)
260 x(ii)=v1_kcc(ll,l,k,j,i)
261 x_low(ii)=weitor_low(0,j,i,1)
262 x_up(ii)=weitor_up(0,j,i,1)
263 write (iout,'(2a4,2i4)') restyp(itor2typ(j)),
264 & restyp(itor2typ(i)),k,ii
268 else if (mask_tor(0,j,i,1).gt.2) then
269 do k=1,nterm_kcc(j,i)
270 do l=1,nterm_kcc_Tb(j,i)
271 do ll=1,nterm_kcc_Tb(j,i)
274 x(ii)=v1_kcc(ll,l,k,j,i)
275 x_low(ii)=weitor_low(0,j,i,1)
276 x_up(ii)=weitor_up(0,j,i,1)
277 write (iout,'(2a4,2i4)') restyp(itor2typ(j)),
278 & restyp(itor2typ(i)),k,ii
282 do k=1,nterm_kcc(j,i)
283 do l=1,nterm_kcc_Tb(j,i)
284 do ll=1,nterm_kcc_Tb(j,i)
287 x(ii)=v2_kcc(ll,l,k,j,i)
288 x_low(ii)=weitor_low(0,j,i,1)
289 x_up(ii)=weitor_up(0,j,i,1)
290 write (iout,'(2a4,2i4)') restyp(itor2typ(j)),
291 & restyp(itor2typ(i)),k,ii
298 write (iout,*) "ERROR! TOR_MODE>1 ALLOWED ONLY WITH NEWCORR"
306 nbacktor_var=ntor_var
307 write (iout,*) "nbacktor_var",nbacktor_var
311 write (iout,*) "Indices of optimizable SCCOR torsionals"
312 do i=-nsccortyp,nsccortyp
313 do j=-nsccortyp,nsccortyp
315 c write (iout,*) "mask",i,j,mask_tor(j,i)
316 if (mask_tor(k,j,i,1).eq.1) then
319 x(ii)=weitor(k,j,i,1)
320 x_low(ii)=weitor_low(k,j,i,1)
321 x_up(ii)=weitor_up(k,j,i,1)
322 write (iout,'(i4,2a4,i4)') k,restyp(j),restyp(i),ii
323 else if (mask_tor(k,j,i,1).eq.2) then
324 do l=1,nterm_sccor(j,i)
327 x(ii)=v1sccor(l,k,j,i)
328 x_low(ii)=weitor_low(k,j,i,1)
329 x_up(ii)=weitor_up(k,j,i,1)
330 write (iout,'(i4,2a4,2i4)') k,restyp(j),restyp(i),l,ii
332 else if (mask_tor(k,j,i,1).gt.2) then
333 do l=1,nterm_sccor(j,i)
336 x(ii)=v1sccor(l,k,j,i)
337 x_low(ii)=weitor_low(k,j,i,1)
338 x_up(ii)=weitor_up(k,j,i,1)
339 write (iout,'(i4,2a4,2i4)') k,restyp(j),restyp(i),l,ii
341 do l=1,nterm_sccor(j,i)
344 x(ii)=v2sccor(l,k,j,i)
345 x_low(ii)=weitor_low(k,j,i,1)
346 x_up(ii)=weitor_up(k,j,i,1)
347 write (iout,'(i4,2a4,2i4)') k,restyp(j),restyp(i),l,ii
356 c Added SC sigmas and anisotropies
358 if (mod_side_other) then
360 if (mask_sigma(1,i).gt.0) then
363 x_low(ii)=sigma_low(1,i)
364 x_up(ii)=sigma_up(1,i)
365 c write (iout,*) "sig0",i,ii,x(ii), x_low(ii),x_up(ii)
367 if (mask_sigma(2,i).gt.0) then
370 x_low(ii)=sigma_low(2,i)
371 x_up(ii)=sigma_up(2,i)
372 c write (iout,*) "sigi",i,ii,x(ii), x_low(ii),x_up(ii)
374 if (mask_sigma(3,i).gt.0) then
377 x_low(ii)=sigma_low(3,i)
378 x_up(ii)=sigma_up(3,i)
379 c write (iout,*) "chi",i,ii,x(ii), x_low(ii),x_up(ii)
381 if (mask_sigma(4,i).gt.0) then
384 x_low(ii)=sigma_low(4,i)
385 x_up(ii)=sigma_up(4,i)
386 c write (iout,*) "alp",i,ii,x(ii), x_low(ii),x_up(ii)
396 if (mask_elec(i,j,1).gt.0) then
399 x_low(ii)=epp_low(i,j)
402 if (mask_elec(i,j,2).gt.0) then
405 x_low(ii)=rpp_low(i,j)
408 if (mask_elec(i,j,3).gt.0) then
411 x_low(ii)=elpp6_low(i,j)
412 x_up(ii)=elpp6_up(i,j)
414 if (mask_elec(i,j,4).gt.0) then
417 x_low(ii)=elpp3_low(i,j)
418 x_up(ii)=elpp3_up(i,j)
425 C Define the SC-p interaction constants
429 if (mask_scp(0,1,1)+mask_scp(0,2,1)+mask_scp(0,1,2)+
430 & mask_scp(0,2,2) .gt. 0) then
432 if (mask_scp(0,j,1).gt.0) then
435 x_low(ii)=epscp_low(0,j)
436 x_up(ii)=epscp_up(0,j)
438 if (mask_scp(0,j,2).gt.0) then
441 x_low(ii)=rscp_low(0,j)
442 x_up(ii)=rscp_up(0,j)
448 if (mask_scp(i,j,1).gt.0) then
451 x_low(ii)=epscp_low(i,j)
452 x_up(ii)=epscp_up(i,j)
454 if (mask_scp(i,j,2).gt.0) then
457 x_low(ii)=rscp_low(i,j)
458 x_up(ii)=rscp_up(i,j)
467 c Compute boundaries, if defined relative to starting values
469 write (iout,*) "i",i," x_low",x_low(i)," x_up",x_up(i)
470 if (x_up(i).eq.0.0d0 .and. x_low(i).gt.0.0d0) then
474 else if (x_up(i).eq.-1.0d0 .and. x_low(i).gt.0.0d0) then
476 if (x(i).gt.0.0d0) then
477 x_low(i)=x(i)*(1.0d0-xchuj)
478 x_up(i)=x(i)*(1.0d0+xchuj)
480 x_low(i)=x(i)*(1.0d0+xchuj)
481 x_up(i)=x(i)*(1.0d0-xchuj)
486 c Base of heat capacity
489 if (nbeta(2,iprot).gt.0) then
491 x(ii)=heatbase(iprot)
499 write (iout,*) "Number of optimizable parameters:",nvarr
501 write (iout,*) "Initial variables and boundaries:"
502 write (iout,'(3x,3a10)') "x","x_low","x_up"
504 write (iout,'(i3,3f10.5)') i,x(i),x_low(i),x_up(i)
509 c------------------------------------------------------------------------------
510 subroutine x2w(nvarr,x)
513 include "DIMENSIONS.ZSCOPT"
514 include "COMMON.CONTROL"
515 include "COMMON.NAMES"
516 include "COMMON.WEIGHTS"
517 include "COMMON.INTERACT"
518 include "COMMON.ENERGIES"
519 include "COMMON.FFIELD"
520 include "COMMON.TORSION"
521 include "COMMON.IOUNITS"
522 include "COMMON.VMCPAR"
523 include "COMMON.CLASSES"
524 include "COMMON.WEIGHTDER"
525 include "COMMON.SCCOR"
527 integer i,j,k,kk,l,ll,ii,nvarr,iti,itj,iprot,iblock
528 double precision rri,v0ij,si
529 double precision x(nvarr)
530 double precision epsij,rrij,sigeps,sigt1sq,sigt2sq,sigii1,sigii2,
531 & ratsig1,ratsig2,rsum_max,dwa16,r_augm
532 double precision eps_temp(ntyp,ntyp)
533 double precision ftune_eps,ftune_epsprim
534 logical in_sc_pair_list
535 external in_sc_pair_list
536 integer itor2typ(3) /10,9,20/
540 if (lprint) write (iout,*) "x2w: nvar",nvarr
541 if (lprint) write (iout,*) "x",(x(i),i=1,nvarr)
544 if (imask(i).eq.1) then
576 write (iout,*) "Weights:"
578 write (iout,'(a,f10.5)') wname(i),ww(i)
589 if (.not. in_sc_pair_list(iti,j)) then
590 eps(iti,j)=eps_orig(iti,j)
591 & +0.5d0*(x(ii)-eps_orig(iti,iti))
592 eps(j,iti)=eps(iti,j)
605 if (mod_fourier(nloctyp).gt.0) then
610 if (mod_fourier(i).gt.0) then
614 if (mask_bnew1(k,j,i).gt.0) then
623 if (mask_bnew2(k,j,i).gt.0) then
632 if (mask_ccnew(k,j,i).gt.0) then
641 if (mask_ddnew(k,j,i).gt.0) then
649 if (mask_e0new(j,i).gt.0) then
658 if (mask_eenew(l,k,j,i).gt.0) then
672 bnew1(j,1,-i)= bnew1(j,1,i)
673 bnew1(j,2,-i)=-bnew1(j,2,i)
674 bnew2(j,1,-i)= bnew2(j,1,i)
675 bnew2(j,2,-i)=-bnew2(j,2,i)
678 c ccnew(j,1,i)=ccnew(j,1,i)/2
679 c ccnew(j,2,i)=ccnew(j,2,i)/2
680 c ddnew(j,1,i)=ddnew(j,1,i)/2
681 c ddnew(j,2,i)=ddnew(j,2,i)/2
682 ccnew(j,1,-i)=ccnew(j,1,i)
683 ccnew(j,2,-i)=-ccnew(j,2,i)
684 ddnew(j,1,-i)=ddnew(j,1,i)
685 ddnew(j,2,-i)=-ddnew(j,2,i)
687 e0new(1,-i)= e0new(1,i)
688 e0new(2,-i)=-e0new(2,i)
689 e0new(3,-i)=-e0new(3,i)
691 eenew(kk,1,1,-i)= eenew(kk,1,1,i)
692 eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
693 eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
694 eenew(kk,2,2,-i)= eenew(kk,2,2,i)
698 write (iout,'(a)') "Coefficients of the multibody terms"
699 do i=-nloctyp+1,nloctyp-1
700 write (iout,*) "Type: ",onelet(iloctyp(i))
701 write (iout,*) "Coefficients of the expansion of B1"
703 write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
705 write (iout,*) "Coefficients of the expansion of B2"
707 write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
709 write (iout,*) "Coefficients of the expansion of C"
710 write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
711 write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
712 write (iout,*) "Coefficients of the expansion of D"
713 write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
714 write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
715 write (iout,*) "Coefficients of the expansion of E"
716 write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
719 write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
728 if (mod_fourier(i).gt.0) then
731 if (mask_fourier(j,i).gt.0) then
738 write (iout,*) "Fourier coefficients of cumulants"
739 write (iout,*) 'Type',i
740 write (iout,'(a,i2,a,f10.5)') ('b(',j,')=',b(j,i),j=1,13)
751 c B1tilde(1,i) = b(3,i)
752 c B1tilde(2,i) =-b(5,i)
759 CCold(1,1,-i)= b(7,i)
760 CCold(2,2,-i)=-b(7,i)
761 CCold(2,1,-i)=-b(9,i)
762 CCold(1,2,-i)=-b(9,i)
763 c Ctilde(1,1,i)=b(7,i)
764 c Ctilde(1,2,i)=b(9,i)
765 c Ctilde(2,1,i)=-b(9,i)
766 c Ctilde(2,2,i)=b(7,i)
771 DDold(1,1,-i)= b(6,i)
772 DDold(2,2,-i)=-b(6,i)
773 DDold(2,1,-i)=-b(8,i)
774 DDold(1,2,-i)=-b(8,i)
775 c Dtilde(1,1,i)=b(6,i)
776 c Dtilde(1,2,i)=b(8,i)
777 c Dtilde(2,1,i)=-b(8,i)
778 c Dtilde(2,2,i)=b(6,i)
779 EEold(1,1,i)= b(10,i)+b(11,i)
780 EEold(2,2,i)=-b(10,i)+b(11,i)
781 EEold(2,1,i)= b(12,i)-b(13,i)
782 EEold(1,2,i)= b(12,i)+b(13,i)
783 EEold(1,1,i)= b(10,i)+b(11,i)
784 EEold(2,2,i)=-b(10,i)+b(11,i)
785 EEold(2,1,i)= b(12,i)-b(13,i)
786 EEold(1,2,i)= b(12,i)+b(13,i)
787 EEold(1,1,-i)= b(10,i)+b(11,i)
788 EEold(2,2,-i)=-b(10,i)+b(11,i)
789 EEold(2,1,-i)=-b(12,i)+b(13,i)
790 EEold(1,2,-i)=-b(12,i)-b(13,i)
798 &"Coefficients of the cumulants (independent of valence angles)"
799 do i=-nloctyp+1,nloctyp-1
800 write (iout,*) 'Type ',onelet(iloctyp(i))
802 write(iout,'(2f10.5)') B(3,i),B(5,i)
804 write(iout,'(2f10.5)') B(2,i),B(4,i)
807 write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
811 write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
815 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
823 if (mod_tor .or. mod_fourier(nloctyp).gt.0) then
824 c write (iout,*) "w2x: mod_tor",mod_tor," mod_fourier",
825 c & mod_fourier(nloctyp)
829 do i=-ntortyp+1,ntortyp-1
830 do j=-ntortyp+1,ntortyp-1
831 if (tor_mode.eq.0) then
833 if (mask_tor(0,j,i,iblock).eq.1) then
835 weitor(0,j,i,iblock)=x(ii)
836 else if (mask_tor(0,j,i,iblock).eq.2) then
837 do k=1,nterm(j,i,iblock)
839 v1(k,j,i,iblock)=x(ii)
841 else if (mask_tor(0,j,i,iblock).gt.2) then
842 do k=1,nterm(j,i,iblock)
844 v1(k,j,i,iblock)=x(ii)
846 do k=1,nterm(j,i,iblock)
848 v2(k,j,i,iblock)=x(ii)
853 do k=1,nterm(i,j,iblock)
854 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
855 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
856 v0ij=v0ij+si*v1(k,i,j,iblock)
859 do k=1,nlor(i,j,iblock)
860 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
863 v0(-i,-j,iblock)=v0ij
865 else if (tor_mode.eq.1) then
866 if (mask_tor(0,j,i,1).eq.1) then
868 weitor(0,j,i,1)=x(ii)
869 else if (mask_tor(0,j,i,1).eq.2) then
870 do k=1,nterm_kcc(j,i)
871 do l=1,nterm_kcc_Tb(j,i)
872 do ll=1,nterm_kcc_Tb(j,i)
874 v1_kcc(ll,l,k,j,i)=x(ii)
878 else if (mask_tor(0,j,i,1).gt.2) then
879 do k=1,nterm_kcc(j,i)
880 do l=1,nterm_kcc_Tb(j,i)
881 do ll=1,nterm_kcc_Tb(j,i)
883 v1_kcc(ll,l,k,j,i)=x(ii)
887 do k=1,nterm_kcc(j,i)
888 do l=1,nterm_kcc_Tb(j,i)
889 do ll=1,nterm_kcc_Tb(j,i)
891 v2_kcc(ll,l,k,j,i)=x(ii)
898 c Compute torsional coefficients from local energy surface expansion coeffs
899 if (mask_tor(0,j,i,1).eq.1) then
901 weitor(0,j,i,1)=x(ii)
902 c else if (mask_tor(0,j,i,1).eq.2) then
904 do k=1,nterm_kcc_Tb(j,i)
905 do l=1,nterm_kcc_Tb(j,i)
906 v1_kcc(k,l,1,i,j)=bnew1(k,1,i)*bnew2(l,1,j)
907 & +bnew1(k,2,i)*bnew2(l,2,j)
908 v2_kcc(k,l,1,i,j)=bnew1(k,1,i)*bnew2(l,2,j)
909 & +bnew1(k,2,i)*bnew2(l,1,j)
912 do k=1,nterm_kcc_Tb(j,i)
913 do l=1,nterm_kcc_Tb(j,i)
914 v1_kcc(k,l,2,i,j)=-0.25*(ccnew(k,1,i)*ddnew(l,1,j)
915 & -ccnew(k,2,i)*ddnew(l,2,j))
916 v2_kcc(k,l,2,i,j)=-0.25*(ccnew(k,2,i)*ddnew(l,1,j)
917 & +ccnew(k,1,i)*ddnew(l,2,j))
923 write (iout,*) "ERROR: TOR_MODE>1 ALLOWED ONLY WITH NEWCORR"
933 do i=-nsccortyp,nsccortyp
934 do j=-nsccortyp,nsccortyp
936 if (mask_tor(k,j,i,1).eq.1) then
938 weitor(k,j,i,1)=x(ii)
939 else if (mask_tor(k,j,i,1).eq.2) then
942 v1sccor(l,k,j,i)=x(ii)
944 else if (mask_tor(k,j,i,1).gt.2) then
945 do l=1,nterm_sccor(j,i)
947 v1sccor(l,k,j,i)=x(ii)
949 do l=1,nterm_sccor(j,i)
951 v2sccor(l,k,j,i)=x(ii)
959 c SC sigmas and anisotropies as parameters
960 if (mod_side_other) then
962 c write (iout,*) "mask",i,(mask_sigma(k,i),k=1,4)
965 if (mask_sigma(1,i).gt.0) then
968 c write(iout,*) "sig0",i,ii,x(ii)
970 if (mask_sigma(2,i).gt.0) then
973 c write(iout,*) "sigi",i,ii,x(ii)
975 if (mask_sigma(3,i).gt.0) then
978 c write(iout,*) "chip",i,ii,x(ii)
980 if (mask_sigma(4,i).gt.0) then
983 c write(iout,*) "alp",i,ii,x(ii)
988 if (mod_side.or.mod_side_other) then
989 if (lprint) write(iout,'(/3a,2i3)') 'Potential is ',
990 & potname(ipot),', exponents are ',expon,2*expon
991 C-----------------------------------------------------------------------
992 C Calculate the "working" parameters of SC interactions.
993 dwa16=2.0d0**(1.0d0/6.0d0)
996 chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
1001 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1002 sigma(j,i)=sigma(i,j)
1003 rs0(i,j)=dwa16*sigma(i,j)
1007 if (lprint) call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1008 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1009 & 'Working parameters of the SC interactions:',
1010 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1015 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1024 c sigeps=dsign(1.0D0,epsij)
1025 aa(i,j)=ftune_eps(epsij)*rrij*rrij
1030 sigt1sq=sigma0(i)**2
1031 sigt2sq=sigma0(j)**2
1034 ratsig1=sigt2sq/sigt1sq
1035 ratsig2=1.0D0/ratsig1
1036 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1037 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1038 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1042 sigmaii(i,j)=rsum_max
1043 sigmaii(j,i)=rsum_max
1044 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max)
1046 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1047 augm(i,j)=epsij*r_augm**(2*expon)
1048 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1055 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1056 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1057 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1068 if (mask_elec(i,j,1).gt.0) then
1071 if (epp(i,j).lt.0.0d0) epp(i,j)=0.0d0
1074 if (mask_elec(i,j,2).gt.0) then
1076 rpp(i,j)=dabs(x(ii))
1077 rpp(j,i)=dabs(x(ii))
1079 if (mask_elec(i,j,3).gt.0) then
1082 if (elpp6(i,j).gt.0.0d0) elpp6(i,j)=0.0d0
1083 elpp6(j,i)=elpp6(i,j)
1085 if (mask_elec(i,j,4).gt.0) then
1091 app (i,j)=epp(i,j)*rri*rri
1092 app (j,i)=epp(j,i)*rri*rri
1093 bpp (i,j)=-2.0D0*epp(i,j)*rri
1094 bpp (j,i)=-2.0D0*epp(i,j)*rri
1095 ael6(i,j)=elpp6(i,j)*4.2D0**6
1096 ael6(j,i)=elpp6(i,j)*4.2D0**6
1097 ael3(i,j)=elpp3(i,j)*4.2D0**3
1098 ael3(j,i)=elpp3(i,j)*4.2D0**3
1102 write (iout,'(/a)') 'Electrostatic interaction constants:'
1103 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1104 & 'IT','JT','APP','BPP','AEL6','AEL3'
1107 write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1108 & ael6(i,j),ael3(i,j)
1115 C Define the SC-p interaction constants
1119 if (mask_scp(0,1,1)+mask_scp(0,2,1)+mask_scp(0,1,2)+
1120 & mask_scp(0,2,2) .gt. 0) then
1122 if (mask_scp(0,j,1).gt.0) then
1128 if (mask_scp(0,j,2).gt.0) then
1138 if (mask_scp(i,j,1).gt.0) then
1142 if (mask_scp(i,j,2).gt.0) then
1150 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1151 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1152 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1153 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1157 write (iout,*) "Parameters of SC-p interactions:"
1159 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1160 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1167 c write (iout,*) "x2w: nvar",nvarr
1169 c Base of heat capacity
1172 if (nbeta(2,iprot).gt.0) then
1174 heatbase(iprot)=x(ii)
1178 if (ii.ne.nvarr) then
1179 write (iout,*) "!!!!! SEVERE ERROR IN X2W !!!!! (chuj nastapil)"
1180 write (iout,*) "NVAR=",nvarr," LAST VARIABLE INDEX=",ii
1181 stop "CHUJ ABSOLUTNY!!!!"
1186 c---------------------------------------------------------------------------
1187 double precision function ftune_eps(x)
1189 double precision x,y
1190 double precision delta /1.0d-5/
1191 if (dabs(x).lt.delta) then
1193 ftune_eps=delta*(3.0d0+6.0d0*y**2-y**4)/8
1199 c---------------------------------------------------------------------------
1200 double precision function ftune_epsprim(x)
1202 double precision x,y
1203 double precision delta /1.0d-5/
1204 if (dabs(x).lt.delta) then
1206 ftune_epsprim=1.5d0*y-0.5*y**3
1208 ftune_epsprim=dsign(1.0d0,x)
1212 c---------------------------------------------------------------------------
1213 double precision function ftune_epsbis(x)
1215 double precision x,y
1216 double precision delta /1.0d-5/
1217 if (dabs(x).lt.delta) then
1219 ftune_epsbis=(1.5d0-1.5d0*y**3)/delta
1225 c---------------------------------------------------------------------------
1226 logical function in_sc_pair_list(iti,itj)
1228 include "DIMENSIONS"
1229 include "DIMENSIONS.ZSCOPT"
1230 include "COMMON.WEIGHTS"
1233 if (ityp_psc(1,i).eq.iti .and. ityp_psc(2,i).eq.itj .or.
1234 & ityp_psc(1,i).eq.itj .and. ityp_psc(2,i).eq.iti) then
1235 in_sc_pair_list=.true.
1239 in_sc_pair_list=.false.