& iscpstart(maxres,maxint_gr),iscpend(maxres,maxint_gr),
& iatsc_s,iatsc_e,iatel_s,iatel_e,iatel_s_vdw,iatel_e_vdw,
& iatscp_s,iatscp_e,ispp,iscp
+C 3/26/20 Interaction lists
+ integer newcontlisti(200*maxres),newcontlistj(200*maxres),
+ & newcontlistppi(200*maxres),newcontlistppj(200*maxres),
+ & newcontlistscpi(200*maxres),newcontlistscpj(200*maxres),
+ & g_listscsc_start,g_listscsc_end,g_listpp_start,g_listpp_end,
+ & g_listscp_start,g_listscp_end
+ common /interact_list/newcontlisti,newcontlistj,g_listscsc_start,
+ & g_listscsc_end,newcontlistppi,newcontlistppj,g_listpp_start,
+ & g_listpp_end,newcontlistscpi,newcontlistscpj,g_listscp_start,
+ & g_listscp_end
C 12/1/95 Array EPS included in the COMMON block.
double precision eps,epslip,sigma,sigmaii,rs0,chi,chip,alp,
& sigma0,sigii,
& dvmax,damax,edriftmax
integer n_timestep,ntwx,ntwe,lang,count_reset_moment,
& count_reset_vel,ntime_split,ntime_split0,
- & maxtime_split
+ & maxtime_split,itime_mat,imatupdate
logical large,print_compon,tbf,rest,reset_moment,reset_vel,
& rattle,mdpdb,RESPA,preminim
common /mdpar/ v_ini,d_time,d_time0,t_bath,
& tau_bath,dvmax,damax,n_timestep,mdpdb,
- & ntime_split,ntime_split0,maxtime_split,
+ & ntime_split,ntime_split0,maxtime_split,itime_mat,imatupdate,
& ntwx,ntwe,lang,large,print_compon,tbf,rest,preminim,
& reset_moment,reset_vel,count_reset_moment,count_reset_vel,
& rattle,RESPA
& .and. mod(itime,count_reset_vel).eq.0) then
call random_vel
write(iout,'(a,f20.2)')
- & "Velocities reset to random values, time",totT
+ & "Velocities reset to random values, time",totT
do i=0,2*nres
do j=1,3
d_t_old(j,i)=d_t(j,i)
enddo
enddo
endif
- if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
+ if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
call inertia_tensor
call vcm_vel(vcm)
do j=1,3
call kinetic(EK)
kinetic_T=2.0d0/(dimen3*Rb)*EK
scalfac=dsqrt(T_bath/kinetic_T)
- write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
+ write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
do i=0,2*nres
do j=1,3
d_t_old(j,i)=scalfac*d_t(j,i)
stop
#endif
endif
+ itime_mat=itime
if (ntwe.ne.0) then
if (mod(itime,ntwe).eq.0) then
call statout(itime)
write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
& (d_t(j,i+nres),j=1,3)
enddo
+ endif
+ if (lang.eq.0) then
c Zeroing the total angular momentum of the system
write(iout,*) "Calling the zero-angular
& momentum subroutine"
- endif
call inertia_tensor
c Getting the potential energy and forces and velocities and accelerations
call vcm_vel(vcm)
if(me.eq.king.or..not.out1file)then
write (iout,*) "vcm right after adjustment:"
write (iout,*) (vcm(j),j=1,3)
+ write (iout,*) "Initial velocities after adjustment"
+ do i=0,nres
+ write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
+ & (d_t(j,i+nres),j=1,3)
+ enddo
+ call flush(iout)
+ endif
endif
- call flush(iout)
write (iout,*) "init_MD before initial structure REST ",rest
if (.not.rest) then
122 continue
c write (iout,*) "PotE-homology",potE
totE=EK+potE
itime=0
+ itime_mat=itime
if (ntwe.ne.0) call statout(itime)
if(me.eq.king.or..not.out1file)
& write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
!
#define WLOS
#ifdef WLOS
+ if (nnt.eq.1) then
+ d_t(:,0)=d_t(:,1)
+ endif
do i=1,nres
if (itype(i).eq.10 .or. itype(i).eq.ntyp1) then
do j=1,3
enddo
end if
enddo
+ d_t(:,nres)=0.0d0
d_t(:,nct)=0.0d0
+ d_t(:,2*nres)=0.0d0
+ if (nnt.gt.1) then
+ d_t(:,0)=d_t(:,1)
+ d_t(:,1)=0.0d0
+ endif
c d_a(:,0)=d_a(:,1)
c d_a(:,1)=0.0d0
c write (iout,*) "Shifting accelerations"
- do ichain=1,nchain
+ do ichain=2,nchain
c write (iout,*) "ichain",chain_border1(1,ichain)-1,
c & chain_border1(1,ichain)
d_t(:,chain_border1(1,ichain)-1)=d_t(:,chain_border1(1,ichain))
matmult.o readrtns_CSA.o parmread.o gen_rand_conf.o printmat.o map.o \
pinorm.o randgens.o rescode.o intcor.o timing.o misc.o \
cart2intgrad.o checkder_p.o contact_cp econstr_local.o econstr_qlike.o \
- econstrq-PMF.o PMFprocess.o energy_p_new_barrier.o \
+ econstrq-PMF.o PMFprocess.o energy_p_new_barrier.o make_xx_list \
energy_p_new-sep_barrier.o gradient_p.o minimize_p.o sumsld.o \
cored.o rmdd.o geomout.o readpdb.o regularize.o thread.o fitsq.o mcm.o \
mc.o bond_move.o refsys.o check_sc_distr.o check_bond.o contact.o \
NEWCORR5D: CPPFLAGS = -DPROCOR -DLINUX -DPGI -DAMD64 -DUNRES -DISNAN -DMP -DMPI \
-DSPLITELE -DLANG0 -DNEWCORR -DCORRCD -DFIVEDIAG -DLBFGS #-DMYGAUSS #-DTIMING
-NEWCORR5D: BIN = ~/bin/unres_ifort_MPICH-okeanos_SC-HCD5.exe
+NEWCORR5D: BIN = ~/bin/unres_ifort_MPICH-okeanos_SC-HCD5-l.exe
NEWCORR5D: ${object_lbfgs} ${object} fdisy.o fdiag.o machpd.o kinetic_CASC.o xdrf/libxdrf.a
gcc -o compinfo compinfo.c
./compinfo | true
include 'COMMON.SPLITELE'
include 'COMMON.TORCNSTR'
include 'COMMON.SAXS'
+ include 'COMMON.MD'
double precision evdw,evdw1,evdw2,evdw2_14,ees,eel_loc,
& eello_turn3,eello_turn4,edfadis,estr,ehpb,ebe,ethetacnstr,
& escloc,etors,edihcnstr,etors_d,esccor,ecorr,ecorr5,ecorr6,eturn6,
edfanei=0.0d0
edfabet=0.0d0
#endif
+ if (nfgtasks.gt.1) then
+ call MPI_Bcast(itime_mat,1,MPI_INT,king,FG_COMM,IERROR)
+ endif
+ if (mod(itime_mat,imatupdate).eq.0) call make_SCp_inter_list
+ if (mod(itime_mat,imatupdate).eq.0) call make_SCSC_inter_list
+ if (mod(itime_mat,imatupdate).eq.0) call make_pp_inter_list
c print *,'Processor',myrank,' calling etotal ipot=',ipot
c print *,'Processor',myrank,' nnt=',nnt,' nct=',nct
#else
#endif
double precision gg(3)
double precision evdw,evdwij
- integer i,j,k,itypi,itypj,itypi1,num_conti,iint
+ integer i,j,k,itypi,itypj,itypi1,num_conti,iint,icont
double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij,
& sigij,r0ij,rcut,sqrij,sss1,sssgrad1
double precision fcont,fprimcont
double precision sscale,sscagrad
c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
evdw=0.0D0
- do i=iatsc_s,iatsc_e
+c do i=iatsc_s,iatsc_e
+ do icont=g_listscsc_start,g_listscsc_end
+ i=newcontlisti(icont)
+ j=newcontlistj(icont)
itypi=iabs(itype(i))
if (itypi.eq.ntyp1) cycle
itypi1=iabs(itype(i+1))
C
C Calculate SC interaction energy.
C
- do iint=1,nint_gr(i)
+c do iint=1,nint_gr(i)
cd write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
cd & 'iend=',iend(i,iint)
- do j=istart(i,iint),iend(i,iint)
+c do j=istart(i,iint),iend(i,iint)
itypj=iabs(itype(j))
if (itypj.eq.ntyp1) cycle
xj=c(1,nres+j)-xi
endif
endif
#endif
- enddo ! j
- enddo ! iint
+c enddo ! j
+c enddo ! iint
C Change 12/1/95
#ifdef FOURBODY
num_cont(i)=num_conti
include 'COMMON.SPLITELE'
double precision gg(3)
double precision evdw,evdwij
- integer i,j,k,itypi,itypj,itypi1,iint
+ integer i,j,k,itypi,itypj,itypi1,iint,icont
double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij,
& fac_augm,e_augm,r_inv_ij,r_shift_inv,sss1,sssgrad1
logical scheck
double precision sscale,sscagrad
c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
evdw=0.0D0
- do i=iatsc_s,iatsc_e
+c do i=iatsc_s,iatsc_e
+ do icont=g_listscsc_start,g_listscsc_end
+ i=newcontlisti(icont)
+ j=newcontlistj(icont)
itypi=iabs(itype(i))
if (itypi.eq.ntyp1) cycle
itypi1=iabs(itype(i+1))
C
C Calculate SC interaction energy.
C
- do iint=1,nint_gr(i)
- do j=istart(i,iint),iend(i,iint)
+c do iint=1,nint_gr(i)
+c do j=istart(i,iint),iend(i,iint)
itypj=iabs(itype(j))
if (itypj.eq.ntyp1) cycle
xj=c(1,nres+j)-xi
cgrad gvdwc(l,k)=gvdwc(l,k)+gg(l)
cgrad enddo
cgrad enddo
- enddo ! j
- enddo ! iint
+c enddo ! j
+c enddo ! iint
enddo ! i
do i=1,nct
do j=1,3
integer icall
common /srutu/ icall
double precision evdw
- integer itypi,itypj,itypi1,iint,ind
+ integer itypi,itypj,itypi1,iint,ind,icont
double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi,
& sss1,sssgrad1
double precision sscale,sscagrad
lprn=.false.
c endif
ind=0
- do i=iatsc_s,iatsc_e
+c do i=iatsc_s,iatsc_e
+ do icont=g_listscsc_start,g_listscsc_end
+ i=newcontlisti(icont)
+ j=newcontlistj(icont)
itypi=iabs(itype(i))
if (itypi.eq.ntyp1) cycle
itypi1=iabs(itype(i+1))
C
C Calculate SC interaction energy.
C
- do iint=1,nint_gr(i)
- do j=istart(i,iint),iend(i,iint)
+c do iint=1,nint_gr(i)
+c do j=istart(i,iint),iend(i,iint)
ind=ind+1
itypj=iabs(itype(j))
if (itypj.eq.ntyp1) cycle
C Calculate the angular part of the gradient and sum add the contributions
C to the appropriate components of the Cartesian gradient.
call sc_grad
- enddo ! j
- enddo ! iint
+! enddo ! j
+! enddo ! iint
enddo ! i
c stop
return
logical lprn
integer xshift,yshift,zshift,subchap
double precision evdw
- integer itypi,itypj,itypi1,iint,ind
+ integer itypi,itypj,itypi1,iint,ind,icont
double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi
double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij,
& sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe,
C do xshift=-1,1
C do yshift=-1,1
C do zshift=-1,1
- do i=iatsc_s,iatsc_e
+c do i=iatsc_s,iatsc_e
+ do icont=g_listscsc_start,g_listscsc_end
+ i=newcontlisti(icont)
+ j=newcontlistj(icont)
itypi=iabs(itype(i))
if (itypi.eq.ntyp1) cycle
itypi1=iabs(itype(i+1))
C
C Calculate SC interaction energy.
C
- do iint=1,nint_gr(i)
- do j=istart(i,iint),iend(i,iint)
+c do iint=1,nint_gr(i)
+c do j=istart(i,iint),iend(i,iint)
IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN
c write(iout,*) "PRZED ZWYKLE", evdwij
c call sc_grad_scale(sss)
call sc_grad
ENDIF ! dyn_ss
- enddo ! j
- enddo ! iint
+c enddo ! j
+c enddo ! iint
enddo ! i
C enddo ! zshift
C enddo ! yshift
common /srutu/ icall
logical lprn
double precision evdw
- integer itypi,itypj,itypi1,iint,ind
+ integer itypi,itypj,itypi1,iint,ind,icont
double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,r0ij,
& xi,yi,zi,fac_augm,e_augm
double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij,
lprn=.false.
c if (icall.eq.0) lprn=.true.
ind=0
- do i=iatsc_s,iatsc_e
+c do i=iatsc_s,iatsc_e
+ do icont=g_listscsc_start,g_listscsc_end
+ i=newcontlisti(icont)
+ j=newcontlistj(icont)
itypi=iabs(itype(i))
if (itypi.eq.ntyp1) cycle
itypi1=iabs(itype(i+1))
C
C Calculate SC interaction energy.
C
- do iint=1,nint_gr(i)
- do j=istart(i,iint),iend(i,iint)
+c do iint=1,nint_gr(i)
+c do j=istart(i,iint),iend(i,iint)
ind=ind+1
itypj=iabs(itype(j))
if (itypj.eq.ntyp1) cycle
C Calculate angular part of the gradient.
c call sc_grad_scale(sss)
call sc_grad
- enddo ! j
- enddo ! iint
+c enddo ! j
+c enddo ! iint
enddo ! i
end
C-----------------------------------------------------------------------------
dimension gg(3)
cd print *,'Entering Esoft_sphere nnt=',nnt,' nct=',nct
evdw=0.0D0
- do i=iatsc_s,iatsc_e
+c do i=iatsc_s,iatsc_e
+ do icont=g_listscsc_start,g_listscsc_end
+ i=newcontlisti(icont)
+ j=newcontlistj(icont)
itypi=iabs(itype(i))
if (itypi.eq.ntyp1) cycle
itypi1=iabs(itype(i+1))
C
C Calculate SC interaction energy.
C
- do iint=1,nint_gr(i)
+c do iint=1,nint_gr(i)
cd write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
cd & 'iend=',iend(i,iint)
- do j=istart(i,iint),iend(i,iint)
+c do j=istart(i,iint),iend(i,iint)
itypj=iabs(itype(j))
if (itypj.eq.ntyp1) cycle
xj=c(1,nres+j)-xi
cgrad gvdwc(l,k)=gvdwc(l,k)+gg(l)
cgrad enddo
cgrad enddo
- enddo ! j
- enddo ! iint
+c enddo ! j
+c enddo ! iint
enddo ! i
return
end
c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
c
CTU KURWA
- do i=iatel_s,iatel_e
+c do i=iatel_s,iatel_e
+ do icont=g_listpp_start,g_listpp_end
+ i=newcontlistppi(icont)
+ j=newcontlistppj(icont)
C do i=75,75
c if (i.le.1) cycle
if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
num_conti=num_cont_hb(i)
#endif
C I TU KURWA
- do j=ielstart(i),ielend(i)
+c do j=ielstart(i),ielend(i)
C do j=16,17
C write (iout,*) i,j
C if (j.le.1) cycle
c & .or.itype(j-1).eq.ntyp1
&) cycle
call eelecij(i,j,ees,evdw1,eel_loc)
- enddo ! j
+c enddo ! j
#ifdef FOURBODY
num_cont_hb(i)=num_conti
#endif
C do xshift=-1,1
C do yshift=-1,1
C do zshift=-1,1
- do i=iatscp_s,iatscp_e
+c do i=iatscp_s,iatscp_e
+ do icont=g_listscp_start,g_listscp_end
+ i=newcontlistscpi(icont)
+ j=newcontlistscpj(icont)
if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
iteli=itel(i)
xi=0.5D0*(c(1,i)+c(1,i+1))
C xi=xi+xshift*boxxsize
C yi=yi+yshift*boxysize
C zi=zi+zshift*boxzsize
- do iint=1,nscp_gr(i)
+c do iint=1,nscp_gr(i)
- do j=iscpstart(i,iint),iscpend(i,iint)
+c do j=iscpstart(i,iint),iscpend(i,iint)
if (itype(j).eq.ntyp1) cycle
itypj=iabs(itype(j))
C Uncomment following three lines for SC-p interactions
gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
enddo
- enddo
+c enddo
- enddo ! iint
+c enddo ! iint
enddo ! i
C enddo !zshift
C enddo !yshift
include 'COMMON.SPLITELE'
integer xshift,yshift,zshift
double precision ggg(3)
- integer i,iint,j,k,iteli,itypj,subchap
+ integer i,iint,j,k,iteli,itypj,subchap,icont
double precision xi,yi,zi,xj,yj,zj,rrij,sss1,sssgrad1,
& fac,e1,e2,rij
double precision evdw2,evdw2_14,evdwij
C do yshift=-1,1
C do zshift=-1,1
if (energy_dec) write (iout,*) "escp:",r_cut_int,rlamb
- do i=iatscp_s,iatscp_e
+c do i=iatscp_s,iatscp_e
+ do icont=g_listscp_start,g_listscp_end
+ i=newcontlistscpi(icont)
+ j=newcontlistscpj(icont)
if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
iteli=itel(i)
xi=0.5D0*(c(1,i)+c(1,i+1))
c & (zi.lt.((zshift-0.5d0)*boxzsize))) then
c go to 136
c endif
- do iint=1,nscp_gr(i)
+c do iint=1,nscp_gr(i)
- do j=iscpstart(i,iint),iscpend(i,iint)
+c do j=iscpstart(i,iint),iscpend(i,iint)
itypj=iabs(itype(j))
if (itypj.eq.ntyp1) cycle
C Uncomment following three lines for SC-p interactions
gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
enddo
c endif !endif for sscale cutoff
- enddo ! j
+c enddo ! j
- enddo ! iint
+c enddo ! iint
enddo ! i
c enddo !zshift
c enddo !yshift
weights_(17)=wbond
weights_(18)=scal14
weights_(21)=wsccor
- weights_(22)=wtube
+ weights_(22)=wliptran
+ weights_(25)=wtube
weights_(26)=wsaxs
weights_(28)=wdfa_dist
weights_(29)=wdfa_tor
wbond=weights(17)
scal14=weights(18)
wsccor=weights(21)
- wtube=weights(22)
+ wliptran=weights(22)
+ wtube=weights(25)
wsaxs=weights(26)
wdfa_dist=weights_(28)
wdfa_tor=weights_(29)
C print *,"przed lipidami"
if (wliptran.gt.0) then
call Eliptransfer(eliptran)
+ else
+ eliptran=0.0d0
endif
C print *,"za lipidami"
if (AFMlog.gt.0) then
include 'COMMON.SBRIDGE'
include 'COMMON.NAMES'
include 'COMMON.IOUNITS'
+ include 'COMMON.SPLITELE'
#ifdef FOURBODY
include 'COMMON.CONTACTS'
include 'COMMON.CONTMAT'
double precision evdw,evdwij
integer i,j,k,itypi,itypj,itypi1,num_conti,iint
double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij,
- & sigij,r0ij,rcut
+ & sigij,r0ij,rcut,sqrij,sss1,sssgrad1
double precision fcont,fprimcont
+ double precision sscale,sscagrad
c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
evdw=0.0D0
do i=iatsc_s,iatsc_e
C Change 12/1/95 to calculate four-body interactions
rij=xj*xj+yj*yj+zj*zj
rrij=1.0D0/rij
+ sqrij=dsqrt(rij)
+ sss1=sscale(sqrij,r_cut_int)
+ if (sss1.eq.0.0d0) cycle
+ sssgrad1=sscagrad(sqrij,r_cut_int)
+
c write (iout,*)'i=',i,' j=',j,' itypi=',itypi,' itypj=',itypj
eps0ij=eps(itypi,itypj)
fac=rrij**expon2
cd & restyp(itypi),i,restyp(itypj),j,a(itypi,itypj),
cd & bb(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,epsi,sigm,
cd & (c(k,i),k=1,3),(c(k,j),k=1,3)
- evdw=evdw+evdwij
+ evdw=evdw+sss1*evdwij
C
C Calculate the components of the gradient in DC and X
C
- fac=-rrij*(e1+evdwij)
+ fac=-rrij*(e1+evdwij)*sss1
+ & +evdwij*sssgrad1/sqrij/expon
gg(1)=xj*fac
gg(2)=yj*fac
gg(3)=zj*fac
include 'COMMON.INTERACT'
include 'COMMON.IOUNITS'
include 'COMMON.NAMES'
+ include 'COMMON.SPLITELE'
double precision gg(3)
double precision evdw,evdwij
integer i,j,k,itypi,itypj,itypi1,iint
double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij,
- & fac_augm,e_augm,r_inv_ij,r_shift_inv
+ & fac_augm,e_augm,r_inv_ij,r_shift_inv,sss1,sssgrad1
logical scheck
+ double precision sscale,sscagrad
c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
evdw=0.0D0
do i=iatsc_s,iatsc_e
e_augm=augm(itypi,itypj)*fac_augm
r_inv_ij=dsqrt(rrij)
rij=1.0D0/r_inv_ij
+ sss1=sscale(rij,r_cut_int)
+ if (sss1.eq.0.0d0) cycle
+ sssgrad1=sscagrad(rij,r_cut_int)
r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
fac=r_shift_inv**expon
C have you changed here?
cd & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
cd & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
cd & (c(k,i),k=1,3),(c(k,j),k=1,3)
- evdw=evdw+evdwij
+ evdw=evdw+evdwij*sss1
C
C Calculate the components of the gradient in DC and X
C
fac=-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2)
+ & +evdwij*sssgrad1*r_inv_ij/expon
gg(1)=xj*fac
gg(2)=yj*fac
gg(3)=zj*fac
include 'COMMON.INTERACT'
include 'COMMON.IOUNITS'
include 'COMMON.CALC'
+ include 'COMMON.SPLITELE'
integer icall
common /srutu/ icall
double precision evdw
integer itypi,itypj,itypi1,iint,ind
- double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi
+ double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi,
+ & sss1,sssgrad1
+ double precision sscale,sscagrad
c double precision rrsave(maxdim)
logical lprn
evdw=0.0D0
cd rrij=rrsave(ind)
cd endif
rij=dsqrt(rrij)
+ sss1=sscale(1.0d0/rij,r_cut_int)
+ if (sss1.eq.0.0d0) cycle
+ sssgrad1=sscagrad(1.0d0/rij,r_cut_int)
C Calculate the angle-dependent terms of energy & contributions to derivatives.
call sc_angular
C Calculate whole angle-dependent part of epsilon and contributions
eps2der=evdwij*eps3rt
eps3der=evdwij*eps2rt
evdwij=evdwij*eps2rt*eps3rt
- evdw=evdw+evdwij
+ evdw=evdw+sss1*evdwij
if (lprn) then
sigm=dabs(aa/bb)**(1.0D0/6.0D0)
epsi=bb**2/aa
fac=-expon*(e1+evdwij)
sigder=fac/sigsq
fac=rrij*fac
+ & +evdwij*sssgrad1/sss1*rij
C Calculate radial part of the gradient
gg(1)=xj*fac
gg(2)=yj*fac
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
rij=dsqrt(rrij)
sss=sscale(1.0d0/rij,r_cut_int)
- sssgrad=sscagrad(1.0d0/rij,r_cut_int)
-
c write (iout,'(a7,4f8.3)')
c & "ssscale",sss,((1.0d0/rij)/sigma(itypi,itypj)),r_cut,rlamb
- if (sss.gt.0.0d0) then
+ if (sss.eq.0.0d0) cycle
+ sssgrad=sscagrad(1.0d0/rij,r_cut_int)
C Calculate angle-dependent terms of energy and contributions to their
C derivatives.
call sc_angular
& evdwij
endif
- if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
- & 'evdw',i,j,evdwij
+ if (energy_dec) write (iout,'(a,2i5,3f10.5)')
+ & 'r sss evdw',i,j,rij,sss,evdwij
C Calculate gradient components.
e1=e1*eps1*eps2rt**2*eps3rt**2
fac=rij*fac
c print '(2i4,6f8.4)',i,j,sss,sssgrad*
c & evdwij,fac,sigma(itypi,itypj),expon
- fac=fac+evdwij/sss*sssgrad*rij
+ fac=fac+evdwij*sssgrad/sss*rij
c fac=0.0d0
C Calculate the radial part of the gradient
gg_lipi(3)=eps1*(eps2rt*eps2rt)
- &*(eps3rt*eps3rt)*sss/2.0d0*(faclip*faclip*
- & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj))
- &+faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj)))
+ & *(eps3rt*eps3rt)*sss/2.0d0*(faclip*faclip*
+ & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj))
+ & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj)))
gg_lipj(3)=ssgradlipj*gg_lipi(3)
gg_lipi(3)=gg_lipi(3)*ssgradlipi
C gg_lipi(3)=0.0d0
gg(2)=yj*fac
gg(3)=zj*fac
C Calculate angular part of the gradient.
+c call sc_grad_scale(sss)
call sc_grad
- endif
ENDIF ! dyn_ss
enddo ! j
enddo ! iint
include 'COMMON.INTERACT'
include 'COMMON.IOUNITS'
include 'COMMON.CALC'
+ include 'COMMON.SPLITELE'
integer xshift,yshift,zshift,subchap
integer icall
common /srutu/ icall
& xi,yi,zi,fac_augm,e_augm
double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij,
& sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe,
- & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip
+ & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip,sssgrad1
double precision dist,sscale,sscagrad,sscagradlip,sscalelip
evdw=0.0D0
c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
dzj=dc_norm(3,nres+j)
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
rij=dsqrt(rrij)
+ sss=sscale(1.0d0/rij,r_cut_int)
+ if (sss.eq.0.0d0) cycle
+ sssgrad=sscagrad(1.0d0/rij,r_cut_int)
C Calculate angle-dependent terms of energy and contributions to their
C derivatives.
call sc_angular
fac=-expon*(e1+evdwij)*rij_shift
sigder=fac*sigder
fac=rij*fac-2*expon*rrij*e_augm
- fac=fac+evdwij/sss*sssgrad/sigma(itypi,itypj)*rij
+ fac=fac+(evdwij+e_augm)*sssgrad/sss*rij
C Calculate the radial part of the gradient
gg(1)=xj*fac
gg(2)=yj*fac
gg(3)=zj*fac
C Calculate angular part of the gradient.
+c call sc_grad_scale(sss)
call sc_grad
enddo ! j
enddo ! iint
end
C-------------------------------------------------------------------------------
subroutine eelecij(i,j,ees,evdw1,eel_loc)
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
#ifdef MPI
include "mpif.h"
include 'COMMON.TIME1'
include 'COMMON.SPLITELE'
include 'COMMON.SHIELD'
- dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
+ double precision ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
& erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
& aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4),gmuij1(4),gmuji1(4),
& gmuij2(4),gmuji2(4)
+ double precision dxi,dyi,dzi
+ double precision dx_normi,dy_normi,dz_normi,aux
+ integer j1,j2,lll,num_conti
common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
& dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
& num_conti,j1,j2
+ integer k,i,j,iteli,itelj,kkk,l,kkll,m,isubchap,ilist,iresshield
+ double precision ael6i,rrmij,rmij,r0ij,fcont,fprimcont,ees0tmp
+ double precision ees,evdw1,eel_loc,aaa,bbb,ael3i
+ double precision dxj,dyj,dzj,dx_normj,dy_normj,dz_normj,xj,yj,zj,
+ & rij,r3ij,r6ij,cosa,cosb,cosg,fac,ev1,ev2,fac3,fac4,
+ & evdwij,el1,el2,eesij,ees0ij,facvdw,facel,fac1,ecosa,
+ & ecosb,ecosg,ury,urz,vry,vrz,facr,a22der,a23der,a32der,
+ & a33der,eel_loc_ij,cosa4,wij,cosbg1,cosbg2,ees0pij,
+ & ees0pij1,ees0mij,ees0mij1,fac3p,ees0mijp,ees0pijp,
+ & ecosa1,ecosb1,ecosg1,ecosa2,ecosb2,ecosg2,ecosap,ecosbp,
+ & ecosgp,ecosam,ecosbm,ecosgm,ghalf,rlocshield
+ double precision a22,a23,a32,a33,geel_loc_ij,geel_loc_ji
+ double precision dist_init,xj_safe,yj_safe,zj_safe,
+ & xj_temp,yj_temp,zj_temp,dist_temp,xmedi,ymedi,zmedi
+ double precision sscale,sscagrad,scalar
+
c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
#ifdef MOMENT
double precision scal_el /1.0d0/
C zj=zj-zmedi
rij=xj*xj+yj*yj+zj*zj
- sss=sscale(sqrt(rij),r_cut_int)
- sssgrad=sscagrad(sqrt(rij),r_cut_int)
+ sss=sscale(dsqrt(rij),r_cut_int)
+ if (sss.eq.0.0d0) return
+ sssgrad=sscagrad(dsqrt(rij),r_cut_int)
c if (sss.gt.0.0d0) then
rrmij=1.0D0/rij
rij=dsqrt(rij)
fac_shield(i)=1.0
fac_shield(j)=1.0
eesij=(el1+el2)
- ees=ees+eesij
+ ees=ees+eesij*sss
endif
evdw1=evdw1+evdwij*sss
cd write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
cd & xmedi,ymedi,zmedi,xj,yj,zj
if (energy_dec) then
- write (iout,'(a6,2i5,0pf7.3,2i5,3e11.3)')
- &'evdw1',i,j,evdwij
- &,iteli,itelj,aaa,evdw1,sss
- write (iout,'(a6,2i5,0pf7.3,2f8.3)') 'ees',i,j,eesij,
- &fac_shield(i),fac_shield(j)
+ write (iout,'(a6,2i5,0pf7.3,2i5,e11.3,3f10.5)')
+ & 'evdw1',i,j,evdwij,iteli,itelj,aaa,evdw1,sss,rij
+ write (iout,'(a6,2i5,0pf7.3,2f8.3)') 'ees',i,j,eesij,
+ & fac_shield(i),fac_shield(j)
endif
C
*
* Radial derivatives. First process both termini of the fragment (i,j)
*
- ggg(1)=facel*xj
- ggg(2)=facel*yj
- ggg(3)=facel*zj
+ aux=facel*sss+rmij*sssgrad*eesij
+ ggg(1)=aux*xj
+ ggg(2)=aux*yj
+ ggg(3)=aux*zj
if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
& (shield_mode.gt.0)) then
C print *,i,j
iresshield=shield_list(ilist,j)
do k=1,3
rlocshield=grad_shield_side(k,ilist,j)*eesij/fac_shield(j)
- & *2.0
+ & *2.0*sss
gshieldx(k,iresshield)=gshieldx(k,iresshield)+
& rlocshield
- & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)*2.0
+ & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)*2.0*sss
gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
C & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)
do k=1,3
gshieldc(k,i)=gshieldc(k,i)+
- & grad_shield(k,i)*eesij/fac_shield(i)*2.0
+ & grad_shield(k,i)*eesij/fac_shield(i)*2.0*sss
gshieldc(k,j)=gshieldc(k,j)+
- & grad_shield(k,j)*eesij/fac_shield(j)*2.0
+ & grad_shield(k,j)*eesij/fac_shield(j)*2.0*sss
gshieldc(k,i-1)=gshieldc(k,i-1)+
- & grad_shield(k,i)*eesij/fac_shield(i)*2.0
+ & grad_shield(k,i)*eesij/fac_shield(i)*2.0*sss
gshieldc(k,j-1)=gshieldc(k,j-1)+
- & grad_shield(k,j)*eesij/fac_shield(j)*2.0
+ & grad_shield(k,j)*eesij/fac_shield(j)*2.0*sss
enddo
endif
cgrad gelc(l,k)=gelc(l,k)+ggg(l)
cgrad enddo
cgrad enddo
- if (sss.gt.0.0) then
- ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj
- ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj
- ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj
- else
- ggg(1)=0.0
- ggg(2)=0.0
- ggg(3)=0.0
- endif
+ facvdw=facvdw+sssgrad*rmij*evdwij
+ ggg(1)=facvdw*xj
+ ggg(2)=facvdw*yj
+ ggg(3)=facvdw*zj
c do k=1,3
c ghalf=0.5D0*ggg(k)
c gvdwpp(k,i)=gvdwpp(k,i)+ghalf
cgrad enddo
#else
C MARYSIA
- facvdw=(ev1+evdwij)*sss
+ facvdw=(ev1+evdwij)
facel=(el1+eesij)
fac1=fac
- fac=-3*rrmij*(facvdw+facvdw+facel)
+ fac=-3*rrmij*(facvdw+facvdw+facel)*sss
+ & +(evdwij+eesij)*sssgrad*rrmij
erij(1)=xj*rmij
erij(2)=yj*rmij
erij(3)=zj*rmij
cd & (dcosg(k),k=1,3)
do k=1,3
ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))*
- & fac_shield(i)**2*fac_shield(j)**2
+ & fac_shield(i)**2*fac_shield(j)**2*sss
enddo
c do k=1,3
c ghalf=0.5D0*ggg(k)
do k=1,3
gelc(k,i)=gelc(k,i)
& +((ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
- & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1))
+ & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1))*sss
& *fac_shield(i)**2*fac_shield(j)**2
gelc(k,j)=gelc(k,j)
& +((ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
- & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1))
+ & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1))*sss
& *fac_shield(i)**2*fac_shield(j)**2
gelc_long(k,j)=gelc_long(k,j)+ggg(k)
gelc_long(k,i)=gelc_long(k,i)-ggg(k)
C fac_shield(j)=0.6
endif
eel_loc_ij=eel_loc_ij
- & *fac_shield(i)*fac_shield(j)
+ & *fac_shield(i)*fac_shield(j)*sss
c if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
c & 'eelloc',i,j,eel_loc_ij
C Now derivative over eel_loc
& +a23*gmuij1(2)
& +a32*gmuij1(3)
& +a33*gmuij1(4))
- & *fac_shield(i)*fac_shield(j)
+ & *fac_shield(i)*fac_shield(j)*sss
c write(iout,*) "derivative over thatai"
c write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3),
c & a33*gmuij1(4)
& +a33*gmuij2(4)
gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
& geel_loc_ij*wel_loc
- & *fac_shield(i)*fac_shield(j)
+ & *fac_shield(i)*fac_shield(j)*sss
c Derivative over j residue
geel_loc_ji=a22*gmuji1(1)
gloc(nphi+j,icg)=gloc(nphi+j,icg)+
& geel_loc_ji*wel_loc
- & *fac_shield(i)*fac_shield(j)
+ & *fac_shield(i)*fac_shield(j)*sss
geel_loc_ji=
& +a22*gmuji2(1)
c & a33*gmuji2(4)
gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+
& geel_loc_ji*wel_loc
- & *fac_shield(i)*fac_shield(j)
+ & *fac_shield(i)*fac_shield(j)*sss
#endif
cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
& gel_loc_loc(i-1)=gel_loc_loc(i-1)+
& (a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j)
& +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j))
- & *fac_shield(i)*fac_shield(j)
+ & *fac_shield(i)*fac_shield(j)*sss
gel_loc_loc(j-1)=gel_loc_loc(j-1)+
& (a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j)
& +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j))
- & *fac_shield(i)*fac_shield(j)
+ & *fac_shield(i)*fac_shield(j)*sss
C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2)
+ aux=eel_loc_ij/sss*sssgrad*rmij
+ ggg(1)=aux*xj
+ ggg(2)=aux*yj
+ ggg(3)=aux*zj
do l=1,3
- ggg(l)=(agg(l,1)*muij(1)+
+ ggg(l)=ggg(l)+(agg(l,1)*muij(1)+
& agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4))
- & *fac_shield(i)*fac_shield(j)
+ & *fac_shield(i)*fac_shield(j)*sss
gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l)
gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l)
cgrad ghalf=0.5d0*ggg(l)
do l=1,3
gel_loc(l,i)=gel_loc(l,i)+(aggi(l,1)*muij(1)+
& aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4))
- & *fac_shield(i)*fac_shield(j)
+ & *fac_shield(i)*fac_shield(j)*sss
gel_loc(l,i+1)=gel_loc(l,i+1)+(aggi1(l,1)*muij(1)+
& aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4))
- & *fac_shield(i)*fac_shield(j)
+ & *fac_shield(i)*fac_shield(j)*sss
gel_loc(l,j)=gel_loc(l,j)+(aggj(l,1)*muij(1)+
& aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4))
- & *fac_shield(i)*fac_shield(j)
+ & *fac_shield(i)*fac_shield(j)*sss
gel_loc(l,j1)=gel_loc(l,j1)+(aggj1(l,1)*muij(1)+
& aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4))
- & *fac_shield(i)*fac_shield(j)
+ & *fac_shield(i)*fac_shield(j)*sss
enddo
ENDIF
C fac_shield(j)=0.6d0
endif
ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij)
- & *fac_shield(i)*fac_shield(j)
+ & *fac_shield(i)*fac_shield(j)*sss
ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij)
- & *fac_shield(i)*fac_shield(j)
+ & *fac_shield(i)*fac_shield(j)*sss
C Diagnostics. Comment out or remove after debugging!
c ees0p(num_conti,i)=0.5D0*fac3*ees0pij
c ees0m(num_conti,i)=0.5D0*fac3*ees0mij
gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k)
enddo
gggp(1)=gggp(1)+ees0pijp*xj
+ & +ees0p(num_conti,i)/sss*rmij*xj*sssgrad
gggp(2)=gggp(2)+ees0pijp*yj
+ & +ees0p(num_conti,i)/sss*rmij*yj*sssgrad
gggp(3)=gggp(3)+ees0pijp*zj
+ & +ees0p(num_conti,i)/sss*rmij*zj*sssgrad
gggm(1)=gggm(1)+ees0mijp*xj
+ & +ees0m(num_conti,i)/sss*rmij*xj*sssgrad
gggm(2)=gggm(2)+ees0mijp*yj
+ & +ees0m(num_conti,i)/sss*rmij*yj*sssgrad
gggm(3)=gggm(3)+ees0mijp*zj
+ & +ees0m(num_conti,i)/sss*rmij*zj*sssgrad
C Derivatives due to the contact function
gacont_hbr(1,num_conti,i)=fprimcont*xj
gacont_hbr(2,num_conti,i)=fprimcont*yj
gacontp_hb1(k,num_conti,i)=!ghalfp
& +(ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i))
& + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
- & *fac_shield(i)*fac_shield(j)
+ & *fac_shield(i)*fac_shield(j)*sss
gacontp_hb2(k,num_conti,i)=!ghalfp
& +(ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j))
& + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
- & *fac_shield(i)*fac_shield(j)
+ & *fac_shield(i)*fac_shield(j)*sss
gacontp_hb3(k,num_conti,i)=gggp(k)
- & *fac_shield(i)*fac_shield(j)
+ & *fac_shield(i)*fac_shield(j)*sss
gacontm_hb1(k,num_conti,i)=!ghalfm
& +(ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i))
& + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
- & *fac_shield(i)*fac_shield(j)
+ & *fac_shield(i)*fac_shield(j)*sss
gacontm_hb2(k,num_conti,i)=!ghalfm
& +(ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j))
& + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
- & *fac_shield(i)*fac_shield(j)
+ & *fac_shield(i)*fac_shield(j)*sss
gacontm_hb3(k,num_conti,i)=gggm(k)
- & *fac_shield(i)*fac_shield(j)
+ & *fac_shield(i)*fac_shield(j)*sss
enddo
C Diagnostics. Comment out or remove after debugging!
C peptide-group centers and side chains and its gradient in virtual-bond and
C side-chain vectors.
C
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
include 'COMMON.GEO'
include 'COMMON.VAR'
include 'COMMON.CONTROL'
include 'COMMON.SPLITELE'
integer xshift,yshift,zshift
- dimension ggg(3)
+ double precision ggg(3)
+ integer i,iint,j,k,iteli,itypj,subchap
+ double precision xi,yi,zi,xj,yj,zj,rrij,sss1,sssgrad1,
+ & fac,e1,e2,rij
+ double precision evdw2,evdw2_14,evdwij
+ double precision xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,
+ & dist_temp, dist_init
+ double precision sscale,sscagrad
evdw2=0.0D0
evdw2_14=0.0d0
c print *,boxxsize,boxysize,boxzsize,'wymiary pudla'
endif
evdwij=e1+e2
evdw2=evdw2+evdwij*sss
- if (energy_dec) write (iout,'(a6,2i5,0pf7.3,2i3,3e11.3)')
- & 'evdw2',i,j,evdwij,iteli,itypj,fac,aad(itypj,iteli),
+ if (energy_dec) write (iout,'(a6,2i5,3f7.3,2i3,3e11.3)')
+ & 'evdw2',i,j,1.0d0/dsqrt(rrij),sss,
+ & evdwij,iteli,itypj,fac,aad(itypj,iteli),
& bad(itypj,iteli)
C
C Calculate contributions to the gradient in the virtual-bond and SC vectors.
include 'COMMON.DERIV'
include 'COMMON.SPLITELE'
include 'COMMON.VAR'
+ include 'COMMON.MD'
c Common blocks from the diagonalization routines
integer IR,IW,IP,IJK,IPK,IDAF,NAV,IODA,KDIAG,ICORFL,IXDR
integer i,idumm,j,k,l,ichir1,ichir2,iblock,m
call proc_proc(rr,i)
#endif
#endif
-
+ itime_mat=0.
kdiag=0
icorfl=0
iw=2
include 'COMMON.IOUNITS'
include 'COMMON.LAGRANGE.5diag'
include 'COMMON.INTERACT'
+ include 'COMMON.VAR'
integer ndim
double precision forces(3*ndim),accel(3,0:maxres2),rs(ndim),
& xsolv(ndim),d_a_vec(6*nres)
innt=chain_border(1,ichain)
inct=chain_border(2,ichain)
do i=iposc,iposc+n-1
- rs(i)=forces(3*(i-1)+j)
+ rs(i-iposc+1)=forces(3*(i-1)+j)
enddo
+#ifdef DEBUG
+ write (iout,*) "j",j," chain",ichain
+ write (iout,*) "rs"
+ write (iout,'(f10.5)') (rs(i),i=1,n)
+#endif
call FDISYS (n,DM(iposc),DU1(iposc),DU2(iposc),rs,xsolv)
+#ifdef DEBUG
+ write (iout,*) "xsolv"
+ write (iout,'(f10.5)') (xsolv(i),i=1,n)
+#endif
ind=1
do i=innt,inct
if (itype(i).eq.10)then
enddo
#ifdef DEBUG
write (iout,*) "d_a_vec"
- write (iout,'(3f10.5)') (d_a_vec(j),j=1,dimen3)
+ write (iout,'(3f10.5)') (d_a_vec(j),j=1,3*(nct-nnt+nside))
#endif
return
end
--- /dev/null
+ subroutine make_SCSC_inter_list
+ implicit none
+ include "DIMENSIONS"
+#ifdef MPI
+ include 'mpif.h'
+ include "COMMON.SETUP"
+#endif
+ include "COMMON.CHAIN"
+ include "COMMON.INTERACT"
+ include "COMMON.SPLITELE"
+ include "COMMON.IOUNITS"
+ double precision xi,yi,zi,xj,yj,zj,xj_safe,yj_safe,zj_safe,
+ & xj_temp,yj_temp,zj_temp
+ double precision dist_init, dist_temp,r_buff_list
+ integer contlisti(200*maxres),contlistj(200*maxres)
+! integer :: newcontlisti(200*nres),newcontlistj(200*nres)
+ integer i,j,itypi,itypj,subchap,xshift,yshift,zshift,iint,
+ & ilist_sc,g_ilist_sc
+ integer displ(0:max_fg_procs),i_ilist_sc(0:max_fg_procs),ierr
+! print *,"START make_SC"
+#ifdef DEBUG
+ write (iout,*) "make_SCSC_inter_list"
+#endif
+ r_buff_list=5.0d0
+ ilist_sc=0
+ do i=iatsc_s,iatsc_e
+ itypi=iabs(itype(i))
+ if (itypi.eq.ntyp1) cycle
+ xi=c(1,nres+i)
+ yi=c(2,nres+i)
+ zi=c(3,nres+i)
+ xi=dmod(xi,boxxsize)
+ if (xi.lt.0) xi=xi+boxxsize
+ yi=dmod(yi,boxysize)
+ if (yi.lt.0) yi=yi+boxysize
+ zi=dmod(zi,boxzsize)
+ if (zi.lt.0) zi=zi+boxzsize
+ do iint=1,nint_gr(i)
+ do j=istart(i,iint),iend(i,iint)
+ itypj=iabs(itype(j))
+ if (itypj.eq.ntyp1) cycle
+ xj=c(1,nres+j)
+ yj=c(2,nres+j)
+ zj=c(3,nres+j)
+ xj=dmod(xj,boxxsize)
+ if (xj.lt.0) xj=xj+boxxsize
+ yj=dmod(yj,boxysize)
+ if (yj.lt.0) yj=yj+boxysize
+ zj=dmod(zj,boxzsize)
+ if (zj.lt.0) zj=zj+boxzsize
+ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ subchap=0
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
+ xj=xj_safe+xshift*boxxsize
+ yj=yj_safe+yshift*boxysize
+ zj=zj_safe+zshift*boxzsize
+ dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ if(dist_temp.lt.dist_init) then
+ dist_init=dist_temp
+ xj_temp=xj
+ yj_temp=yj
+ zj_temp=zj
+ subchap=1
+ endif
+ enddo
+ enddo
+ enddo
+ if (subchap.eq.1) then
+ xj=xj_temp-xi
+ yj=yj_temp-yi
+ zj=zj_temp-zi
+ else
+ xj=xj_safe-xi
+ yj=yj_safe-yi
+ zj=zj_safe-zi
+ endif
+! r_buff_list is a read value for a buffer
+ if (sqrt(dist_init).le.(r_cut_int+r_buff_list)) then
+! Here the list is created
+ ilist_sc=ilist_sc+1
+! this can be substituted by cantor and anti-cantor
+ contlisti(ilist_sc)=i
+ contlistj(ilist_sc)=j
+
+ endif
+ enddo
+ enddo
+ enddo
+! call MPI_Reduce(ilist_sc,g_ilist_sc,1,&
+! MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
+! call MPI_Gather(newnss,1,MPI_INTEGER,&
+! i_newnss,1,MPI_INTEGER,king,FG_COMM,IERR)
+#ifdef MPI
+#ifdef DEBUG
+ write (iout,*) "before MPIREDUCE",ilist_sc
+ do i=1,ilist_sc
+ write (iout,*) i,contlisti(i),contlistj(i)
+ enddo
+#endif
+ if (nfgtasks.gt.1)then
+
+ call MPI_Reduce(ilist_sc,g_ilist_sc,1,
+ & MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
+c write(iout,*) "before bcast",g_ilist_sc
+ call MPI_Gather(ilist_sc,1,MPI_INTEGER,
+ & i_ilist_sc,1,MPI_INTEGER,king,FG_COMM,IERR)
+ displ(0)=0
+ do i=1,nfgtasks-1,1
+ displ(i)=i_ilist_sc(i-1)+displ(i-1)
+ enddo
+! write(iout,*) "before gather",displ(0),displ(1)
+ call MPI_Gatherv(contlisti,ilist_sc,MPI_INTEGER,
+ & newcontlisti,i_ilist_sc,displ,MPI_INTEGER,
+ & king,FG_COMM,IERR)
+ call MPI_Gatherv(contlistj,ilist_sc,MPI_INTEGER,
+ & newcontlistj,i_ilist_sc,displ,MPI_INTEGER,
+ & king,FG_COMM,IERR)
+ call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM,IERR)
+! write(iout,*) "before bcast",g_ilist_sc
+! call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM)
+ call MPI_Bcast(newcontlisti,g_ilist_sc,MPI_INT,king,FG_COMM,
+ & IERR)
+ call MPI_Bcast(newcontlistj,g_ilist_sc,MPI_INT,king,FG_COMM,
+ & IERR)
+! call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM)
+
+ else
+#endif
+ g_ilist_sc=ilist_sc
+
+ do i=1,ilist_sc
+ newcontlisti(i)=contlisti(i)
+ newcontlistj(i)=contlistj(i)
+ enddo
+#ifdef MPI
+ endif
+#endif
+#ifdef DEBUG
+ write (iout,*) "after GATHERV",g_ilist_sc
+ do i=1,g_ilist_sc
+ write (iout,*) i,newcontlisti(i),newcontlistj(i)
+ enddo
+#endif
+ call int_bounds(g_ilist_sc,g_listscsc_start,g_listscsc_end)
+ return
+ end
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ subroutine make_SCp_inter_list
+ implicit none
+ include "DIMENSIONS"
+#ifdef MPI
+ include 'mpif.h'
+ include "COMMON.SETUP"
+#endif
+ include "COMMON.CHAIN"
+ include "COMMON.INTERACT"
+ include "COMMON.SPLITELE"
+ include "COMMON.IOUNITS"
+ double precision xi,yi,zi,xj,yj,zj,xj_safe,yj_safe,zj_safe,
+ & xj_temp,yj_temp,zj_temp
+ double precision dist_init, dist_temp,r_buff_list
+ integer contlistscpi(200*maxres),contlistscpj(200*maxres)
+! integer :: newcontlistscpi(200*nres),newcontlistscpj(200*nres)
+ integer i,j,itypi,itypj,subchap,xshift,yshift,zshift,iint,
+ & ilist_scp,g_ilist_scp
+ integer displ(0:max_fg_procs),i_ilist_scp(0:max_fg_procs),ierr
+ integer contlistscpi_f(200*maxres),contlistscpj_f(200*maxres)
+ integer ilist_scp_first,ifirstrun,g_ilist_sc
+! print *,"START make_SC"
+#ifdef DEBUG
+ write (iout,*) "make_SCp_inter_list"
+#endif
+ r_buff_list=5.0
+ ilist_scp=0
+ ilist_scp_first=0
+ do i=iatscp_s,iatscp_e
+ if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
+ xi=0.5D0*(c(1,i)+c(1,i+1))
+ yi=0.5D0*(c(2,i)+c(2,i+1))
+ zi=0.5D0*(c(3,i)+c(3,i+1))
+ xi=mod(xi,boxxsize)
+ if (xi.lt.0) xi=xi+boxxsize
+ yi=mod(yi,boxysize)
+ if (yi.lt.0) yi=yi+boxysize
+ zi=mod(zi,boxzsize)
+ if (zi.lt.0) zi=zi+boxzsize
+
+ do iint=1,nscp_gr(i)
+
+ do j=iscpstart(i,iint),iscpend(i,iint)
+ itypj=iabs(itype(j))
+ if (itypj.eq.ntyp1) cycle
+! Uncomment following three lines for SC-p interactions
+! xj=c(1,nres+j)-xi
+! yj=c(2,nres+j)-yi
+! zj=c(3,nres+j)-zi
+! Uncomment following three lines for Ca-p interactions
+! xj=c(1,j)-xi
+! yj=c(2,j)-yi
+! zj=c(3,j)-zi
+ xj=c(1,j)
+ yj=c(2,j)
+ zj=c(3,j)
+ xj=mod(xj,boxxsize)
+ if (xj.lt.0) xj=xj+boxxsize
+ yj=mod(yj,boxysize)
+ if (yj.lt.0) yj=yj+boxysize
+ zj=mod(zj,boxzsize)
+ if (zj.lt.0) zj=zj+boxzsize
+ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ subchap=0
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
+ xj=xj_safe+xshift*boxxsize
+ yj=yj_safe+yshift*boxysize
+ zj=zj_safe+zshift*boxzsize
+ dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ if(dist_temp.lt.dist_init) then
+ dist_init=dist_temp
+ xj_temp=xj
+ yj_temp=yj
+ zj_temp=zj
+ subchap=1
+ endif
+ enddo
+ enddo
+ enddo
+ if (subchap.eq.1) then
+ xj=xj_temp-xi
+ yj=yj_temp-yi
+ zj=zj_temp-zi
+ else
+ xj=xj_safe-xi
+ yj=yj_safe-yi
+ zj=zj_safe-zi
+ endif
+#ifdef DEBUG
+ ! r_buff_list is a read value for a buffer
+ if ((sqrt(dist_init).le.(r_cut_int)).and.(ifirstrun.eq.0))
+ & then
+! Here the list is created
+ ilist_scp_first=ilist_scp_first+1
+! this can be substituted by cantor and anti-cantor
+ contlistscpi_f(ilist_scp_first)=i
+ contlistscpj_f(ilist_scp_first)=j
+ endif
+#endif
+! r_buff_list is a read value for a buffer
+ if (sqrt(dist_init).le.(r_cut_int+r_buff_list)) then
+! Here the list is created
+ ilist_scp=ilist_scp+1
+! this can be substituted by cantor and anti-cantor
+ contlistscpi(ilist_scp)=i
+ contlistscpj(ilist_scp)=j
+ endif
+ enddo
+ enddo
+ enddo
+#ifdef MPI
+#ifdef DEBUG
+ write (iout,*) "before MPIREDUCE",ilist_scp
+ do i=1,ilist_scp
+ write (iout,*) i,contlistscpi(i),contlistscpj(i)
+ enddo
+#endif
+ if (nfgtasks.gt.1)then
+
+ call MPI_Reduce(ilist_scp,g_ilist_scp,1,
+ & MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
+c write(iout,*) "before bcast",g_ilist_sc
+ call MPI_Gather(ilist_scp,1,MPI_INTEGER,
+ & i_ilist_scp,1,MPI_INTEGER,king,FG_COMM,IERR)
+ displ(0)=0
+ do i=1,nfgtasks-1,1
+ displ(i)=i_ilist_scp(i-1)+displ(i-1)
+ enddo
+! write(iout,*) "before gather",displ(0),displ(1)
+ call MPI_Gatherv(contlistscpi,ilist_scp,MPI_INTEGER,
+ & newcontlistscpi,i_ilist_scp,displ,MPI_INTEGER,
+ & king,FG_COMM,IERR)
+ call MPI_Gatherv(contlistscpj,ilist_scp,MPI_INTEGER,
+ & newcontlistscpj,i_ilist_scp,displ,MPI_INTEGER,
+ & king,FG_COMM,IERR)
+ call MPI_Bcast(g_ilist_scp,1,MPI_INT,king,FG_COMM,IERR)
+! write(iout,*) "before bcast",g_ilist_sc
+! call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM)
+ call MPI_Bcast(newcontlistscpi,g_ilist_scp,MPI_INT,king,FG_COMM,
+ & IERR)
+ call MPI_Bcast(newcontlistscpj,g_ilist_scp,MPI_INT,king,FG_COMM,
+ & IERR)
+! call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM)
+ else
+#endif
+ g_ilist_scp=ilist_scp
+
+ do i=1,ilist_scp
+ newcontlistscpi(i)=contlistscpi(i)
+ newcontlistscpj(i)=contlistscpj(i)
+ enddo
+#ifdef MPI
+ endif
+#endif
+#ifdef DEBUG
+ write (iout,*) "after MPIREDUCE",g_ilist_scp
+ do i=1,g_ilist_scp
+ write (iout,*) i,newcontlistscpi(i),newcontlistscpj(i)
+ enddo
+
+! if (ifirstrun.eq.0) ifirstrun=1
+! do i=1,ilist_scp_first
+! do j=1,g_ilist_scp
+! if ((newcontlistscpi(j).eq.contlistscpi_f(i)).and.&
+! (newcontlistscpj(j).eq.contlistscpj_f(i))) go to 126
+! enddo
+! print *,itime_mat,"ERROR matrix needs updating"
+! print *,contlistscpi_f(i),contlistscpj_f(i)
+! 126 continue
+! enddo
+#endif
+ call int_bounds(g_ilist_scp,g_listscp_start,g_listscp_end)
+
+ return
+ end
+!-----------------------------------------------------------------------------
+ subroutine make_pp_inter_list
+ implicit none
+ include "DIMENSIONS"
+#ifdef MPI
+ include 'mpif.h'
+ include "COMMON.SETUP"
+#endif
+ include "COMMON.CHAIN"
+ include "COMMON.INTERACT"
+ include "COMMON.SPLITELE"
+ include "COMMON.IOUNITS"
+ double precision xi,yi,zi,xj,yj,zj,xj_safe,yj_safe,zj_safe,
+ & xj_temp,yj_temp,zj_temp
+ double precision xmedj,ymedj,zmedj
+ double precision dist_init, dist_temp,r_buff_list,dxi,dyi,dzi,
+ & xmedi,ymedi,zmedi
+ double precision dx_normi,dy_normi,dz_normi,dxj,dyj,dzj,
+ & dx_normj,dy_normj,dz_normj
+ integer contlistppi(200*maxres),contlistppj(200*maxres)
+! integer :: newcontlistppi(200*nres),newcontlistppj(200*nres)
+ integer i,j,itypi,itypj,subchap,xshift,yshift,zshift,iint,
+ & ilist_pp,g_ilist_pp
+ integer displ(0:max_fg_procs),i_ilist_pp(0:max_fg_procs),ierr
+! print *,"START make_SC"
+#ifdef DEBUG
+ write (iout,*) "make_pp_inter_list"
+#endif
+ ilist_pp=0
+ r_buff_list=5.0
+ do i=iatel_s,iatel_e
+ if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
+ dxi=dc(1,i)
+ dyi=dc(2,i)
+ dzi=dc(3,i)
+ dx_normi=dc_norm(1,i)
+ dy_normi=dc_norm(2,i)
+ dz_normi=dc_norm(3,i)
+ xmedi=c(1,i)+0.5d0*dxi
+ ymedi=c(2,i)+0.5d0*dyi
+ zmedi=c(3,i)+0.5d0*dzi
+ xmedi=dmod(xmedi,boxxsize)
+ if (xmedi.lt.0) xmedi=xmedi+boxxsize
+ ymedi=dmod(ymedi,boxysize)
+ if (ymedi.lt.0) ymedi=ymedi+boxysize
+ zmedi=dmod(zmedi,boxzsize)
+ if (zmedi.lt.0) zmedi=zmedi+boxzsize
+ do j=ielstart(i),ielend(i)
+! write (iout,*) i,j,itype(i),itype(j)
+ if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1) cycle
+
+! 1,j)
+ dxj=dc(1,j)
+ dyj=dc(2,j)
+ dzj=dc(3,j)
+ dx_normj=dc_norm(1,j)
+ dy_normj=dc_norm(2,j)
+ dz_normj=dc_norm(3,j)
+! xj=c(1,j)+0.5D0*dxj-xmedi
+! yj=c(2,j)+0.5D0*dyj-ymedi
+! zj=c(3,j)+0.5D0*dzj-zmedi
+ xj=c(1,j)+0.5D0*dxj
+ yj=c(2,j)+0.5D0*dyj
+ zj=c(3,j)+0.5D0*dzj
+ xj=mod(xj,boxxsize)
+ if (xj.lt.0) xj=xj+boxxsize
+ yj=mod(yj,boxysize)
+ if (yj.lt.0) yj=yj+boxysize
+ zj=mod(zj,boxzsize)
+ if (zj.lt.0) zj=zj+boxzsize
+
+ dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
+ xj=xj_safe+xshift*boxxsize
+ yj=yj_safe+yshift*boxysize
+ zj=zj_safe+zshift*boxzsize
+ dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
+ if(dist_temp.lt.dist_init) then
+ dist_init=dist_temp
+ xj_temp=xj
+ yj_temp=yj
+ zj_temp=zj
+ endif
+ enddo
+ enddo
+ enddo
+
+ if (sqrt(dist_init).le.(r_cut_int+r_buff_list)) then
+! Here the list is created
+ ilist_pp=ilist_pp+1
+! this can be substituted by cantor and anti-cantor
+ contlistppi(ilist_pp)=i
+ contlistppj(ilist_pp)=j
+ endif
+ enddo
+ enddo
+! enddo
+#ifdef MPI
+#ifdef DEBUG
+ write (iout,*) "before MPIREDUCE",ilist_pp
+ do i=1,ilist_pp
+ write (iout,*) i,contlistppi(i),contlistppj(i)
+ enddo
+#endif
+ if (nfgtasks.gt.1)then
+
+ call MPI_Reduce(ilist_pp,g_ilist_pp,1,
+ & MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
+! write(iout,*) "before bcast",g_ilist_sc
+ call MPI_Gather(ilist_pp,1,MPI_INTEGER,
+ & i_ilist_pp,1,MPI_INTEGER,king,FG_COMM,IERR)
+ displ(0)=0
+ do i=1,nfgtasks-1,1
+ displ(i)=i_ilist_pp(i-1)+displ(i-1)
+ enddo
+! write(iout,*) "before gather",displ(0),displ(1)
+ call MPI_Gatherv(contlistppi,ilist_pp,MPI_INTEGER,
+ & newcontlistppi,i_ilist_pp,displ,MPI_INTEGER,
+ & king,FG_COMM,IERR)
+ call MPI_Gatherv(contlistppj,ilist_pp,MPI_INTEGER,
+ & newcontlistppj,i_ilist_pp,displ,MPI_INTEGER,
+ & king,FG_COMM,IERR)
+ call MPI_Bcast(g_ilist_pp,1,MPI_INT,king,FG_COMM,IERR)
+! write(iout,*) "before bcast",g_ilist_sc
+! call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM)
+ call MPI_Bcast(newcontlistppi,g_ilist_pp,MPI_INT,king,FG_COMM,
+ & IERR)
+ call MPI_Bcast(newcontlistppj,g_ilist_pp,MPI_INT,king,FG_COMM,
+ & IERR)
+
+! call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM)
+
+ else
+#endif
+ g_ilist_pp=ilist_pp
+
+ do i=1,ilist_pp
+ newcontlistppi(i)=contlistppi(i)
+ newcontlistppj(i)=contlistppj(i)
+ enddo
+#ifdef MPI
+ endif
+#endif
+ call int_bounds(g_ilist_pp,g_listpp_start,g_listpp_end)
+#ifdef DEBUG
+ write (iout,*) "after MPIREDUCE",g_ilist_pp
+ do i=1,g_ilist_pp
+ write (iout,*) i,newcontlistppi(i),newcontlistppj(i)
+ enddo
+#endif
+ return
+ end
timlim=60.0D0*timlim
safety = 60.0d0*safety
modecalc=0
+ call readi(controlcard,"INTER_LIST_UPDATE",imatupdate,100)
call reada(controlcard,"T_BATH",t_bath,300.0d0)
minim=(index(controlcard,'MINIMIZE').gt.0)
dccart=(index(controlcard,'CART').gt.0)
& "Initial time step of numerical integration:",d_time,
& " natural units"
write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
+ write (iout,'(a60,f10.5,a)') "Cutoff on interactions",r_cut_int,
+ & " A"
+ write(iout,'(a60,i5)')"Frequency of updating interaction list",
+ & imatupdate
if (RESPA) then
write (iout,'(2a,i4,a)')
& "A-MTS algorithm used; initial time step for fast-varying",
c & " n",n," iposc",iposc,iposc+n-1
innt=chain_border(1,ichain)
inct=chain_border(2,ichain)
+c diagnostics
+c innt=chain_border(1,1)
+c inct=chain_border(2,1)
do i=innt,inct
vvec(ind+1)=v_work(j,i)
ind=ind+1
innt=chain_border(1,ichain)
inct=chain_border(2,ichain)
iposc=iposd_chain(ichain)
+c for debugging only
+c innt=chain_border(1,1)
+c inct=chain_border(2,1)
+c iposc=iposd_chain(1)
c write (iout,*)"stochastic_force ichain=",ichain," innt",innt,
c & " inct",inct," iposc",iposc
do j=1,3
projects / unres.git / commitdiff