8 use control_data, only:maxterm_sccor
10 !-----------------------------------------------------------------------------
11 ! Max. number of residue types and parameters in expressions for
12 ! virtual-bond angle bending potentials
13 ! integer,parameter :: maxthetyp=3
14 ! integer,parameter :: maxthetyp1=maxthetyp+1
16 ! & maxtheterm2=6,maxtheterm3=4,maxsingle=6,maxdouble=4,
17 ! & mmaxtheterm=maxtheterm)
18 !-----------------------------------------------------------------------------
19 ! Max. number of types of dihedral angles & multiplicity of torsional barriers
20 ! and the number of terms in double torsionals
21 ! integer,parameter :: maxlor=3,maxtermd_1=8,maxtermd_2=8
22 ! parameter (maxtor=4,maxterm=10)
23 !-----------------------------------------------------------------------------
24 ! Max number of torsional terms in SCCOR
25 !el integer,parameter :: maxterm_sccor=6
26 !-----------------------------------------------------------------------------
27 character(len=1),dimension(:),allocatable :: secstruc !(maxres)
28 !-----------------------------------------------------------------------------
31 !-----------------------------------------------------------------------------
33 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
34 !-----------------------------------------------------------------------------
36 !-----------------------------------------------------------------------------
37 subroutine write_rbank(jlee,adif,nft)
40 use geometry_data, only: nres,rad2deg
41 ! implicit real*8 (a-h,o-z)
42 ! include 'DIMENSIONS'
43 ! include 'COMMON.IOUNITS'
44 ! include 'COMMON.CSA'
45 ! include 'COMMON.BANK'
46 ! include 'COMMON.CHAIN'
47 ! include 'COMMON.GEO'
49 integer :: nft,i,k,j,l,jlee
52 open(icsa_rbank,file=csa_rbank,status="unknown")
53 write (icsa_rbank,900) jlee,nbank,nstep,nft,icycle,adif
55 write (icsa_rbank,952) k,rene(k),rrmsn(k),rpncn(k)
58 write (icsa_rbank,850) (rad2deg*rvar(i,l,j,k),i=1,4)
65 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
67 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
71 end subroutine write_rbank
72 !-----------------------------------------------------------------------------
73 subroutine read_rbank(jlee,adif)
76 use geometry_data, only: nres,deg2rad
78 ! implicit real*8 (a-h,o-z)
79 ! include 'DIMENSIONS'
81 ! include 'COMMON.IOUNITS'
82 ! include 'COMMON.CSA'
83 ! include 'COMMON.BANK'
84 ! include 'COMMON.CHAIN'
85 ! include 'COMMON.GEO'
86 ! include 'COMMON.SETUP'
87 character(len=80) :: karta
89 integer :: nbankr,nstepr,nftr,icycler,kk,k,j,l,i,&
90 ierror,ierrcode,jlee,jleer
93 open(icsa_rbank,file=csa_rbank,status="old")
94 read (icsa_rbank,901) jleer,nbankr,nstepr,nftr,icycler,adif
95 print *,jleer,nbankr,nstepr,nftr,icycler,adif
96 ! print *, 'adif from read_rbank ',adif
98 if(nbankr.ne.nbank) then
99 write (iout,*) 'ERROR in READ_BANK: NBANKR',nbankr,&
101 call mpi_abort(mpi_comm_world,ierror,ierrcode)
103 if(jleer.ne.jlee) then
104 write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,&
106 call mpi_abort(mpi_comm_world,ierror,ierrcode)
112 read (icsa_rbank,'(a80)') karta
113 write(iout,*) "READ_RBANK: kk=",kk
115 ! if (index(karta,"*").gt.0) then
116 ! write (iout,*) "***** Stars in bankr ***** k=",k,
120 ! read (30,850) (rdummy,i=1,4)
125 call reada(karta,"total E",rene(kk),1.0d20)
126 call reada(karta,"rmsd from N",rrmsn(kk),0.0d0)
127 call reada(karta,"%NC",rpncn(kk),0.0d0)
128 write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
129 "%NC",bpncn(kk),ibank(kk)
130 ! read (icsa_rbank,953) kdummy,rene(kk),rrmsn(kk),rpncn(kk)
133 read (icsa_rbank,850) (rvar(i,l,j,kk),i=1,4)
134 ! write (iout,850) (rvar(i,l,j,kk),i=1,4)
136 rvar(i,l,j,kk)=deg2rad*rvar(i,l,j,kk)
142 !d write (*,*) "read_rbank ******************* kk",kk,
144 if (kk.lt.nbankr) nbankr=kk
149 !d write (*,850) (rvar(i,l,j,kk),i=1,4)
156 901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
157 953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2)
160 end subroutine read_rbank
161 !-----------------------------------------------------------------------------
162 subroutine write_bank(jlee,nft)
165 use control_data, only: vdisulf
166 use geometry_data, only: nres,rad2deg
167 ! implicit real*8 (a-h,o-z)
168 ! include 'DIMENSIONS'
169 ! include 'COMMON.IOUNITS'
170 ! include 'COMMON.CSA'
171 ! include 'COMMON.BANK'
172 ! include 'COMMON.CHAIN'
173 ! include 'COMMON.GEO'
174 ! include 'COMMON.SBRIDGE'
175 ! include 'COMMON.CONTROL'
176 character(len=7) :: chtmp
177 character(len=40) :: chfrm
180 integer :: nft,k,l,i,j,jlee
182 open(icsa_bank,file=csa_bank,status="unknown")
183 write (icsa_bank,900) jlee,nbank,nstep,nft,icycle,cutdif
184 write (icsa_bank,902) nglob_csa, eglob_csa
185 open (igeom,file=intname,status='UNKNOWN')
187 write (icsa_bank,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
188 if (vdisulf) write (icsa_bank,'(101i4)') &
189 bvar_nss(k),((bvar_ss(j,i,k),j=1,2),i=1,bvar_nss(k))
192 write (icsa_bank,850) (rad2deg*bvar(i,l,j,k),i=1,4)
195 if (bvar_nss(k).le.9) then
196 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
197 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
199 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
200 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
201 write (igeom,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
202 bvar_ss(2,i,k),i=10,bvar_nss(k))
204 write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
205 write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
206 write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
207 write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
212 if (nstep/200.gt.ilastnstep) then
214 ilastnstep=(ilastnstep+1)*1.5
215 write(chfrm,'(a2,i1,a1)') '(i',int(dlog10(dble(nstep))+1),')'
216 write(chtmp,chfrm) nstep
217 open(icsa_int,file=prefix(:ilen(prefix)) &
218 //'_'//chtmp(:ilen(chtmp))//'.int',status='UNKNOWN')
220 if (bvar_nss(k).le.9) then
221 write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
222 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
224 write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
225 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
226 write (icsa_int,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
227 bvar_ss(2,i,k),i=10,bvar_nss(k))
229 write (icsa_int,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
230 write (icsa_int,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
231 write (icsa_int,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
232 write (icsa_int,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
240 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
242 902 format (1x,'nglob_csa =',i4,' eglob_csa =',1pe14.5)
243 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
247 end subroutine write_bank
248 !-----------------------------------------------------------------------------
249 subroutine write_bank_reminimized(jlee,nft)
252 use geometry_data, only: nres,rad2deg
254 ! implicit real*8 (a-h,o-z)
255 ! include 'DIMENSIONS'
256 ! include 'COMMON.IOUNITS'
257 ! include 'COMMON.CSA'
258 ! include 'COMMON.BANK'
259 ! include 'COMMON.CHAIN'
260 ! include 'COMMON.GEO'
261 ! include 'COMMON.SBRIDGE'
263 integer :: nft,i,l,j,k,jlee
265 open(icsa_bank_reminimized,file=csa_bank_reminimized,&
267 write (icsa_bank_reminimized,900) &
268 jlee,nbank,nstep,nft,icycle,cutdif
269 open (igeom,file=intname,status='UNKNOWN')
271 write (icsa_bank_reminimized,952) k,bene(k),brmsn(k),&
275 write (icsa_bank_reminimized,850) (rad2deg*bvar(i,l,j,k),i=1,4)
279 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
280 nss,(ihpb(i),jhpb(i),i=1,nss)
282 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
283 nss,(ihpb(i),jhpb(i),i=1,9)
284 write (igeom,'(3X,11(1X,2I3))') (ihpb(i),jhpb(i),i=10,nss)
286 write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
287 write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
288 write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
289 write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
291 close(icsa_bank_reminimized)
296 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
298 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
302 end subroutine write_bank_reminimized
303 !-----------------------------------------------------------------------------
304 subroutine read_bank(jlee,nft,cutdifr)
307 use control_data, only: vdisulf
308 use geometry_data, only: nres,deg2rad
310 ! implicit real*8 (a-h,o-z)
311 ! include 'DIMENSIONS'
312 ! include 'COMMON.IOUNITS'
313 ! include 'COMMON.CSA'
314 ! include 'COMMON.BANK'
315 ! include 'COMMON.CHAIN'
316 ! include 'COMMON.GEO'
317 ! include 'COMMON.CONTROL'
318 ! include 'COMMON.SBRIDGE'
319 character(len=80) :: karta
323 integer :: nft,kk,k,l,i,j,jlee
324 real(kind=8) :: cutdifr
326 open(icsa_bank,file=csa_bank,status="old")
327 read (icsa_bank,901) jlee,nbank,nstep,nft,icycle,cutdifr
328 read (icsa_bank,902) nglob_csa, eglob_csa
329 ! if(jleer.ne.jlee) then
330 ! write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,
332 ! call mpi_abort(mpi_comm_world,ierror,ierrcode)
337 read (icsa_bank,'(a80)') karta
338 write(iout,*) "READ_BANK: kk=",kk
340 ! if (index(karta,"*").gt.0) then
341 ! write (iout,*) "***** Stars in bank ***** k=",k,
345 ! read (33,850) (rdummy,i=1,4)
350 call reada(karta,"total E",bene(kk),1.0d20)
351 call reada(karta,"rmsd from N",brmsn(kk),0.0d0)
352 call reada(karta,"%NC",bpncn(kk),0.0d0)
353 read (karta(ilen(karta)-1:),*,end=111,err=111) ibank(kk)
357 write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
358 "%NC",bpncn(kk),ibank(kk)
359 ! read (icsa_bank,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
361 read (icsa_bank,'(101i4)') &
362 bvar_nss(kk),((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
363 bvar_ns(kk)=ns-2*bvar_nss(kk)
364 write(iout,*) 'read SSBOND',bvar_nss(kk),&
365 ((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
366 !d write(iout,*) 'read CYS #free ', bvar_ns(kk)
370 do while( iss(i).ne.bvar_ss(1,j,kk)-nres .and. &
371 iss(i).ne.bvar_ss(2,j,kk)-nres .and. &
375 if (j.gt.bvar_nss(kk)) then
380 !d write(iout,*)'read CYS free',(bvar_s(l,kk),l=1,bvar_ns(kk))
384 read (icsa_bank,850) (bvar(i,l,j,kk),i=1,4)
385 ! write (iout,850) (bvar(i,l,j,kk),i=1,4)
387 bvar(i,l,j,kk)=deg2rad*bvar(i,l,j,kk)
394 if (kk.lt.nbank) nbank=kk
395 !d write (*,*) "read_bank ******************* kk",kk,
401 !d write (*,850) (bvar(i,l,j,kk),i=1,4)
407 ! read (33,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
410 ! read (33,850) (bvar(i,l,j,k),i=1,4)
412 ! bvar(i,l,j,k)=deg2rad*bvar(i,l,j,k)
420 952 format (1x,'#',i4,' total E ',f12.3,' rmsd from N ',f8.3,i5)
421 901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
422 902 format (1x,11x,i4,12x,1pe14.5)
423 953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2,i5)
426 end subroutine read_bank
427 !-----------------------------------------------------------------------------
428 subroutine write_bank1(jlee)
431 use geometry_data, only: nres,rad2deg
432 ! implicit real*8 (a-h,o-z)
433 ! include 'DIMENSIONS'
434 ! include 'COMMON.IOUNITS'
435 ! include 'COMMON.CSA'
436 ! include 'COMMON.BANK'
437 ! include 'COMMON.CHAIN'
438 ! include 'COMMON.GEO'
440 integer :: k,i,l,j,jlee
442 #if defined(AIX) || defined(PGI)
443 open(icsa_bank1,file=csa_bank1,position="append")
445 open(icsa_bank1,file=csa_bank1,access="append")
447 write (icsa_bank1,900) jlee,nbank,nstep,cutdif
449 write (icsa_bank1,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
452 write (icsa_bank1,850) (rad2deg*bvar(i,l,j,k),i=1,4)
458 900 format (4x,"jlee =",i5,3x,"nbank =",i5,3x,"nstep =",i10,f15.5)
459 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
463 end subroutine write_bank1
464 !-----------------------------------------------------------------------------
466 !-----------------------------------------------------------------------------
467 ! subroutine cartprint
469 ! use geometry_data, only: c
470 ! use energy_data, only: itype
471 ! implicit real*8 (a-h,o-z)
472 ! include 'DIMENSIONS'
473 ! include 'COMMON.CHAIN'
474 ! include 'COMMON.INTERACT'
475 ! include 'COMMON.NAMES'
476 ! include 'COMMON.IOUNITS'
481 ! write (iout,110) restyp(itype(i,1)),i,c(1,i),c(2,i),&
482 ! c(3,i),c(1,nres+i),c(2,nres+i),c(3,nres+i)
484 ! 100 format (//' alpha-carbon coordinates ',&
485 ! ' centroid coordinates'/ &
486 ! ' ', 6X,'X',11X,'Y',11X,'Z',&
487 ! 10X,'X',11X,'Y',11X,'Z')
488 ! 110 format (a,'(',i3,')',6f12.5)
490 ! end subroutine cartprint
491 !-----------------------------------------------------------------------------
493 !-----------------------------------------------------------------------------
494 subroutine secstrp2dihc
498 ! implicit real*8 (a-h,o-z)
499 ! include 'DIMENSIONS'
500 ! include 'COMMON.GEO'
501 ! include 'COMMON.BOUNDS'
502 ! include 'COMMON.CHAIN'
503 ! include 'COMMON.TORCNSTR'
504 ! include 'COMMON.IOUNITS'
505 !el character(len=1),dimension(nres) :: secstruc !(maxres)
506 !el COMMON/SECONDARYS/secstruc
507 character(len=80) :: line
512 integer :: i,ii,lenpre
514 allocate(secstruc(nres))
516 !dr call getenv_loc('SECPREDFIL',secpred)
518 secpred=prefix(:lenpre)//'.spred'
520 #if defined(WINIFL) || defined(WINPGI)
521 open(isecpred,file=secpred,status='old',readonly,shared)
522 #elif (defined CRAY) || (defined AIX)
523 open(isecpred,file=secpred,status='old',action='read')
525 open(isecpred,file=secpred,status='old')
527 open(isecpred,file=secpred,status='old',action='read')
529 ! read secondary structure prediction from JPRED here!
530 ! read(isecpred,'(A80)',err=100,end=100) line
531 ! read(line,'(f10.3)',err=110) ftors
532 read(isecpred,'(f10.3)',err=110) ftors(1)
534 write (iout,*) 'FTORS factor =',ftors(1)
535 ! initialize secstruc to any
542 call read_secstr_pred(isecpred,iout,errflag)
544 write(iout,*)'There is a problem with the list of secondary-',&
545 'structure prediction'
548 ! 8/13/98 Set limits to generating the dihedral angles
557 if ( secstruc(i) .eq. 'H') then
558 ! Helix restraints for this residue
561 phi0(ii) = 45.0D0*deg2rad
562 drange(ii)= 5.0D0*deg2rad
563 phibound(1,i) = phi0(ii)-drange(ii)
564 phibound(2,i) = phi0(ii)+drange(ii)
565 else if (secstruc(i) .eq. 'E') then
566 ! strand restraints for this residue
569 phi0(ii) = 180.0D0*deg2rad
570 drange(ii)= 5.0D0*deg2rad
571 phibound(1,i) = phi0(ii)-drange(ii)
572 phibound(2,i) = phi0(ii)+drange(ii)
574 ! no restraints for this residue
575 ndih_nconstr=ndih_nconstr+1
576 idih_nconstr(ndih_nconstr)=i
580 ! deallocate(secstruc)
583 write(iout,'(A30,A80)')'Error reading file SECPRED',secpred
584 ! deallocate(secstruc)
587 write(iout,'(A20)')'Error reading FTORS'
588 ! deallocate(secstruc)
590 end subroutine secstrp2dihc
591 !-----------------------------------------------------------------------------
592 subroutine read_secstr_pred(jin,jout,errors)
594 ! implicit real*8 (a-h,o-z)
595 ! INCLUDE 'DIMENSIONS'
596 ! include 'COMMON.IOUNITS'
597 ! include 'COMMON.CHAIN'
598 !el character(len=1),dimension(nres) :: secstruc !(maxres)
599 !el COMMON/SECONDARYS/secstruc
601 character(len=80) :: line,line1 !,ucase
602 logical :: errflag,errors,blankline
605 integer :: jin,jout,iseq,ipos,ipos1,iend,il,&
608 read (jin,'(a)') line
609 write (jout,'(2a)') '> ',line(1:78)
611 ! Remember that we number full residues starting from 2, then, iseq=1 and iseq=nres
612 ! correspond to the end-groups. ADD to the secondary structure prediction "-" for the
613 ! end-groups in the input file "*.spred"
616 do while (index(line1,'$END').eq.0)
617 ! Override commented lines.
620 do while (.not.blankline)
622 call mykey(line,line1,ipos,blankline,errflag)
623 if (errflag) write (jout,'(2a)') &
624 'Error when reading sequence in line: ',line
625 errors=errors .or. errflag
626 if (.not. blankline .and. .not. errflag) then
629 !el if (iseq.le.maxres) then
630 if (line1(1:1).eq.'-' ) then
631 secstruc(iseq)=line1(1:1)
632 else if ( ( ucase(line1(1:1)).eq.'E' ) .or. &
633 ( ucase(line1(1:1)).eq.'H' ) ) then
634 secstruc(iseq)=ucase(line1(1:1))
637 write (jout,1010) line1(1:1), iseq
642 !el write (jout,1000) iseq,maxres
645 do while (ipos1.le.iend)
650 !el if (iseq.le.maxres) then
651 if (line1(ipos1-1:ipos1-1).eq.'-' ) then
652 secstruc(iseq)=line1(ipos1-1:ipos1-1)
653 else if((ucase(line1(ipos1-1:ipos1-1)).eq.'E').or. &
654 (ucase(line1(ipos1-1:ipos1-1)).eq.'H') ) then
655 secstruc(iseq)=ucase(line1(ipos1-1:ipos1-1))
658 write (jout,1010) line1(ipos1-1:ipos1-1), iseq
663 !el write (jout,1000) iseq,maxres
670 read (jin,'(a)') line
671 write (jout,'(2a)') '> ',line(1:78)
675 !d write (jout,'(10a8)') (sequence(i),i=1,iseq-1)
677 !d check whether the found length of the chain is correct.
678 length_of_chain=iseq-1
679 if (length_of_chain .ne. nres) then
681 write (jout,'(a,i4,a,i4,a)') &
682 'Error: the number of labels specified in $SEC_STRUC_PRED (' &
683 ,length_of_chain,') does not match with the number of residues (' &
688 1000 format('Error - the number of residues (',i4,&
689 ') has exceeded maximum (',i4,').')
690 1010 format ('Error - unrecognized secondary structure label',a4,&
693 end subroutine read_secstr_pred
695 !-----------------------------------------------------------------------------
697 !-----------------------------------------------------------------------------
702 use control_data, only:maxterm !,maxtor
706 use control, only: getenv_loc
708 ! Read the parameters of the probability distributions of the virtual-bond
709 ! valence angles and the side chains and energy parameters.
711 ! Important! Energy-term weights ARE NOT read here; they are read from the
712 ! main input file instead, because NO defaults have yet been set for these
715 ! implicit real*8 (a-h,o-z)
716 ! include 'DIMENSIONS'
721 ! include 'COMMON.IOUNITS'
722 ! include 'COMMON.CHAIN'
723 ! include 'COMMON.INTERACT'
724 ! include 'COMMON.GEO'
725 ! include 'COMMON.LOCAL'
726 ! include 'COMMON.TORSION'
727 ! include 'COMMON.SCCOR'
728 ! include 'COMMON.SCROT'
729 ! include 'COMMON.FFIELD'
730 ! include 'COMMON.NAMES'
731 ! include 'COMMON.SBRIDGE'
732 ! include 'COMMON.MD'
733 ! include 'COMMON.SETUP'
734 character(len=1) :: t1,t2,t3
735 character(len=1) :: onelett(4) = (/"G","A","P","D"/)
736 character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/)
737 logical :: lprint,LaTeX,SPLIT_FOURIERTOR
738 real(kind=8),dimension(3,3,maxlob) :: blower !(3,3,maxlob)
739 real(kind=8),dimension(13) :: buse
740 character(len=3) :: lancuch !,ucase
742 integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm,jj
743 integer :: maxinter,junk,kk,ii,ncatprotparm,nkcctyp
744 real(kind=8) :: v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,si,&
745 dwa16,rjunk,akl,v0ij,rri,epsij,rrij,sigeps,sigt1sq,&
746 sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,&
747 res1,epsijlip,epspeptube,epssctube,sigmapeptube, &
749 integer :: ichir1,ichir2,ijunk
752 ! real(kind=8),dimension(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) :: v1_el,v2_el !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
753 !el allocate(v1_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
754 !el allocate(v2_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
756 ! For printing parameters after they are read set the following in the UNRES
759 ! setenv PRINT_PARM YES
761 ! To print parameters in LaTeX format rather than as ASCII tables:
765 call getenv_loc("PRINT_PARM",lancuch)
766 lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
767 call getenv_loc("LATEX",lancuch)
768 LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
770 dwa16=2.0d0**(1.0d0/6.0d0)
772 ! Assign virtual-bond length
775 vblinv2=vblinv*vblinv
777 ! Read the virtual-bond parameters, masses, and moments of inertia
778 ! and Stokes' radii of the peptide group and side chains
780 allocate(dsc(ntyp1)) !(ntyp1)
781 allocate(dsc_inv(ntyp1)) !(ntyp1)
782 allocate(nbondterm_nucl(ntyp_molec(2))) !(ntyp)
783 allocate(vbldsc0_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
784 allocate(aksc_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
785 allocate(nbondterm(ntyp)) !(ntyp)
786 allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
787 allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
788 allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
789 allocate(long_r_sidechain(ntyp))
790 allocate(short_r_sidechain(ntyp))
795 allocate(msc(ntyp+1)) !(ntyp+1)
796 allocate(isc(ntyp+1)) !(ntyp+1)
797 allocate(restok(ntyp+1)) !(ntyp+1)
799 read (ibond,*) vbldp0,akp,mp,ip,pstok
802 read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
803 dsc(i) = vbldsc0(1,i)
807 dsc_inv(i)=1.0D0/dsc(i)
816 allocate(msc(ntyp+1,5)) !(ntyp+1)
817 allocate(isc(ntyp+1,5)) !(ntyp+1)
818 allocate(restok(ntyp+1,5)) !(ntyp+1)
820 read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp(1),ip(1),pstok(1)
822 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
823 j=1,nbondterm(i)),msc(i,1),isc(i,1),restok(i,1)
824 dsc(i) = vbldsc0(1,i)
828 dsc_inv(i)=1.0D0/dsc(i)
832 read (ibond_nucl,*) vbldp0_nucl,akp_nucl,mp(2),ip(2),pstok(2)
835 read (ibond_nucl,*) vbldsc0_nucl(1,i),aksc_nucl(1,i),msc(i,2),isc(i,2),restok(i,2)
836 ! dsc(i) = vbldsc0_nucl(1,i)
840 ! dsc_inv(i)=1.0D0/dsc(i)
843 ! read (ibond_nucl,*) junk,vbldp0_nucl,akp_nucl,rjunk,mp(2),ip(2),pstok(2)
844 ! do i=1,ntyp_molec(2)
845 ! read (ibond_nucl,*) nbondterm_nucl(i),(vbldsc0_nucl(j,i),&
846 ! aksc_nucl(j,i),abond0_nucl(j,i),&
847 ! j=1,nbondterm_nucl(i)),msc(i,2),isc(i,2),restok(i,2)
848 ! dsc(i) = vbldsc0(1,i)
852 ! dsc_inv(i)=1.0D0/dsc(i)
857 write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
858 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',&
860 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp(1),ip(1),pstok(1)
862 write (iout,'(a10,i3,6f10.5)') restyp(i,1),nbondterm(i),&
863 vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i,1),isc(i,1),restok(i,1)
865 write (iout,'(13x,3f10.5)') &
866 vbldsc0(j,i),aksc(j,i),abond0(j,i)
870 if (oldion.eq.1) then
872 read(iion,*) msc(i,5),restok(i,5)
873 print *,msc(i,5),restok(i,5)
877 read (iion,*) ncatprotparm
878 allocate(catprm(ncatprotparm,4))
880 read (iion,*) (catprm(i,k),i=1,ncatprotparm)
884 ! read (iion,*) (vcatprm(k),k=1,ncatprotpram)
885 !----------------------------------------------------
886 allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
887 allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp)) !(-ntyp:ntyp)
888 allocate(athet(2,-ntyp:ntyp,-1:1,-1:1))
889 allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
890 allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
891 allocate(gthet(3,-ntyp:ntyp)) !(3,-ntyp:ntyp)
901 allocate(liptranene(ntyp))
902 !C reading lipid parameters
903 write (iout,*) "iliptranpar",iliptranpar
905 read(iliptranpar,*) pepliptran
908 read(iliptranpar,*) liptranene(i)
909 print *,liptranene(i)
915 ! Read the parameters of the probability distribution/energy expression
916 ! of the virtual-bond valence angles theta
919 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),&
920 (bthet(j,i,1,1),j=1,2)
921 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
922 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
923 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
927 athet(1,i,1,-1)=athet(1,i,1,1)
928 athet(2,i,1,-1)=athet(2,i,1,1)
929 bthet(1,i,1,-1)=-bthet(1,i,1,1)
930 bthet(2,i,1,-1)=-bthet(2,i,1,1)
931 athet(1,i,-1,1)=-athet(1,i,1,1)
932 athet(2,i,-1,1)=-athet(2,i,1,1)
933 bthet(1,i,-1,1)=bthet(1,i,1,1)
934 bthet(2,i,-1,1)=bthet(2,i,1,1)
938 athet(1,i,-1,-1)=athet(1,-i,1,1)
939 athet(2,i,-1,-1)=-athet(2,-i,1,1)
940 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
941 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
942 athet(1,i,-1,1)=athet(1,-i,1,1)
943 athet(2,i,-1,1)=-athet(2,-i,1,1)
944 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
945 bthet(2,i,-1,1)=bthet(2,-i,1,1)
946 athet(1,i,1,-1)=-athet(1,-i,1,1)
947 athet(2,i,1,-1)=athet(2,-i,1,1)
948 bthet(1,i,1,-1)=bthet(1,-i,1,1)
949 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
954 polthet(j,i)=polthet(j,-i)
957 gthet(j,i)=gthet(j,-i)
965 'Parameters of the virtual-bond valence angles:'
966 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',&
967 ' ATHETA0 ',' A1 ',' A2 ',&
970 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
971 a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
973 write (iout,'(/a/9x,5a/79(1h-))') &
974 'Parameters of the expression for sigma(theta_c):',&
975 ' ALPH0 ',' ALPH1 ',' ALPH2 ',&
976 ' ALPH3 ',' SIGMA0C '
978 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
979 (polthet(j,i),j=0,3),sigc0(i)
981 write (iout,'(/a/9x,5a/79(1h-))') &
982 'Parameters of the second gaussian:',&
983 ' THETA0 ',' SIGMA0 ',' G1 ',&
986 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,theta0(i),&
987 sig0(i),(gthet(j,i),j=1,3)
991 'Parameters of the virtual-bond valence angles:'
992 write (iout,'(/a/9x,5a/79(1h-))') &
993 'Coefficients of expansion',&
994 ' theta0 ',' a1*10^2 ',' a2*10^2 ',&
995 ' b1*10^1 ',' b2*10^1 '
997 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),&
998 a0thet(i),(100*athet(j,i,1,1),j=1,2),&
999 (10*bthet(j,i,1,1),j=1,2)
1001 write (iout,'(/a/9x,5a/79(1h-))') &
1002 'Parameters of the expression for sigma(theta_c):',&
1003 ' alpha0 ',' alph1 ',' alph2 ',&
1004 ' alhp3 ',' sigma0c '
1006 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i,1),&
1007 (polthet(j,i),j=0,3),sigc0(i)
1009 write (iout,'(/a/9x,5a/79(1h-))') &
1010 'Parameters of the second gaussian:',&
1011 ' theta0 ',' sigma0*10^2 ',' G1*10^-1',&
1014 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),theta0(i),&
1015 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
1021 ! Read the parameters of Utheta determined from ab initio surfaces
1022 ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
1024 IF (tor_mode.eq.0) THEN
1025 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,&
1026 ntheterm3,nsingle,ndouble
1027 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
1029 !----------------------------------------------------
1030 allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1031 allocate(aa0thet(-nthetyp-1:nthetyp+1,&
1032 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1033 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1034 allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,&
1035 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1036 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1037 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1038 allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1039 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1040 allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1041 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1042 allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1043 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1044 allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1045 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1046 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1047 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1048 allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1049 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1050 allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1051 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1052 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1053 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1055 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
1057 ithetyp(i)=-ithetyp(-i)
1060 aa0thet(:,:,:,:)=0.0d0
1061 aathet(:,:,:,:,:)=0.0d0
1062 bbthet(:,:,:,:,:,:)=0.0d0
1063 ccthet(:,:,:,:,:,:)=0.0d0
1064 ddthet(:,:,:,:,:,:)=0.0d0
1065 eethet(:,:,:,:,:,:)=0.0d0
1066 ffthet(:,:,:,:,:,:,:)=0.0d0
1067 ggthet(:,:,:,:,:,:,:)=0.0d0
1069 ! VAR:iblock means terminally blocking group 1=non-proline 2=proline
1071 ! VAR:ntethtyp is type of theta potentials type currently 0=glycine
1072 ! VAR:1=non-glicyne non-proline 2=proline
1073 ! VAR:negative values for D-aminoacid
1075 do j=-nthetyp,nthetyp
1076 do k=-nthetyp,nthetyp
1077 read (ithep,'(6a)',end=111,err=111) res1
1078 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
1079 ! VAR: aa0thet is variable describing the average value of Foureir
1080 ! VAR: expansion series
1081 ! VAR: aathet is foureir expansion in theta/2 angle for full formula
1082 ! VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
1083 !ondens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
1084 read (ithep,*,end=111,err=111) &
1085 (aathet(l,i,j,k,iblock),l=1,ntheterm)
1086 read (ithep,*,end=111,err=111) &
1087 ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1088 (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1089 (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1090 (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1092 read (ithep,*,end=111,err=111) &
1093 (((ffthet(llll,lll,ll,i,j,k,iblock),&
1094 ffthet(lll,llll,ll,i,j,k,iblock),&
1095 ggthet(llll,lll,ll,i,j,k,iblock),&
1096 ggthet(lll,llll,ll,i,j,k,iblock),&
1097 llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
1102 ! For dummy ends assign glycine-type coefficients of theta-only terms; the
1103 ! coefficients of theta-and-gamma-dependent terms are zero.
1104 ! IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
1105 ! RECOMENTDED AFTER VERSION 3.3)
1109 ! aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
1110 ! aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
1112 ! aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
1113 ! aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
1116 ! aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
1118 ! aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
1121 ! AND COMMENT THE LOOPS BELOW
1125 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
1126 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
1128 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
1129 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
1132 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1134 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1139 ! Substitution for D aminoacids from symmetry.
1142 do j=-nthetyp,nthetyp
1143 do k=-nthetyp,nthetyp
1144 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
1146 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
1150 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
1151 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
1152 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
1153 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
1159 ffthet(llll,lll,ll,i,j,k,iblock)= &
1160 ffthet(llll,lll,ll,-i,-j,-k,iblock)
1161 ffthet(lll,llll,ll,i,j,k,iblock)= &
1162 ffthet(lll,llll,ll,-i,-j,-k,iblock)
1163 ggthet(llll,lll,ll,i,j,k,iblock)= &
1164 -ggthet(llll,lll,ll,-i,-j,-k,iblock)
1165 ggthet(lll,llll,ll,i,j,k,iblock)= &
1166 -ggthet(lll,llll,ll,-i,-j,-k,iblock)
1175 ! Control printout of the coefficients of virtual-bond-angle potentials
1178 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
1183 write (iout,'(//4a)') &
1184 'Type ',onelett(i),onelett(j),onelett(k)
1185 write (iout,'(//a,10x,a)') " l","a[l]"
1186 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
1187 write (iout,'(i2,1pe15.5)') &
1188 (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
1190 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))') &
1191 "b",l,"c",l,"d",l,"e",l
1193 write (iout,'(i2,4(1pe15.5))') m,&
1194 bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
1195 ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
1199 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))') &
1200 "f+",l,"f-",l,"g+",l,"g-",l
1203 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
1204 ffthet(n,m,l,i,j,k,iblock),&
1205 ffthet(m,n,l,i,j,k,iblock),&
1206 ggthet(n,m,l,i,j,k,iblock),&
1207 ggthet(m,n,l,i,j,k,iblock)
1218 !C here will be the apropriate recalibrating for D-aminoacid
1219 read (ithep,*,end=121,err=121) nthetyp
1220 allocate(nbend_kcc_Tb(-nthetyp:nthetyp))
1221 allocate(v1bend_chyb(0:36,-nthetyp:nthetyp))
1222 do i=-nthetyp+1,nthetyp-1
1223 read (ithep,*,end=121,err=121) nbend_kcc_Tb(i)
1224 do j=0,nbend_kcc_Tb(i)
1225 read (ithep,*,end=121,err=121) ijunk,v1bend_chyb(j,i)
1229 write (iout,'(a)') &
1230 "Parameters of the valence-only potentials"
1231 do i=-nthetyp+1,nthetyp-1
1232 write (iout,'(2a)') "Type ",toronelet(i)
1233 do k=0,nbend_kcc_Tb(i)
1234 write(iout,'(i5,f15.5)') k,v1bend_chyb(k,i)
1240 write (2,*) "Start reading THETA_PDB",ithep_pdb
1242 ! write (2,*) 'i=',i
1243 read (ithep_pdb,*,err=111,end=111) &
1244 a0thet(i),(athet(j,i,1,1),j=1,2),&
1245 (bthet(j,i,1,1),j=1,2)
1246 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
1247 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
1248 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
1249 sigc0(i)=sigc0(i)**2
1252 athet(1,i,1,-1)=athet(1,i,1,1)
1253 athet(2,i,1,-1)=athet(2,i,1,1)
1254 bthet(1,i,1,-1)=-bthet(1,i,1,1)
1255 bthet(2,i,1,-1)=-bthet(2,i,1,1)
1256 athet(1,i,-1,1)=-athet(1,i,1,1)
1257 athet(2,i,-1,1)=-athet(2,i,1,1)
1258 bthet(1,i,-1,1)=bthet(1,i,1,1)
1259 bthet(2,i,-1,1)=bthet(2,i,1,1)
1262 a0thet(i)=a0thet(-i)
1263 athet(1,i,-1,-1)=athet(1,-i,1,1)
1264 athet(2,i,-1,-1)=-athet(2,-i,1,1)
1265 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
1266 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
1267 athet(1,i,-1,1)=athet(1,-i,1,1)
1268 athet(2,i,-1,1)=-athet(2,-i,1,1)
1269 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
1270 bthet(2,i,-1,1)=bthet(2,-i,1,1)
1271 athet(1,i,1,-1)=-athet(1,-i,1,1)
1272 athet(2,i,1,-1)=athet(2,-i,1,1)
1273 bthet(1,i,1,-1)=bthet(1,-i,1,1)
1274 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
1275 theta0(i)=theta0(-i)
1279 polthet(j,i)=polthet(j,-i)
1282 gthet(j,i)=gthet(j,-i)
1285 write (2,*) "End reading THETA_PDB"
1289 !--------------- Reading theta parameters for nucleic acid-------
1290 read (ithep_nucl,*,err=111,end=111) nthetyp_nucl,ntheterm_nucl,&
1291 ntheterm2_nucl,ntheterm3_nucl,nsingle_nucl,ndouble_nucl
1292 nntheterm_nucl=max0(ntheterm_nucl,ntheterm2_nucl,ntheterm3_nucl)
1293 allocate(ithetyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
1294 allocate(aa0thet_nucl(nthetyp_nucl+1,&
1295 nthetyp_nucl+1,nthetyp_nucl+1))
1296 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1297 allocate(aathet_nucl(ntheterm_nucl+1,nthetyp_nucl+1,&
1298 nthetyp_nucl+1,nthetyp_nucl+1))
1299 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1300 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1301 allocate(bbthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1302 nthetyp_nucl+1,nthetyp_nucl+1))
1303 allocate(ccthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1304 nthetyp_nucl+1,nthetyp_nucl+1))
1305 allocate(ddthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1306 nthetyp_nucl+1,nthetyp_nucl+1))
1307 allocate(eethet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1308 nthetyp_nucl+1,nthetyp_nucl+1))
1309 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1310 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1311 allocate(ffthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1312 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1313 allocate(ggthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1314 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1316 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1317 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1319 read (ithep_nucl,*,err=111,end=111) (ithetyp_nucl(i),i=1,ntyp1_molec(2))
1321 aa0thet_nucl(:,:,:)=0.0d0
1322 aathet_nucl(:,:,:,:)=0.0d0
1323 bbthet_nucl(:,:,:,:,:)=0.0d0
1324 ccthet_nucl(:,:,:,:,:)=0.0d0
1325 ddthet_nucl(:,:,:,:,:)=0.0d0
1326 eethet_nucl(:,:,:,:,:)=0.0d0
1327 ffthet_nucl(:,:,:,:,:,:)=0.0d0
1328 ggthet_nucl(:,:,:,:,:,:)=0.0d0
1333 read (ithep_nucl,'(3a)',end=111,err=111) t1,t2,t3
1334 read (ithep_nucl,*,end=111,err=111) aa0thet_nucl(i,j,k)
1335 read (ithep_nucl,*,end=111,err=111)(aathet_nucl(l,i,j,k),l=1,ntheterm_nucl)
1336 read (ithep_nucl,*,end=111,err=111) &
1337 (((bbthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1338 (ccthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1339 (ddthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1340 (eethet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl)),ll=1,ntheterm2_nucl)
1341 read (ithep_nucl,*,end=111,err=111) &
1342 (((ffthet_nucl(llll,lll,ll,i,j,k),ffthet_nucl(lll,llll,ll,i,j,k), &
1343 ggthet_nucl(llll,lll,ll,i,j,k),ggthet_nucl(lll,llll,ll,i,j,k), &
1344 llll=1,lll-1),lll=2,ndouble_nucl),ll=1,ntheterm3_nucl)
1349 !-------------------------------------------
1350 allocate(nlob(ntyp1)) !(ntyp1)
1351 allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
1352 allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
1353 allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
1360 gaussc(:,:,:,:)=0.0D0
1364 ! Read the parameters of the probability distribution/energy expression
1365 ! of the side chains.
1368 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1372 dsc_inv(i)=1.0D0/dsc(i)
1383 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1384 ((blower(k,l,1),l=1,k),k=1,3)
1385 censc(1,1,-i)=censc(1,1,i)
1386 censc(2,1,-i)=censc(2,1,i)
1387 censc(3,1,-i)=-censc(3,1,i)
1389 read (irotam,*,end=112,err=112) bsc(j,i)
1390 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1391 ((blower(k,l,j),l=1,k),k=1,3)
1392 censc(1,j,-i)=censc(1,j,i)
1393 censc(2,j,-i)=censc(2,j,i)
1394 censc(3,j,-i)=-censc(3,j,i)
1395 ! BSC is amplitude of Gaussian
1402 akl=akl+blower(k,m,j)*blower(l,m,j)
1406 if (((k.eq.3).and.(l.ne.3)) &
1407 .or.((l.eq.3).and.(k.ne.3))) then
1408 gaussc(k,l,j,-i)=-akl
1409 gaussc(l,k,j,-i)=-akl
1411 gaussc(k,l,j,-i)=akl
1412 gaussc(l,k,j,-i)=akl
1421 write (iout,'(/a)') 'Parameters of side-chain local geometry'
1424 if (nlobi.gt.0) then
1426 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i,1),&
1427 ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
1428 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
1429 'log h',(bsc(j,i),j=1,nlobi)
1430 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
1431 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
1433 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
1434 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
1437 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
1438 write (iout,'(a,f10.4,4(16x,f10.4))') &
1439 'Center ',(bsc(j,i),j=1,nlobi)
1440 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),&
1449 ! Read scrot parameters for potentials determined from all-atom AM1 calculations
1450 ! added by Urszula Kozlowska 07/11/2007
1452 !el Maximum number of SC local term fitting function coefficiants
1453 !el integer,parameter :: maxsccoef=65
1455 allocate(sc_parmin(65,ntyp)) !(maxsccoef,ntyp)
1458 read (irotam,*,end=112,err=112)
1460 read (irotam,*,end=112,err=112)
1463 read(irotam,*,end=112,err=112) sc_parmin(j,i)
1467 !---------reading nucleic acid parameters for rotamers-------------------
1468 allocate(sc_parmin_nucl(9,ntyp_molec(2))) !(maxsccoef,ntyp)
1469 do i=1,ntyp_molec(2)
1470 read (irotam_nucl,*,end=112,err=112)
1472 read(irotam_nucl,*,end=112,err=112) sc_parmin_nucl(j,i)
1478 write (iout,*) "Base rotamer parameters"
1479 do i=1,ntyp_molec(2)
1480 write (iout,'(a)') restyp(i,2)
1481 write (iout,'(i5,f10.5)') (i,sc_parmin_nucl(j,i),j=1,9)
1486 ! Read the parameters of the probability distribution/energy expression
1487 ! of the side chains.
1489 write (2,*) "Start reading ROTAM_PDB"
1491 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1495 dsc_inv(i)=1.0D0/dsc(i)
1506 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1507 ((blower(k,l,1),l=1,k),k=1,3)
1509 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
1510 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1511 ((blower(k,l,j),l=1,k),k=1,3)
1518 akl=akl+blower(k,m,j)*blower(l,m,j)
1528 write (2,*) "End reading ROTAM_PDB"
1534 !C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
1535 !C interaction energy of the Gly, Ala, and Pro prototypes.
1537 read (ifourier,*) nloctyp
1538 SPLIT_FOURIERTOR = nloctyp.lt.0
1539 nloctyp = iabs(nloctyp)
1540 !C allocate(b1(2,nres)) !(2,-maxtor:maxtor)
1541 !C allocate(b2(2,nres)) !(2,-maxtor:maxtor)
1542 !C allocate(b1tilde(2,nres)) !(2,-maxtor:maxtor)
1543 !C allocate(ctilde(2,2,nres))
1544 !C allocate(dtilde(2,2,nres)) !(2,2,-maxtor:maxtor)
1545 !C allocate(gtb1(2,nres))
1546 !C allocate(gtb2(2,nres))
1547 !C allocate(cc(2,2,nres))
1548 !C allocate(dd(2,2,nres))
1549 !C allocate(ee(2,2,nres))
1550 !C allocate(gtcc(2,2,nres))
1551 !C allocate(gtdd(2,2,nres))
1552 !C allocate(gtee(2,2,nres))
1555 allocate(itype2loc(-ntyp1:ntyp1))
1556 allocate(iloctyp(-nloctyp:nloctyp))
1557 allocate(bnew1(3,2,-nloctyp:nloctyp))
1558 allocate(bnew2(3,2,-nloctyp:nloctyp))
1559 allocate(ccnew(3,2,-nloctyp:nloctyp))
1560 allocate(ddnew(3,2,-nloctyp:nloctyp))
1561 allocate(e0new(3,-nloctyp:nloctyp))
1562 allocate(eenew(2,2,2,-nloctyp:nloctyp))
1563 allocate(bnew1tor(3,2,-nloctyp:nloctyp))
1564 allocate(bnew2tor(3,2,-nloctyp:nloctyp))
1565 allocate(ccnewtor(3,2,-nloctyp:nloctyp))
1566 allocate(ddnewtor(3,2,-nloctyp:nloctyp))
1567 allocate(e0newtor(3,-nloctyp:nloctyp))
1568 allocate(eenewtor(2,2,2,-nloctyp:nloctyp))
1570 read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp)
1571 read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1)
1572 itype2loc(ntyp1)=nloctyp
1573 iloctyp(nloctyp)=ntyp1
1575 itype2loc(-i)=-itype2loc(i)
1578 allocate(iloctyp(-nloctyp:nloctyp))
1579 allocate(itype2loc(-ntyp1:ntyp1))
1586 iloctyp(-i)=-iloctyp(i)
1588 !c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1589 !c write (iout,*) "nloctyp",nloctyp,
1590 !c & " iloctyp",(iloctyp(i),i=0,nloctyp)
1591 !c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1592 !c write (iout,*) "nloctyp",nloctyp,
1593 !c & " iloctyp",(iloctyp(i),i=0,nloctyp)
1596 !c write (iout,*) "NEWCORR",i
1597 read (ifourier,*,end=115,err=115)
1600 read (ifourier,*,end=115,err=115) bnew1(ii,j,i)
1603 !c write (iout,*) "NEWCORR BNEW1"
1604 !c write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2)
1607 read (ifourier,*,end=115,err=115) bnew2(ii,j,i)
1610 !c write (iout,*) "NEWCORR BNEW2"
1611 !c write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2)
1613 read (ifourier,*,end=115,err=115) ccnew(kk,1,i)
1614 read (ifourier,*,end=115,err=115) ccnew(kk,2,i)
1616 !c write (iout,*) "NEWCORR CCNEW"
1617 !c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
1619 read (ifourier,*,end=115,err=115) ddnew(kk,1,i)
1620 read (ifourier,*,end=115,err=115) ddnew(kk,2,i)
1622 !c write (iout,*) "NEWCORR DDNEW"
1623 !c write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2)
1627 read (ifourier,*,end=115,err=115) eenew(ii,jj,kk,i)
1631 !c write (iout,*) "NEWCORR EENEW1"
1632 !c write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
1634 read (ifourier,*,end=115,err=115) e0new(ii,i)
1636 !c write (iout,*) (e0new(ii,i),ii=1,3)
1638 !c write (iout,*) "NEWCORR EENEW"
1641 ccnew(ii,1,i)=ccnew(ii,1,i)/2
1642 ccnew(ii,2,i)=ccnew(ii,2,i)/2
1643 ddnew(ii,1,i)=ddnew(ii,1,i)/2
1644 ddnew(ii,2,i)=ddnew(ii,2,i)/2
1649 bnew1(ii,1,-i)= bnew1(ii,1,i)
1650 bnew1(ii,2,-i)=-bnew1(ii,2,i)
1651 bnew2(ii,1,-i)= bnew2(ii,1,i)
1652 bnew2(ii,2,-i)=-bnew2(ii,2,i)
1655 !c ccnew(ii,1,i)=ccnew(ii,1,i)/2
1656 !c ccnew(ii,2,i)=ccnew(ii,2,i)/2
1657 !c ddnew(ii,1,i)=ddnew(ii,1,i)/2
1658 !c ddnew(ii,2,i)=ddnew(ii,2,i)/2
1659 ccnew(ii,1,-i)=ccnew(ii,1,i)
1660 ccnew(ii,2,-i)=-ccnew(ii,2,i)
1661 ddnew(ii,1,-i)=ddnew(ii,1,i)
1662 ddnew(ii,2,-i)=-ddnew(ii,2,i)
1664 e0new(1,-i)= e0new(1,i)
1665 e0new(2,-i)=-e0new(2,i)
1666 e0new(3,-i)=-e0new(3,i)
1668 eenew(kk,1,1,-i)= eenew(kk,1,1,i)
1669 eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
1670 eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
1671 eenew(kk,2,2,-i)= eenew(kk,2,2,i)
1675 write (iout,'(a)') "Coefficients of the multibody terms"
1676 do i=-nloctyp+1,nloctyp-1
1677 write (iout,*) "Type: ",onelet(iloctyp(i))
1678 write (iout,*) "Coefficients of the expansion of B1"
1680 write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
1682 write (iout,*) "Coefficients of the expansion of B2"
1684 write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
1686 write (iout,*) "Coefficients of the expansion of C"
1687 write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
1688 write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
1689 write (iout,*) "Coefficients of the expansion of D"
1690 write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
1691 write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
1692 write (iout,*) "Coefficients of the expansion of E"
1693 write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
1696 write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
1701 IF (SPLIT_FOURIERTOR) THEN
1703 !c write (iout,*) "NEWCORR TOR",i
1704 read (ifourier,*,end=115,err=115)
1707 read (ifourier,*,end=115,err=115) bnew1tor(ii,j,i)
1710 !c write (iout,*) "NEWCORR BNEW1 TOR"
1711 !c write (iout,*) ((bnew1tor(ii,j,i),ii=1,3),j=1,2)
1714 read (ifourier,*,end=115,err=115) bnew2tor(ii,j,i)
1717 !c write (iout,*) "NEWCORR BNEW2 TOR"
1718 !c write (iout,*) ((bnew2tor(ii,j,i),ii=1,3),j=1,2)
1720 read (ifourier,*,end=115,err=115) ccnewtor(kk,1,i)
1721 read (ifourier,*,end=115,err=115) ccnewtor(kk,2,i)
1723 !c write (iout,*) "NEWCORR CCNEW TOR"
1724 !c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
1726 read (ifourier,*,end=115,err=115) ddnewtor(kk,1,i)
1727 read (ifourier,*,end=115,err=115) ddnewtor(kk,2,i)
1729 !c write (iout,*) "NEWCORR DDNEW TOR"
1730 !c write (iout,*) ((ddnewtor(ii,j,i),ii=1,3),j=1,2)
1734 read (ifourier,*,end=115,err=115) eenewtor(ii,jj,kk,i)
1738 !c write (iout,*) "NEWCORR EENEW1 TOR"
1739 !c write(iout,*)(((eenewtor(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
1741 read (ifourier,*,end=115,err=115) e0newtor(ii,i)
1743 !c write (iout,*) (e0newtor(ii,i),ii=1,3)
1745 !c write (iout,*) "NEWCORR EENEW TOR"
1748 ccnewtor(ii,1,i)=ccnewtor(ii,1,i)/2
1749 ccnewtor(ii,2,i)=ccnewtor(ii,2,i)/2
1750 ddnewtor(ii,1,i)=ddnewtor(ii,1,i)/2
1751 ddnewtor(ii,2,i)=ddnewtor(ii,2,i)/2
1756 bnew1tor(ii,1,-i)= bnew1tor(ii,1,i)
1757 bnew1tor(ii,2,-i)=-bnew1tor(ii,2,i)
1758 bnew2tor(ii,1,-i)= bnew2tor(ii,1,i)
1759 bnew2tor(ii,2,-i)=-bnew2tor(ii,2,i)
1762 ccnewtor(ii,1,-i)=ccnewtor(ii,1,i)
1763 ccnewtor(ii,2,-i)=-ccnewtor(ii,2,i)
1764 ddnewtor(ii,1,-i)=ddnewtor(ii,1,i)
1765 ddnewtor(ii,2,-i)=-ddnewtor(ii,2,i)
1767 e0newtor(1,-i)= e0newtor(1,i)
1768 e0newtor(2,-i)=-e0newtor(2,i)
1769 e0newtor(3,-i)=-e0newtor(3,i)
1771 eenewtor(kk,1,1,-i)= eenewtor(kk,1,1,i)
1772 eenewtor(kk,1,2,-i)=-eenewtor(kk,1,2,i)
1773 eenewtor(kk,2,1,-i)=-eenewtor(kk,2,1,i)
1774 eenewtor(kk,2,2,-i)= eenewtor(kk,2,2,i)
1778 write (iout,'(a)') &
1779 "Single-body coefficients of the torsional potentials"
1780 do i=-nloctyp+1,nloctyp-1
1781 write (iout,*) "Type: ",onelet(iloctyp(i))
1782 write (iout,*) "Coefficients of the expansion of B1tor"
1784 write (iout,'(3hB1(,i1,1h),3f10.5)') &
1785 j,(bnew1tor(k,j,i),k=1,3)
1787 write (iout,*) "Coefficients of the expansion of B2tor"
1789 write (iout,'(3hB2(,i1,1h),3f10.5)') &
1790 j,(bnew2tor(k,j,i),k=1,3)
1792 write (iout,*) "Coefficients of the expansion of Ctor"
1793 write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3)
1794 write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3)
1795 write (iout,*) "Coefficients of the expansion of Dtor"
1796 write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3)
1797 write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3)
1798 write (iout,*) "Coefficients of the expansion of Etor"
1799 write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3)
1802 write (iout,'(1hE,2i1,2f10.5)') &
1803 j,k,(eenewtor(l,j,k,i),l=1,2)
1809 do i=-nloctyp+1,nloctyp-1
1812 bnew1tor(ii,j,i)=bnew1(ii,j,i)
1817 bnew2tor(ii,j,i)=bnew2(ii,j,i)
1821 ccnewtor(ii,1,i)=ccnew(ii,1,i)
1822 ccnewtor(ii,2,i)=ccnew(ii,2,i)
1823 ddnewtor(ii,1,i)=ddnew(ii,1,i)
1824 ddnewtor(ii,2,i)=ddnew(ii,2,i)
1827 ENDIF !SPLIT_FOURIER_TOR
1829 allocate(ccold(2,2,-nloctyp-1:nloctyp+1))
1830 allocate(ddold(2,2,-nloctyp-1:nloctyp+1))
1831 allocate(eeold(2,2,-nloctyp-1:nloctyp+1))
1832 allocate(b(13,-nloctyp-1:nloctyp+1))
1834 write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)"
1836 read (ifourier,*,end=115,err=115)
1837 read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
1839 write (iout,*) 'Type ',onelet(iloctyp(i))
1840 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
1854 !c B1tilde(1,i) = b(3)
1855 !c! B1tilde(2,i) =-b(5)
1856 !c! B1tilde(1,-i) =-b(3)
1857 !c! B1tilde(2,-i) =b(5)
1858 !c! b1tilde(1,i)=0.0d0
1859 !c b1tilde(2,i)=0.0d0
1864 !cc B1tilde(1,i) = b(3,i)
1865 !cc B1tilde(2,i) =-b(5,i)
1866 !c B1tilde(1,-i) =-b(3,i)
1867 !c B1tilde(2,-i) =b(5,i)
1868 !cc b1tilde(1,i)=0.0d0
1869 !cc b1tilde(2,i)=0.0d0
1870 !cc B2(1,i) = b(2,i)
1871 !cc B2(2,i) = b(4,i)
1873 !c B2(2,-i) =-b(4,i)
1877 CCold(1,1,i)= b(7,i)
1878 CCold(2,2,i)=-b(7,i)
1879 CCold(2,1,i)= b(9,i)
1880 CCold(1,2,i)= b(9,i)
1881 CCold(1,1,-i)= b(7,i)
1882 CCold(2,2,-i)=-b(7,i)
1883 CCold(2,1,-i)=-b(9,i)
1884 CCold(1,2,-i)=-b(9,i)
1889 !c Ctilde(1,1,i)= CCold(1,1,i)
1890 !c Ctilde(1,2,i)= CCold(1,2,i)
1891 !c Ctilde(2,1,i)=-CCold(2,1,i)
1892 !c Ctilde(2,2,i)=-CCold(2,2,i)
1897 !c Ctilde(1,1,i)= CCold(1,1,i)
1898 !c Ctilde(1,2,i)= CCold(1,2,i)
1899 !c Ctilde(2,1,i)=-CCold(2,1,i)
1900 !c Ctilde(2,2,i)=-CCold(2,2,i)
1902 !c Ctilde(1,1,i)=0.0d0
1903 !c Ctilde(1,2,i)=0.0d0
1904 !c Ctilde(2,1,i)=0.0d0
1905 !c Ctilde(2,2,i)=0.0d0
1906 DDold(1,1,i)= b(6,i)
1907 DDold(2,2,i)=-b(6,i)
1908 DDold(2,1,i)= b(8,i)
1909 DDold(1,2,i)= b(8,i)
1910 DDold(1,1,-i)= b(6,i)
1911 DDold(2,2,-i)=-b(6,i)
1912 DDold(2,1,-i)=-b(8,i)
1913 DDold(1,2,-i)=-b(8,i)
1918 !c Dtilde(1,1,i)= DD(1,1,i)
1919 !c Dtilde(1,2,i)= DD(1,2,i)
1920 !c Dtilde(2,1,i)=-DD(2,1,i)
1921 !c Dtilde(2,2,i)=-DD(2,2,i)
1923 !c Dtilde(1,1,i)=0.0d0
1924 !c Dtilde(1,2,i)=0.0d0
1925 !c Dtilde(2,1,i)=0.0d0
1926 !c Dtilde(2,2,i)=0.0d0
1927 EEold(1,1,i)= b(10,i)+b(11,i)
1928 EEold(2,2,i)=-b(10,i)+b(11,i)
1929 EEold(2,1,i)= b(12,i)-b(13,i)
1930 EEold(1,2,i)= b(12,i)+b(13,i)
1931 EEold(1,1,-i)= b(10,i)+b(11,i)
1932 EEold(2,2,-i)=-b(10,i)+b(11,i)
1933 EEold(2,1,-i)=-b(12,i)+b(13,i)
1934 EEold(1,2,-i)=-b(12,i)-b(13,i)
1935 write(iout,*) "TU DOCHODZE"
1941 !c ee(2,1,i)=ee(1,2,i)
1946 "Coefficients of the cumulants (independent of valence angles)"
1947 do i=-nloctyp+1,nloctyp-1
1948 write (iout,*) 'Type ',onelet(iloctyp(i))
1950 write(iout,'(2f10.5)') B(3,i),B(5,i)
1952 write(iout,'(2f10.5)') B(2,i),B(4,i)
1955 write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
1959 write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
1963 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1972 ! Read torsional parameters in old format
1974 allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
1976 read (itorp,*,end=113,err=113) ntortyp,nterm_old
1977 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1978 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1980 !el from energy module--------
1981 allocate(v1(nterm_old,ntortyp,ntortyp))
1982 allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
1983 !el---------------------------
1988 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
1994 write (iout,'(/a/)') 'Torsional constants:'
1997 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1998 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
2004 ! Read torsional parameters
2006 IF (TOR_MODE.eq.0) THEN
2007 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
2008 read (itorp,*,end=113,err=113) ntortyp
2009 !el from energy module---------
2010 allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2011 allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2013 allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
2014 allocate(vlor2(maxlor,ntortyp,ntortyp))
2015 allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
2016 allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2018 allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
2019 allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
2020 !el---------------------------
2023 !el---------------------------
2025 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
2027 itortyp(i)=-itortyp(-i)
2029 itortyp(ntyp1)=ntortyp
2030 itortyp(-ntyp1)=-ntortyp
2032 write (iout,*) 'ntortyp',ntortyp
2034 do j=-ntortyp+1,ntortyp-1
2035 read (itorp,*,end=113,err=113) nterm(i,j,iblock),&
2037 nterm(-i,-j,iblock)=nterm(i,j,iblock)
2038 nlor(-i,-j,iblock)=nlor(i,j,iblock)
2041 do k=1,nterm(i,j,iblock)
2042 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),&
2044 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
2045 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
2046 v0ij=v0ij+si*v1(k,i,j,iblock)
2048 ! write(iout,*) i,j,k,iblock,nterm(i,j,iblock) !
2049 ! write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),&!
2050 ! v2(k,-i,-j,iblock),v2(k,i,j,iblock)!
2052 do k=1,nlor(i,j,iblock)
2053 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),&
2054 vlor2(k,i,j),vlor3(k,i,j)
2055 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
2058 v0(-i,-j,iblock)=v0ij
2064 write (iout,'(/a/)') 'Torsional constants:'
2066 do i=-ntortyp,ntortyp
2067 do j=-ntortyp,ntortyp
2068 write (iout,*) 'ityp',i,' jtyp',j
2069 write (iout,*) 'Fourier constants'
2070 do k=1,nterm(i,j,iblock)
2071 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
2074 write (iout,*) 'Lorenz constants'
2075 do k=1,nlor(i,j,iblock)
2076 write (iout,'(3(1pe15.5))') &
2077 vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
2083 !elwrite (iout,'(/a/)') 'Torsional constants:',vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
2085 ! 6/23/01 Read parameters for double torsionals
2087 !el from energy module------------
2088 allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2089 allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2090 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
2091 allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2092 allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2093 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
2094 allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2095 allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2096 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
2097 !---------------------------------
2101 do j=-ntortyp+1,ntortyp-1
2102 do k=-ntortyp+1,ntortyp-1
2103 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
2104 ! write (iout,*) "OK onelett",
2107 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
2108 .or. t3.ne.toronelet(k)) then
2109 write (iout,*) "Error in double torsional parameter file",&
2112 call MPI_Finalize(Ierror)
2114 stop "Error in double torsional parameter file"
2116 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),&
2117 ntermd_2(i,j,k,iblock)
2118 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
2119 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
2120 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,&
2121 ntermd_1(i,j,k,iblock))
2122 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,&
2123 ntermd_1(i,j,k,iblock))
2124 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,&
2125 ntermd_1(i,j,k,iblock))
2126 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,&
2127 ntermd_1(i,j,k,iblock))
2128 ! Martix of D parameters for one dimesional foureir series
2129 do l=1,ntermd_1(i,j,k,iblock)
2130 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
2131 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
2132 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
2133 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
2134 ! write(iout,*) "whcodze" ,
2135 ! & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
2137 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),&
2138 v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
2139 v2s(m,l,i,j,k,iblock),&
2140 m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
2141 ! Martix of D parameters for two dimesional fourier series
2142 do l=1,ntermd_2(i,j,k,iblock)
2144 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
2145 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
2146 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
2147 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
2156 write (iout,*) 'Constants for double torsionals'
2159 do j=-ntortyp+1,ntortyp-1
2160 do k=-ntortyp+1,ntortyp-1
2161 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
2162 ' nsingle',ntermd_1(i,j,k,iblock),&
2163 ' ndouble',ntermd_2(i,j,k,iblock)
2165 write (iout,*) 'Single angles:'
2166 do l=1,ntermd_1(i,j,k,iblock)
2167 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
2168 v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
2169 v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
2170 v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
2173 write (iout,*) 'Pairs of angles:'
2174 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
2175 do l=1,ntermd_2(i,j,k,iblock)
2176 write (iout,'(i5,20f10.5)') &
2177 l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
2180 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
2181 do l=1,ntermd_2(i,j,k,iblock)
2182 write (iout,'(i5,20f10.5)') &
2183 l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
2184 (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
2194 itype2loc(i)=itortyp(i)
2198 ELSE IF (TOR_MODE.eq.1) THEN
2200 !C read valence-torsional parameters
2201 read (itorp,*,end=121,err=121) ntortyp
2203 write (iout,*) "Valence-torsional parameters read in ntortyp",&
2205 read (itorp,*,end=121,err=121) (itortyp(i),i=1,ntyp)
2206 write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp)
2209 itype2loc(i)=itortyp(i)
2213 itortyp(i)=-itortyp(-i)
2215 do i=-ntortyp+1,ntortyp-1
2216 do j=-ntortyp+1,ntortyp-1
2217 !C first we read the cos and sin gamma parameters
2218 read (itorp,'(13x,a)',end=121,err=121) string
2219 write (iout,*) i,j,string
2220 read (itorp,*,end=121,err=121) &
2221 nterm_kcc(j,i),nterm_kcc_Tb(j,i)
2222 !C read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
2223 do k=1,nterm_kcc(j,i)
2224 do l=1,nterm_kcc_Tb(j,i)
2225 do ll=1,nterm_kcc_Tb(j,i)
2226 read (itorp,*,end=121,err=121) ii,jj,kk, &
2227 v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
2235 !c AL 4/8/16: Calculate coefficients from one-body parameters
2237 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
2238 allocate(nterm_kcc(-ntyp1:ntyp1,-ntyp1:ntyp1))
2239 allocate(nterm_kcc_Tb(-ntyp1:ntyp1,-ntyp1:ntyp1))
2240 allocate(v1_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1))
2241 allocate(v2_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1))
2244 print *,i,itortyp(i)
2245 itortyp(i)=itype2loc(i)
2248 "Val-tor parameters calculated from cumulant coefficients ntortyp"&
2250 do i=-ntortyp+1,ntortyp-1
2251 do j=-ntortyp+1,ntortyp-1
2254 do k=1,nterm_kcc_Tb(j,i)
2255 do l=1,nterm_kcc_Tb(j,i)
2256 v1_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,1,j)&
2257 +bnew1tor(k,2,i)*bnew2tor(l,2,j)
2258 v2_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,2,j)&
2259 +bnew1tor(k,2,i)*bnew2tor(l,1,j)
2262 do k=1,nterm_kcc_Tb(j,i)
2263 do l=1,nterm_kcc_Tb(j,i)
2265 v1_kcc(k,l,2,i,j)=-(ccnewtor(k,1,i)*ddnewtor(l,1,j) &
2266 -ccnewtor(k,2,i)*ddnewtor(l,2,j))
2267 v2_kcc(k,l,2,i,j)=-(ccnewtor(k,2,i)*ddnewtor(l,1,j) &
2268 +ccnewtor(k,1,i)*ddnewtor(l,2,j))
2270 v1_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,1,i)*ddnewtor(l,1,j) &
2271 -ccnewtor(k,2,i)*ddnewtor(l,2,j))
2272 v2_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,2,i)*ddnewtor(l,1,j) &
2273 +ccnewtor(k,1,i)*ddnewtor(l,2,j))
2277 !c f(theta,gamma)=-(b21(theta)*b11(theta)+b12(theta)*b22(theta))*cos(gamma)-(b22(theta)*b11(theta)+b21(theta)*b12(theta))*sin(gamma)+(c11(theta)*d11(theta)-c12(theta)*d12(theta))*cos(2*gamma)+(c12(theta)*d11(theta)+c11(theta)*d12(theta))*sin(2*gamma)
2281 write (iout,*) "TOR_MODE>1 only with NEWCORR"
2286 if (tor_mode.gt.0 .and. lprint) then
2287 !c Print valence-torsional parameters
2288 write (iout,'(a)') &
2289 "Parameters of the valence-torsional potentials"
2290 do i=-ntortyp+1,ntortyp-1
2291 do j=-ntortyp+1,ntortyp-1
2292 write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
2293 write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
2294 do k=1,nterm_kcc(j,i)
2295 do l=1,nterm_kcc_Tb(j,i)
2296 do ll=1,nterm_kcc_Tb(j,i)
2297 write (iout,'(3i5,2f15.4)')&
2298 k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
2306 allocate(itortyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
2307 read (itorp_nucl,*,end=113,err=113) ntortyp_nucl
2308 ! print *,"ntortyp_nucl",ntortyp_nucl,ntyp_molec(2)
2309 !el from energy module---------
2310 allocate(nterm_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2311 allocate(nlor_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2313 allocate(vlor1_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
2314 allocate(vlor2_nucl(maxlor,ntortyp_nucl,ntortyp_nucl))
2315 allocate(vlor3_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,maxtor,maxtor)
2316 allocate(v0_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2318 allocate(v1_nucl(maxterm,ntortyp_nucl,ntortyp_nucl))
2319 allocate(v2_nucl(maxterm,ntortyp_nucl,ntortyp_nucl)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
2320 !el---------------------------
2323 !el--------------------
2324 read (itorp_nucl,*,end=113,err=113) &
2325 (itortyp_nucl(i),i=1,ntyp_molec(2))
2326 ! print *,itortyp_nucl(:)
2327 !c write (iout,*) 'ntortyp',ntortyp
2330 read (itorp_nucl,*,end=113,err=113) nterm_nucl(i,j),nlor_nucl(i,j)
2331 ! print *,nterm_nucl(i,j),nlor_nucl(i,j)
2334 do k=1,nterm_nucl(i,j)
2335 read (itorp_nucl,*,end=113,err=113) kk,v1_nucl(k,i,j),v2_nucl(k,i,j)
2336 v0ij=v0ij+si*v1_nucl(k,i,j)
2339 do k=1,nlor_nucl(i,j)
2340 read (itorp,*,end=113,err=113) kk,vlor1_nucl(k,i,j),&
2341 vlor2_nucl(k,i,j),vlor3_nucl(k,i,j)
2342 v0ij=v0ij+vlor1_nucl(k,i,j)/(1+vlor3_nucl(k,i,j)**2)
2348 ! Read of Side-chain backbone correlation parameters
2349 ! Modified 11 May 2012 by Adasko
2352 read (isccor,*,end=119,err=119) nsccortyp
2354 !el from module energy-------------
2355 allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2356 allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
2357 allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
2358 allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp)) !(maxterm_sccor,20,20)
2359 !-----------------------------------
2361 !el from module energy-------------
2362 allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
2364 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
2366 isccortyp(i)=-isccortyp(-i)
2368 iscprol=isccortyp(20)
2369 ! write (iout,*) 'ntortyp',ntortyp
2371 !c maxinter is maximum interaction sites
2372 !el from module energy---------
2373 allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2374 allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
2375 -nsccortyp:nsccortyp))
2376 allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
2377 -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
2378 allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
2379 -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
2380 !-----------------------------------
2384 read (isccor,*,end=119,err=119) &
2385 nterm_sccor(i,j),nlor_sccor(i,j)
2391 nterm_sccor(-i,j)=nterm_sccor(i,j)
2392 nterm_sccor(-i,-j)=nterm_sccor(i,j)
2393 nterm_sccor(i,-j)=nterm_sccor(i,j)
2394 do k=1,nterm_sccor(i,j)
2395 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
2397 if (j.eq.iscprol) then
2398 if (i.eq.isccortyp(10)) then
2399 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
2400 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2402 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
2403 +v2sccor(k,l,i,j)*dsqrt(0.75d0)
2404 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
2405 +v1sccor(k,l,i,j)*dsqrt(0.75d0)
2406 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
2407 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
2408 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
2409 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
2412 if (i.eq.isccortyp(10)) then
2413 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
2414 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2416 if (j.eq.isccortyp(10)) then
2417 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
2418 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
2420 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
2421 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2422 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
2423 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
2424 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
2425 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
2429 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
2430 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
2431 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
2432 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
2435 do k=1,nlor_sccor(i,j)
2436 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
2437 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2438 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
2439 (1+vlor3sccor(k,i,j)**2)
2441 v0sccor(l,i,j)=v0ijsccor
2442 v0sccor(l,-i,j)=v0ijsccor1
2443 v0sccor(l,i,-j)=v0ijsccor2
2444 v0sccor(l,-i,-j)=v0ijsccor3
2450 !el from module energy-------------
2451 allocate(isccortyp(ntyp)) !(-ntyp:ntyp)
2453 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
2454 ! write (iout,*) 'ntortyp',ntortyp
2456 !c maxinter is maximum interaction sites
2457 !el from module energy---------
2458 allocate(nterm_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2459 allocate(v1sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp))
2460 allocate(v2sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
2461 allocate(v0sccor(maxinter,nsccortyp,nsccortyp)) !???(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
2462 !-----------------------------------
2466 read (isccor,*,end=119,err=119) &
2467 nterm_sccor(i,j),nlor_sccor(i,j)
2471 do k=1,nterm_sccor(i,j)
2472 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
2474 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
2477 do k=1,nlor_sccor(i,j)
2478 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
2479 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2480 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
2481 (1+vlor3sccor(k,i,j)**2)
2483 v0sccor(l,i,j)=v0ijsccor !el ,iblock
2492 write (iout,'(/a/)') 'Torsional constants:'
2495 write (iout,*) 'ityp',i,' jtyp',j
2496 write (iout,*) 'Fourier constants'
2497 do k=1,nterm_sccor(i,j)
2498 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
2500 write (iout,*) 'Lorenz constants'
2501 do k=1,nlor_sccor(i,j)
2502 write (iout,'(3(1pe15.5))') &
2503 vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2510 ! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
2511 ! interaction energy of the Gly, Ala, and Pro prototypes.
2514 ! Read electrostatic-interaction parameters
2519 write (iout,'(/a)') 'Electrostatic interaction constants:'
2520 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
2521 'IT','JT','APP','BPP','AEL6','AEL3'
2523 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
2524 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
2525 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
2526 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
2531 app (i,j)=epp(i,j)*rri*rri
2532 bpp (i,j)=-2.0D0*epp(i,j)*rri
2533 ael6(i,j)=elpp6(i,j)*4.2D0**6
2534 ael3(i,j)=elpp3(i,j)*4.2D0**3
2536 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
2542 ! Read side-chain interaction parameters.
2544 !el from module energy - COMMON.INTERACT-------
2545 allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
2546 allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
2547 allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
2549 allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
2550 allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
2551 allocate(epslip(ntyp,ntyp))
2559 !--------------------------------
2561 read (isidep,*,end=117,err=117) ipot,expon
2562 if (ipot.lt.1 .or. ipot.gt.5) then
2563 write (iout,'(2a)') 'Error while reading SC interaction',&
2564 'potential file - unknown potential type.'
2566 call MPI_Finalize(Ierror)
2572 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
2573 ', exponents are ',expon,2*expon
2574 ! goto (10,20,30,30,40) ipot
2576 !----------------------- LJ potential ---------------------------------
2578 ! 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2579 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2580 (sigma0(i),i=1,ntyp)
2582 write (iout,'(/a/)') 'Parameters of the LJ potential:'
2583 write (iout,'(a/)') 'The epsilon array:'
2584 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2585 write (iout,'(/a)') 'One-body parameters:'
2586 write (iout,'(a,4x,a)') 'residue','sigma'
2587 write (iout,'(a3,6x,f10.5)') (restyp(i,1),sigma0(i),i=1,ntyp)
2590 !----------------------- LJK potential --------------------------------
2592 ! 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2593 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2594 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
2596 write (iout,'(/a/)') 'Parameters of the LJK potential:'
2597 write (iout,'(a/)') 'The epsilon array:'
2598 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2599 write (iout,'(/a)') 'One-body parameters:'
2600 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
2601 write (iout,'(a3,6x,2f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2605 !---------------------- GB or BP potential -----------------------------
2608 ! print *,"I AM in SCELE",scelemode
2609 if (scelemode.eq.0) then
2611 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
2613 read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
2614 read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
2615 read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
2616 read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
2618 read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
2621 ! For the GB potential convert sigma'**2 into chi'
2624 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
2628 write (iout,'(/a/)') 'Parameters of the BP potential:'
2629 write (iout,'(a/)') 'The epsilon array:'
2630 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2631 write (iout,'(/a)') 'One-body parameters:'
2632 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',&
2634 write (iout,'(a3,6x,4f10.5)') (restyp(i,1),sigma0(i),sigii(i),&
2635 chip(i),alp(i),i=1,ntyp)
2638 ! print *,ntyp,"NTYP"
2639 allocate(icharge(ntyp1))
2640 ! print *,ntyp,icharge(i)
2642 read (isidep,*) (icharge(i),i=1,ntyp)
2643 print *,ntyp,icharge(i)
2644 ! if(.not.allocated(eps)) allocate(eps(-ntyp
2645 !c write (2,*) "icharge",(icharge(i),i=1,ntyp)
2646 allocate(alphapol(ntyp,ntyp),epshead(ntyp,ntyp),sig0head(ntyp,ntyp))
2647 allocate(sigiso1(ntyp,ntyp),rborn(ntyp,ntyp),sigmap1(ntyp,ntyp))
2648 allocate(sigmap2(ntyp,ntyp),sigiso2(ntyp,ntyp))
2649 allocate(chis(ntyp,ntyp),wquad(ntyp,ntyp),chipp(ntyp,ntyp))
2650 allocate(epsintab(ntyp,ntyp))
2651 allocate(dtail(2,ntyp,ntyp))
2652 allocate(alphasur(4,ntyp,ntyp),alphiso(4,ntyp,ntyp))
2653 allocate(wqdip(2,ntyp,ntyp))
2654 allocate(wstate(4,ntyp,ntyp))
2655 allocate(dhead(2,2,ntyp,ntyp))
2656 allocate(nstate(ntyp,ntyp))
2657 allocate(debaykap(ntyp,ntyp))
2659 if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2660 if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2664 ! write (*,*) "Im in ALAB", i, " ", j
2666 eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), & !6 w tej linii
2667 (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), & !6 w tej linii
2668 chis(i,j),chis(j,i), & !2 w tej linii
2669 nstate(i,j),(wstate(k,i,j),k=1,4), & !5 w tej lini - 1 integer pierwszy
2670 dhead(1,1,i,j),dhead(1,2,i,j),dhead(2,1,i,j),dhead(2,2,i,j),& ! 4 w tej linii
2671 dtail(1,i,j),dtail(2,i,j), & ! 2 w tej lini
2672 epshead(i,j),sig0head(i,j), & ! 2 w tej linii
2673 rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j), & ! 5 w tej linii
2674 alphapol(i,j),alphapol(j,i), & ! 2 w tej linii
2675 (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j),debaykap(i,j) !8 w tej linii
2676 IF ((LaTeX).and.(i.gt.24)) then
2677 write (2,'(2a4,1h&,14(f8.2,1h&),23(f8.2,1h&))') restyp(i,1),restyp(j,1), &
2678 eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), & !6 w tej linii
2679 (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), & !6 w tej linii
2680 chis(i,j),chis(j,i) !2 w tej linii
2682 print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i), wqdip(1,i,j)
2687 IF ((LaTeX).and.(i.gt.24)) then
2688 write (2,'(2a4,1h&,14(f8.2,1h&),23(f8.2,1h&))') restyp(i,1),restyp(j,1), &
2689 dhead(1,1,i,j),dhead(2,1,i,j),& ! 2 w tej linii
2690 dtail(1,i,j),dtail(2,i,j), & ! 2 w tej lini
2691 epshead(i,j),sig0head(i,j), & ! 2 w tej linii
2692 rborn(i,j),rborn(j,i), & ! 3 w tej linii
2693 alphapol(i,j),alphapol(j,i), & ! 2 w tej linii
2694 (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j),debaykap(i,j) !8 w tej linii
2701 sigma(i,j) = sigma(j,i)
2702 sigmap1(i,j)=sigmap1(j,i)
2703 sigmap2(i,j)=sigmap2(j,i)
2704 sigiso1(i,j)=sigiso1(j,i)
2705 sigiso2(i,j)=sigiso2(j,i)
2706 ! print *,"ATU",sigma(j,i),sigma(i,j),i,j
2707 nstate(i,j) = nstate(j,i)
2708 dtail(1,i,j) = dtail(1,j,i)
2709 dtail(2,i,j) = dtail(2,j,i)
2711 alphasur(k,i,j) = alphasur(k,j,i)
2712 wstate(k,i,j) = wstate(k,j,i)
2713 alphiso(k,i,j) = alphiso(k,j,i)
2716 dhead(2,1,i,j) = dhead(1,1,j,i)
2717 dhead(2,2,i,j) = dhead(1,2,j,i)
2718 dhead(1,1,i,j) = dhead(2,1,j,i)
2719 dhead(1,2,i,j) = dhead(2,2,j,i)
2721 epshead(i,j) = epshead(j,i)
2722 sig0head(i,j) = sig0head(j,i)
2725 wqdip(k,i,j) = wqdip(k,j,i)
2728 wquad(i,j) = wquad(j,i)
2729 epsintab(i,j) = epsintab(j,i)
2730 debaykap(i,j)=debaykap(j,i)
2731 ! if (epsintab(i,j).ne.1.0) print *,"WHAT?",i,j,epsintab(i,j)
2738 !--------------------- GBV potential -----------------------------------
2740 ! 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2741 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2742 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
2743 (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
2745 write (iout,'(/a/)') 'Parameters of the GBV potential:'
2746 write (iout,'(a/)') 'The epsilon array:'
2747 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2748 write (iout,'(/a)') 'One-body parameters:'
2749 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',&
2750 's||/s_|_^2',' chip ',' alph '
2751 write (iout,'(a3,6x,5f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2752 sigii(i),chip(i),alp(i),i=1,ntyp)
2755 write(iout,*)"Wrong ipot"
2761 !-----------------------------------------------------------------------
2762 ! Calculate the "working" parameters of SC interactions.
2764 !el from module energy - COMMON.INTERACT-------
2765 allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1))
2766 if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2767 allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)
2768 if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2769 allocate(r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
2770 allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),&
2772 allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),&
2778 if (scelemode.eq.0) then
2787 sc_aa_tube_par(:)=0.0d0
2788 sc_bb_tube_par(:)=0.0d0
2790 !--------------------------------
2795 epslip(i,j)=epslip(j,i)
2798 if (scelemode.eq.0) then
2801 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
2802 sigma(j,i)=sigma(i,j)
2803 rs0(i,j)=dwa16*sigma(i,j)
2808 if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
2809 'Working parameters of the SC interactions:',&
2810 ' a ',' b ',' augm ',' sigma ',' r0 ',&
2815 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2817 ! print *,"SIGMA ZLA?",sigma(i,j)
2825 sigeps=dsign(1.0D0,epsij)
2827 aa_aq(i,j)=epsij*rrij*rrij
2828 ! print *,"ADASKO",epsij,rrij,expon
2829 bb_aq(i,j)=-sigeps*epsij*rrij
2830 aa_aq(j,i)=aa_aq(i,j)
2831 bb_aq(j,i)=bb_aq(i,j)
2832 epsijlip=epslip(i,j)
2833 sigeps=dsign(1.0D0,epsijlip)
2834 epsijlip=dabs(epsijlip)
2835 aa_lip(i,j)=epsijlip*rrij*rrij
2836 bb_lip(i,j)=-sigeps*epsijlip*rrij
2837 aa_lip(j,i)=aa_lip(i,j)
2838 bb_lip(j,i)=bb_lip(i,j)
2839 !C write(iout,*) aa_lip
2840 if ((ipot.gt.2).and. (scelemode.eq.0)) then
2841 sigt1sq=sigma0(i)**2
2842 sigt2sq=sigma0(j)**2
2845 ratsig1=sigt2sq/sigt1sq
2846 ratsig2=1.0D0/ratsig1
2847 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
2848 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
2849 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
2853 ! if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2854 sigmaii(i,j)=rsum_max
2855 sigmaii(j,i)=rsum_max
2857 ! sigmaii(i,j)=r0(i,j)
2858 ! sigmaii(j,i)=r0(i,j)
2860 !d write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
2861 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
2862 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
2863 augm(i,j)=epsij*r_augm**(2*expon)
2864 ! augm(i,j)=0.5D0**(2*expon)*aa(i,j)
2871 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
2872 restyp(i,1),restyp(j,1),aa_aq(i,j),bb_aq(i,j),augm(i,j),&
2873 sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
2878 allocate(eps_nucl(ntyp_molec(2),ntyp_molec(2)))
2879 allocate(sigma_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2880 allocate(elpp6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2881 allocate(elpp3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2882 allocate(elpp63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2883 allocate(elpp32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2884 allocate(chi_nucl(ntyp_molec(2),ntyp_molec(2)),chip_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2885 allocate(ael3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2886 allocate(ael6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2887 allocate(ael32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2888 allocate(ael63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2889 allocate(aa_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2890 allocate(bb_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2891 allocate(r0_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2892 allocate(sigmaii_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2893 allocate(eps_scp_nucl(ntyp_molec(2)),rscp_nucl(ntyp_molec(2))) !(ntyp,2)
2902 read (isidep_nucl,*) ipot_nucl
2903 ! print *,"TU?!",ipot_nucl
2904 if (ipot_nucl.eq.1) then
2905 do i=1,ntyp_molec(2)
2906 do j=i,ntyp_molec(2)
2907 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),elpp6_nucl(i,j),&
2908 elpp3_nucl(i,j), elpp63_nucl(i,j),elpp32_nucl(i,j)
2912 do i=1,ntyp_molec(2)
2913 do j=i,ntyp_molec(2)
2914 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),chi_nucl(i,j),&
2915 chi_nucl(j,i),chip_nucl(i,j),chip_nucl(j,i),&
2916 elpp6_nucl(i,j),elpp3_nucl(i,j),elpp63_nucl(i,j),elpp32_nucl(i,j)
2920 ! rpp(1,1)=2**(1.0/6.0)*5.16158
2921 do i=1,ntyp_molec(2)
2922 do j=i,ntyp_molec(2)
2923 rrij=sigma_nucl(i,j)
2927 epsij=4*eps_nucl(i,j)
2928 sigeps=dsign(1.0D0,epsij)
2930 aa_nucl(i,j)=epsij*rrij*rrij
2931 bb_nucl(i,j)=-sigeps*epsij*rrij
2932 ael3_nucl(i,j)=elpp3_nucl(i,j)*dsqrt(rrij)
2933 ael6_nucl(i,j)=elpp6_nucl(i,j)*rrij
2934 ael63_nucl(i,j)=elpp63_nucl(i,j)*rrij
2935 ael32_nucl(i,j)=elpp32_nucl(i,j)*rrij
2936 sigmaii_nucl(i,j)=sigma_nucl(i,j)/sqrt(1-(chi_nucl(i,j)+chi_nucl(j,i)- &
2937 2*chi_nucl(i,j)*chi_nucl(j,i))/(1-chi_nucl(i,j)*chi_nucl(j,i)))
2940 aa_nucl(i,j)=aa_nucl(j,i)
2941 bb_nucl(i,j)=bb_nucl(j,i)
2942 ael3_nucl(i,j)=ael3_nucl(j,i)
2943 ael6_nucl(i,j)=ael6_nucl(j,i)
2944 ael63_nucl(i,j)=ael63_nucl(j,i)
2945 ael32_nucl(i,j)=ael32_nucl(j,i)
2946 elpp3_nucl(i,j)=elpp3_nucl(j,i)
2947 elpp6_nucl(i,j)=elpp6_nucl(j,i)
2948 elpp63_nucl(i,j)=elpp63_nucl(j,i)
2949 elpp32_nucl(i,j)=elpp32_nucl(j,i)
2950 eps_nucl(i,j)=eps_nucl(j,i)
2951 sigma_nucl(i,j)=sigma_nucl(j,i)
2952 sigmaii_nucl(i,j)=sigmaii_nucl(j,i)
2956 write(iout,*) "tube param"
2957 read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep, &
2958 ccavtubpep,dcavtubpep,tubetranenepep
2959 sigmapeptube=sigmapeptube**6
2960 sigeps=dsign(1.0D0,epspeptube)
2961 epspeptube=dabs(epspeptube)
2962 pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
2963 pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
2964 write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
2966 read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i), &
2967 ccavtub(i),dcavtub(i),tubetranene(i)
2968 sigmasctube=sigmasctube**6
2969 sigeps=dsign(1.0D0,epssctube)
2970 epssctube=dabs(epssctube)
2971 sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
2972 sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
2973 write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
2975 !-----------------READING SC BASE POTENTIALS-----------------------------
2976 allocate(eps_scbase(ntyp_molec(1),ntyp_molec(2)))
2977 allocate(sigma_scbase(ntyp_molec(1),ntyp_molec(2)))
2978 allocate(chi_scbase(ntyp_molec(1),ntyp_molec(2),2))
2979 allocate(chipp_scbase(ntyp_molec(1),ntyp_molec(2),2))
2980 allocate(alphasur_scbase(4,ntyp_molec(1),ntyp_molec(2)))
2981 allocate(sigmap1_scbase(ntyp_molec(1),ntyp_molec(2)))
2982 allocate(sigmap2_scbase(ntyp_molec(1),ntyp_molec(2)))
2983 allocate(chis_scbase(ntyp_molec(1),ntyp_molec(2),2))
2984 allocate(dhead_scbasei(ntyp_molec(1),ntyp_molec(2)))
2985 allocate(dhead_scbasej(ntyp_molec(1),ntyp_molec(2)))
2986 allocate(rborn_scbasei(ntyp_molec(1),ntyp_molec(2)))
2987 allocate(rborn_scbasej(ntyp_molec(1),ntyp_molec(2)))
2988 allocate(wdipdip_scbase(3,ntyp_molec(1),ntyp_molec(2)))
2989 allocate(wqdip_scbase(2,ntyp_molec(1),ntyp_molec(2)))
2990 allocate(alphapol_scbase(ntyp_molec(1),ntyp_molec(2)))
2991 allocate(epsintab_scbase(ntyp_molec(1),ntyp_molec(2)))
2994 do i=1,ntyp_molec(1)
2995 do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
2996 write (*,*) "Im in ", i, " ", j
2997 read(isidep_scbase,*) &
2998 eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
2999 chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
3000 write(*,*) "eps",eps_scbase(i,j)
3001 read(isidep_scbase,*) &
3002 (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
3003 chis_scbase(i,j,1),chis_scbase(i,j,2)
3004 read(isidep_scbase,*) &
3005 dhead_scbasei(i,j), &
3006 dhead_scbasej(i,j), &
3007 rborn_scbasei(i,j),rborn_scbasej(i,j)
3008 read(isidep_scbase,*) &
3009 (wdipdip_scbase(k,i,j),k=1,3), &
3010 (wqdip_scbase(k,i,j),k=1,2)
3011 read(isidep_scbase,*) &
3012 alphapol_scbase(i,j), &
3013 epsintab_scbase(i,j)
3016 allocate(aa_scbase(ntyp_molec(1),ntyp_molec(2)))
3017 allocate(bb_scbase(ntyp_molec(1),ntyp_molec(2)))
3019 do i=1,ntyp_molec(1)
3020 do j=1,ntyp_molec(2)-1
3021 epsij=eps_scbase(i,j)
3022 rrij=sigma_scbase(i,j)
3027 sigeps=dsign(1.0D0,epsij)
3029 aa_scbase(i,j)=epsij*rrij*rrij
3030 bb_scbase(i,j)=-sigeps*epsij*rrij
3033 !-----------------READING PEP BASE POTENTIALS-------------------
3034 allocate(eps_pepbase(ntyp_molec(2)))
3035 allocate(sigma_pepbase(ntyp_molec(2)))
3036 allocate(chi_pepbase(ntyp_molec(2),2))
3037 allocate(chipp_pepbase(ntyp_molec(2),2))
3038 allocate(alphasur_pepbase(4,ntyp_molec(2)))
3039 allocate(sigmap1_pepbase(ntyp_molec(2)))
3040 allocate(sigmap2_pepbase(ntyp_molec(2)))
3041 allocate(chis_pepbase(ntyp_molec(2),2))
3042 allocate(wdipdip_pepbase(3,ntyp_molec(2)))
3045 do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
3046 write (*,*) "Im in ", i, " ", j
3047 read(isidep_pepbase,*) &
3048 eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
3049 chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
3050 write(*,*) "eps",eps_pepbase(j)
3051 read(isidep_pepbase,*) &
3052 (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
3053 chis_pepbase(j,1),chis_pepbase(j,2)
3054 read(isidep_pepbase,*) &
3055 (wdipdip_pepbase(k,j),k=1,3)
3057 allocate(aa_pepbase(ntyp_molec(2)))
3058 allocate(bb_pepbase(ntyp_molec(2)))
3060 do j=1,ntyp_molec(2)-1
3061 epsij=eps_pepbase(j)
3062 rrij=sigma_pepbase(j)
3067 sigeps=dsign(1.0D0,epsij)
3069 aa_pepbase(j)=epsij*rrij*rrij
3070 bb_pepbase(j)=-sigeps*epsij*rrij
3072 !--------------READING SC PHOSPHATE-------------------------------------
3073 allocate(eps_scpho(ntyp_molec(1)))
3074 allocate(sigma_scpho(ntyp_molec(1)))
3075 allocate(chi_scpho(ntyp_molec(1),2))
3076 allocate(chipp_scpho(ntyp_molec(1),2))
3077 allocate(alphasur_scpho(4,ntyp_molec(1)))
3078 allocate(sigmap1_scpho(ntyp_molec(1)))
3079 allocate(sigmap2_scpho(ntyp_molec(1)))
3080 allocate(chis_scpho(ntyp_molec(1),2))
3081 allocate(wqq_scpho(ntyp_molec(1)))
3082 allocate(wqdip_scpho(2,ntyp_molec(1)))
3083 allocate(alphapol_scpho(ntyp_molec(1)))
3084 allocate(epsintab_scpho(ntyp_molec(1)))
3085 allocate(dhead_scphoi(ntyp_molec(1)))
3086 allocate(rborn_scphoi(ntyp_molec(1)))
3087 allocate(rborn_scphoj(ntyp_molec(1)))
3088 allocate(alphi_scpho(ntyp_molec(1)))
3092 do j=1,ntyp_molec(1) ! without U then we will take T for U
3093 write (*,*) "Im in scpho ", i, " ", j
3094 read(isidep_scpho,*) &
3095 eps_scpho(j),sigma_scpho(j),chi_scpho(j,1),&
3096 chi_scpho(j,2),chipp_scpho(j,1),chipp_scpho(j,2)
3097 write(*,*) "eps",eps_scpho(j)
3098 read(isidep_scpho,*) &
3099 (alphasur_scpho(k,j),k=1,4),sigmap1_scpho(j),sigmap2_scpho(j), &
3100 chis_scpho(j,1),chis_scpho(j,2)
3101 read(isidep_scpho,*) &
3102 (wqdip_scpho(k,j),k=1,2),wqq_scpho(j),dhead_scphoi(j)
3103 read(isidep_scpho,*) &
3104 epsintab_scpho(j),alphapol_scpho(j),rborn_scphoi(j),rborn_scphoj(j), &
3108 allocate(aa_scpho(ntyp_molec(1)))
3109 allocate(bb_scpho(ntyp_molec(1)))
3111 do j=1,ntyp_molec(1)
3118 sigeps=dsign(1.0D0,epsij)
3120 aa_scpho(j)=epsij*rrij*rrij
3121 bb_scpho(j)=-sigeps*epsij*rrij
3125 read(isidep_peppho,*) &
3126 eps_peppho,sigma_peppho
3127 read(isidep_peppho,*) &
3128 (alphasur_peppho(k),k=1,4),sigmap1_peppho,sigmap2_peppho
3129 read(isidep_peppho,*) &
3130 (wqdip_peppho(k),k=1,2)
3138 sigeps=dsign(1.0D0,epsij)
3140 aa_peppho=epsij*rrij*rrij
3141 bb_peppho=-sigeps*epsij*rrij
3144 allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
3149 ! Define the SC-p interaction constants (hard-coded; old style)
3152 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
3154 ! aad(i,1)=0.3D0*4.0D0**12
3155 ! Following line for constants currently implemented
3156 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
3157 aad(i,1)=1.5D0*4.0D0**12
3158 ! aad(i,1)=0.17D0*5.6D0**12
3160 ! "Soft" SC-p repulsion
3162 ! Following line for constants currently implemented
3163 ! aad(i,1)=0.3D0*4.0D0**6
3164 ! "Hard" SC-p repulsion
3165 bad(i,1)=3.0D0*4.0D0**6
3166 ! bad(i,1)=-2.0D0*0.17D0*5.6D0**6
3175 ! 8/9/01 Read the SC-p interaction constants from file
3178 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
3181 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
3182 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
3183 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
3184 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
3188 write (iout,*) "Parameters of SC-p interactions:"
3190 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
3191 eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
3196 allocate(aad_nucl(ntyp_molec(2)),bad_nucl(ntyp_molec(2))) !(ntyp,2)
3198 do i=1,ntyp_molec(2)
3199 read (iscpp_nucl,*,end=118,err=118) eps_scp_nucl(i),rscp_nucl(i)
3201 do i=1,ntyp_molec(2)
3202 aad_nucl(i)=dabs(eps_scp_nucl(i))*rscp_nucl(i)**12
3203 bad_nucl(i)=-2*eps_scp_nucl(i)*rscp_nucl(i)**6
3205 r0pp=1.12246204830937298142*5.16158
3211 ! Define the constants of the disulfide bridge
3215 ! Old arbitrary potential - commented out.
3220 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
3221 ! energy surface of diethyl disulfide.
3222 ! A. Liwo and U. Kozlowska, 11/24/03
3240 allocate(alphapolcat(ntyp,ntyp),epsheadcat(ntyp,ntyp),sig0headcat(ntyp,ntyp))
3241 allocate(sigiso1cat(ntyp,ntyp),rborn1cat(ntyp,ntyp),rborn2cat(ntyp,ntyp),sigmap1cat(ntyp,ntyp))
3242 allocate(sigmap2cat(ntyp,ntyp),sigiso2cat(ntyp,ntyp))
3243 allocate(chis1cat(ntyp,ntyp),chis2cat(ntyp,ntyp),wquadcat(ntyp,ntyp),chipp1cat(ntyp,ntyp),chipp2cat(ntyp,ntyp))
3244 allocate(epsintabcat(ntyp,ntyp))
3245 allocate(dtailcat(2,ntyp,ntyp))
3246 allocate(alphasurcat(4,ntyp,ntyp),alphisocat(4,ntyp,ntyp))
3247 allocate(wqdipcat(2,ntyp,ntyp))
3248 allocate(wstatecat(4,ntyp,ntyp))
3249 allocate(dheadcat(2,2,ntyp,ntyp))
3250 allocate(nstatecat(ntyp,ntyp))
3251 allocate(debaykapcat(ntyp,ntyp))
3253 if (.not.allocated(epscat)) allocate (epscat(0:ntyp1,0:ntyp1))
3254 if (.not.allocated(sigmacat)) allocate(sigmacat(0:ntyp1,0:ntyp1))
3255 ! if (.not.allocated(chicat)) allocate(chicat(ntyp1,ntyp1)) !(ntyp,ntyp)
3256 if (.not.allocated(chi1cat)) allocate(chi1cat(ntyp1,ntyp1)) !(ntyp,ntyp)
3257 if (.not.allocated(chi2cat)) allocate(chi2cat(ntyp1,ntyp1)) !(ntyp,ntyp)
3260 allocate (ichargecat(ntyp_molec(5)))
3261 ! i to SC, j to jon, isideocat - nazwa pliku z ktorego czytam parametry
3262 if (oldion.eq.0) then
3263 if (.not.allocated(icharge)) then ! this mean you are oprating in old sc-sc mode
3264 allocate(icharge(1:ntyp1))
3265 read(iion,*) (icharge(i),i=1,ntyp)
3270 do i=1,ntyp_molec(5)
3271 read(iion,*) msc(i,5),restok(i,5),ichargecat(i)
3272 print *,msc(i,5),restok(i,5)
3276 do j=1,ntyp_molec(5)
3278 ! do j=1,ntyp_molec(5)
3279 ! write (*,*) "Im in ALAB", i, " ", j
3281 epscat(i,j),sigmacat(i,j), &
3282 ! chicat(i,j),chicat(j,i),chippcat(i,j),chippcat(j,i), &
3283 chi1cat(i,j),chi2cat(i,j),chipp1cat(i,j),chipp2cat(i,j), &
3285 (alphasurcat(k,i,j),k=1,4),sigmap1cat(i,j),sigmap2cat(i,j),&
3286 ! chiscat(i,j),chiscat(j,i), &
3287 chis1cat(i,j),chis2cat(i,j), &
3289 nstatecat(i,j),(wstatecat(k,i,j),k=1,4), & !5 w tej lini - 1 integer pierwszy
3290 dheadcat(1,1,i,j),dheadcat(1,2,i,j),dheadcat(2,1,i,j),dheadcat(2,2,i,j),&
3291 dtailcat(1,i,j),dtailcat(2,i,j), &
3292 epsheadcat(i,j),sig0headcat(i,j), &
3294 ! rborncat(i,j),rborncat(j,i),&
3295 rborn1cat(i,j),rborn2cat(i,j),&
3296 (wqdipcat(k,i,j),k=1,2), &
3297 alphapolcat(i,j),alphapolcat(j,i), &
3298 (alphisocat(k,i,j),k=1,4),sigiso1cat(i,j),sigiso2cat(i,j),epsintabcat(i,j),debaykapcat(i,j)
3299 ! print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i)
3301 ! write (iout,*) 'i= ', i, ' j= ', j
3302 ! write (iout,*) 'epsi0= ', epscat(1,j)
3303 ! write (iout,*) 'sigma0= ', sigmacat(1,j)
3304 ! write (iout,*) 'chi(i,j)= ', chicat(1,j)
3305 ! write (iout,*) 'chi(j,i)= ', chicat(j,1)
3306 ! write (iout,*) 'chip(i,j)= ', chippcat(1,j)
3307 ! write (iout,*) 'chip(j,i)= ', chippcat(j,1)
3308 ! write (iout,*) 'alphasur1= ', alphasurcat(1,1,j)
3309 ! write (iout,*) 'alphasur2= ', alphasurcat(2,1,j)
3310 ! write (iout,*) 'alphasur3= ', alphasurcat(3,1,j)
3311 ! write (iout,*) 'alphasur4= ', alphasurcat(4,1,j)
3312 ! write (iout,*) 'sig1= ', sigmap1cat(1,j)
3313 ! write (iout,*) 'chis(i,j)= ', chiscat(1,j)
3314 ! write (iout,*) 'chis(j,i)= ', chiscat(j,1)
3315 ! write (iout,*) 'dhead= ', dheadcat(1,1,1,j)
3316 ! write (iout,*) 'a1= ', rborncat(j,1)
3317 ! write (iout,*) 'a2= ', rborncat(1,j)
3318 ! write (iout,*) 'epsin= ', epsintabcat(1,j), epsintabcat(j,1)
3319 ! write (iout,*) 'alphapol1= ', alphapolcat(1,j)
3320 ! write (iout,*) 'w1= ', wqdipcat(1,1,j)
3321 ! write (iout,*) 'w2= ', wqdipcat(2,1,j)
3325 ! If ((i.eq.1).and.(j.eq.27)) then
3326 ! write (iout,*) 'SEP'
3327 ! Write (iout,*) 'w1= ', wqdipcat(1,1,27)
3328 ! Write (iout,*) 'w2= ', wqdipcat(2,1,27)
3333 allocate(aa_aq_cat(-ntyp:ntyp,ntyp),bb_aq_cat(-ntyp:ntyp,ntyp))
3335 do j=1,ntyp_molec(5)
3339 sigeps=dsign(1.0D0,epsij)
3341 aa_aq_cat(i,j)=epsij*rrij*rrij
3342 bb_aq_cat(i,j)=-sigeps*epsij*rrij
3347 do j=1,ntyp_molec(5)
3349 write (iout,*) 'i= ', i, ' j= ', j
3350 write (iout,*) 'epsi0= ', epscat(i,j)
3351 write (iout,*) 'sigma0= ', sigmacat(i,j)
3352 write (iout,*) 'chi1= ', chi1cat(i,j)
3353 write (iout,*) 'chi1= ', chi2cat(i,j)
3354 write (iout,*) 'chip1= ', chipp1cat(1,j)
3355 write (iout,*) 'chip2= ', chipp2cat(1,j)
3356 write (iout,*) 'alphasur1= ', alphasurcat(1,1,j)
3357 write (iout,*) 'alphasur2= ', alphasurcat(2,1,j)
3358 write (iout,*) 'alphasur3= ', alphasurcat(3,1,j)
3359 write (iout,*) 'alphasur4= ', alphasurcat(4,1,j)
3360 write (iout,*) 'sig1= ', sigmap1cat(1,j)
3361 write (iout,*) 'sig2= ', sigmap2cat(1,j)
3362 write (iout,*) 'chis1= ', chis1cat(1,j)
3363 write (iout,*) 'chis1= ', chis2cat(1,j)
3364 write (iout,*) 'nstatecat(i,j)= ', nstatecat(1,j)
3365 write (iout,*) 'wstatecat(k,i,j),k=1= ',wstatecat(1,1,j)
3366 write (iout,*) 'dhead= ', dheadcat(1,1,1,j)
3367 write (iout,*) 'dhead2= ', dheadcat(1,2,1,j)
3368 write (iout,*) 'a1= ', rborn1cat(i,j)
3369 write (iout,*) 'a2= ', rborn2cat(i,j)
3370 write (iout,*) 'epsin= ', epsintabcat(1,j), epsintabcat(j,1)
3371 write (iout,*) 'alphapol1= ', alphapolcat(1,j)
3372 write (iout,*) 'alphapol2= ', alphapolcat(j,1)
3373 write (iout,*) 'w1= ', wqdipcat(1,i,j)
3374 write (iout,*) 'w2= ', wqdipcat(2,i,j)
3375 write (iout,*) 'debaykapcat(i,j)= ', debaykapcat(1,j)
3378 If ((i.eq.1).and.(j.eq.27)) then
3379 write (iout,*) 'SEP'
3380 Write (iout,*) 'w1= ', wqdipcat(1,1,27)
3381 Write (iout,*) 'w2= ', wqdipcat(2,1,27)
3391 write (iout,'(/a)') "Disulfide bridge parameters:"
3392 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
3393 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
3394 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
3395 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
3398 if (shield_mode.gt.0) then
3399 pi=4.0D0*datan(1.0D0)
3400 !C VSolvSphere the volume of solving sphere
3402 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
3403 !C there will be no distinction between proline peptide group and normal peptide
3404 !C group in case of shielding parameters
3405 VSolvSphere=4.0/3.0*pi*(4.50d0)**3
3406 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3
3407 write (iout,*) VSolvSphere,VSolvSphere_div
3408 !C long axis of side chain
3410 long_r_sidechain(i)=vbldsc0(1,i)
3411 ! if (scelemode.eq.0) then
3412 short_r_sidechain(i)=sigma(i,i)/sqrt(2.0)
3413 if (short_r_sidechain(i).eq.0.0) short_r_sidechain(i)=0.2
3415 ! short_r_sidechain(i)=sigma(i,i)
3417 write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
3424 111 write (iout,*) "Error reading bending energy parameters."
3426 112 write (iout,*) "Error reading rotamer energy parameters."
3428 113 write (iout,*) "Error reading torsional energy parameters."
3430 114 write (iout,*) "Error reading double torsional energy parameters."
3432 115 write (iout,*) &
3433 "Error reading cumulant (multibody energy) parameters."
3435 116 write (iout,*) "Error reading electrostatic energy parameters."
3437 117 write (iout,*) "Error reading side chain interaction parameters."
3439 118 write (iout,*) "Error reading SCp interaction parameters."
3441 119 write (iout,*) "Error reading SCCOR parameters"
3443 121 write (iout,*) "Error in Czybyshev parameters"
3446 call MPI_Finalize(Ierror)
3450 end subroutine parmread
3452 !-----------------------------------------------------------------------------
3454 !-----------------------------------------------------------------------------
3455 subroutine printmat(ldim,m,n,iout,key,a)
3458 character(len=3),dimension(n) :: key
3459 real(kind=8),dimension(ldim,n) :: a
3461 integer :: i,j,k,m,iout,nlim
3465 write (iout,1000) (key(k),k=i,nlim)
3467 1000 format (/5x,8(6x,a3))
3468 1020 format (/80(1h-)/)
3470 write (iout,1010) key(j),(a(j,k),k=i,nlim)
3473 1010 format (a3,2x,8(f9.4))
3475 end subroutine printmat
3476 !-----------------------------------------------------------------------------
3478 !-----------------------------------------------------------------------------
3480 ! Read the PDB file and convert the peptide geometry into virtual-chain
3483 use energy_data, only: itype,istype
3487 ! use control, only: rescode,sugarcode
3488 ! implicit real*8 (a-h,o-z)
3489 ! include 'DIMENSIONS'
3490 ! include 'COMMON.LOCAL'
3491 ! include 'COMMON.VAR'
3492 ! include 'COMMON.CHAIN'
3493 ! include 'COMMON.INTERACT'
3494 ! include 'COMMON.IOUNITS'
3495 ! include 'COMMON.GEO'
3496 ! include 'COMMON.NAMES'
3497 ! include 'COMMON.CONTROL'
3498 ! include 'COMMON.DISTFIT'
3499 ! include 'COMMON.SETUP'
3500 integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift,k!,ity!,&
3502 logical :: lprn=.true.,fail
3503 real(kind=8),dimension(3) :: e1,e2,e3
3504 real(kind=8) :: dcj,efree_temp
3505 character(len=3) :: seq,res,res2
3506 character(len=5) :: atom
3507 character(len=80) :: card
3508 real(kind=8),dimension(3,20) :: sccor
3509 integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin !rescode,
3510 integer :: isugar,molecprev,firstion
3511 character*1 :: sugar
3513 real(kind=8),dimension(3) :: ccc
3515 integer,dimension(2,maxres/3) :: hfrag_alloc
3516 integer,dimension(4,maxres/3) :: bfrag_alloc
3517 real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm)
3518 real(kind=8),dimension(:,:), allocatable :: c_temporary
3519 integer,dimension(:,:) , allocatable :: itype_temporary
3520 integer,dimension(:),allocatable :: istype_temp
3527 ! write (2,*) "UNRES_PDB",unres_pdb
3547 !-----------------------------
3548 allocate(hfrag(2,maxres/3)) !(2,maxres/3)
3549 allocate(bfrag(4,maxres/3)) !(4,maxres/3)
3550 if(.not. allocated(istype)) allocate(istype(maxres))
3552 read (ipdbin,'(a80)',end=10) card
3553 write (iout,'(a)') card
3554 if (card(:5).eq.'HELIX') then
3557 read(card(22:25),*) hfrag(1,nhfrag)
3558 read(card(34:37),*) hfrag(2,nhfrag)
3560 if (card(:5).eq.'SHEET') then
3563 read(card(24:26),*) bfrag(1,nbfrag)
3564 read(card(35:37),*) bfrag(2,nbfrag)
3565 !rc----------------------------------------
3566 !rc to be corrected !!!
3567 bfrag(3,nbfrag)=bfrag(1,nbfrag)
3568 bfrag(4,nbfrag)=bfrag(2,nbfrag)
3569 !rc----------------------------------------
3571 if (card(:3).eq.'END') then
3573 else if (card(:3).eq.'TER') then
3578 itype(ires_old,molecule)=ntyp1_molec(molecule)
3579 itype(ires_old-1,molecule)=ntyp1_molec(molecule)
3580 nres_molec(molecule)=nres_molec(molecule)+2
3582 ! write (iout,*) "Chain ended",ires,ishift,ires_old
3585 dc(j,ires)=sccor(j,iii)
3588 call sccenter(ires,iii,sccor)
3594 if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
3595 ! Fish out the ATOM cards.
3596 ! write(iout,*) 'card',card(1:20)
3597 ! print *,"ATU ",card(1:6), CARD(3:6)
3599 if (index(card(1:4),'ATOM').gt.0) then
3600 read (card(12:16),*) atom
3601 ! write (iout,*) "! ",atom," !",ires
3602 ! if (atom.eq.'CA' .or. atom.eq.'CH3') then
3603 read (card(23:26),*) ires
3604 read (card(18:20),'(a3)') res
3605 ! write (iout,*) "ires",ires,ires-ishift+ishift1,
3606 ! & " ires_old",ires_old
3607 ! write (iout,*) "ishift",ishift," ishift1",ishift1
3608 ! write (iout,*) "IRES",ires-ishift+ishift1,ires_old
3609 if (ires-ishift+ishift1.ne.ires_old) then
3610 ! Calculate the CM of the preceding residue.
3611 ! if (ibeg.eq.0) call sccenter(ires,iii,sccor)
3613 ! write (iout,*) "Calculating sidechain center iii",iii
3616 dc(j,ires+ishift1-ishift-1)=sccor(j,iii)
3619 call sccenter(ires_old,iii,sccor)
3623 ! Start new residue.
3624 if (res.eq.'Cl-' .or. res.eq.'Na+') then
3627 else if (ibeg.eq.1) then
3628 write (iout,*) "BEG ires",ires
3630 if (res.ne.'GLY' .and. res.ne. 'ACE') then
3633 nres_molec(molecule)=nres_molec(molecule)+1
3635 ires=ires-ishift+ishift1
3637 ! write (iout,*) "ishift",ishift," ires",ires,&
3638 ! " ires_old",ires_old
3640 else if (ibeg.eq.2) then
3642 ishift=-ires_old+ires-1 !!!!!
3643 ishift1=ishift1-1 !!!!!
3644 ! write (iout,*) "New chain started",ires,ishift,ishift1,"!"
3645 ires=ires-ishift+ishift1
3646 ! print *,ires,ishift,ishift1
3650 ishift=ishift-(ires-ishift+ishift1-ires_old-1)
3651 ires=ires-ishift+ishift1
3654 ! print *,'atom',ires,atom
3655 if (res.eq.'ACE' .or. res.eq.'NHE') then
3658 if (atom.eq.'CA '.or.atom.eq.'N ') then
3660 itype(ires,molecule)=rescode(ires,res,0,molecule)
3662 ! nres_molec(molecule)=nres_molec(molecule)+1
3666 itype(ires,molecule)=rescode(ires,res2,0,molecule)
3667 ! nres_molec(molecule)=nres_molec(molecule)+1
3668 read (card(19:19),'(a1)') sugar
3669 isugar=sugarcode(sugar,ires)
3670 ! if (ibeg.eq.1) then
3674 ! print *,"ires=",ires,istype(ires)
3680 ires=ires-ishift+ishift1
3682 ! write (iout,*) "ires_old",ires_old," ires",ires
3683 if (card(27:27).eq."A" .or. card(27:27).eq."B") then
3686 ! write (2,*) "ires",ires," res ",res!," ity"!,ity
3687 if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
3688 res.eq.'NHE'.and.atom(:2).eq.'HN') then
3689 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3690 ! print *,ires,ishift,ishift1
3691 ! write (iout,*) "backbone ",atom
3693 write (iout,'(2i3,2x,a,3f8.3)') &
3694 ires,itype(ires,1),res,(c(j,ires),j=1,3)
3697 nres_molec(molecule)=nres_molec(molecule)+1
3699 sccor(j,iii)=c(j,ires)
3701 else if (.not.unres_pdb .and. (atom.eq."C1'" .or. &
3702 atom.eq."C2'" .or. atom.eq."C3'" &
3703 .or. atom.eq."C4'" .or. atom.eq."O4'")) then
3704 read(card(31:54),'(3f8.3)') (ccc(j),j=1,3)
3705 !c write (2,'(i5,3f10.5)') ires,(ccc(j),j=1,3)
3706 ! print *,ires,ishift,ishift1
3710 ! sccor(j,iii)=c(j,ires)
3713 c(j,ires)=c(j,ires)+ccc(j)/5.0
3715 print *,counter,molecule
3716 if (counter.eq.5) then
3718 nres_molec(molecule)=nres_molec(molecule)+1
3721 ! sccor(j,iii)=c(j,ires)
3725 ! print *, "ATOM",atom(1:3)
3726 else if (atom.eq."C5'") then
3727 read (card(19:19),'(a1)') sugar
3728 isugar=sugarcode(sugar,ires)
3733 ! print *,ires,istype(ires)
3737 ! print *,"nres_molec(molecule)",nres_molec(molecule),ires
3738 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3739 nres_molec(molecule)=nres_molec(molecule)+1
3740 print *,"nres_molec(molecule)",nres_molec(molecule),ires
3744 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
3746 else if ((atom.eq."C1'").and.unres_pdb) then
3748 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
3749 ! write (*,*) card(23:27),ires,itype(ires,1)
3750 else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. &
3751 atom.ne.'N' .and. atom.ne.'C' .and. &
3752 atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. &
3753 atom.ne.'OXT' .and. atom(:2).ne.'3H' &
3754 .and. atom.ne.'P '.and. &
3755 atom(1:1).ne.'H' .and. &
3756 atom.ne.'OP1' .and. atom.ne.'OP2 '.and. atom.ne.'OP3'&
3758 ! write (iout,*) "sidechain ",atom
3759 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
3760 if ((molecule.ne.2).or.(atom(3:3).ne."'")) then
3761 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
3763 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
3766 ! print *,"IONS",ions,card(1:6)
3767 else if ((ions).and.(card(1:6).eq.'HETATM')) then
3768 if (firstion.eq.0) then
3772 dc(j,ires)=sccor(j,iii)
3775 call sccenter(ires,iii,sccor)
3778 read (card(12:16),*) atom
3779 ! print *,"HETATOM", atom
3780 read (card(18:20),'(a3)') res
3781 if ((atom(1:2).eq.'NA').or.(atom(1:2).eq.'CL').or.&
3782 (atom(1:2).eq.'CA').or.(atom(1:2).eq.'MG') &
3783 .or.(atom(1:2).eq.'K ')) &
3786 if (molecule.ne.5) molecprev=molecule
3788 nres_molec(molecule)=nres_molec(molecule)+1
3789 print *,"HERE",nres_molec(molecule)
3791 itype(ires,molecule)=rescode(ires,res,0,molecule)
3792 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3796 10 write (iout,'(a,i5)') ' Number of residues found: ',ires
3797 if (ires.eq.0) return
3798 ! Calculate dummy residue coordinates inside the "chain" of a multichain
3801 if (((ires_old.ne.ires).and.(molecule.ne.5)) &
3803 nres_molec(molecule)=nres_molec(molecule)-2
3804 print *,'I have',nres, nres_molec(:)
3806 do k=1,4 ! ions are without dummy
3807 if (nres_molec(k).eq.0) cycle
3809 ! write (iout,*) i,itype(i,1)
3810 ! if (itype(i,1).eq.ntyp1) then
3811 ! write (iout,*) "dummy",i,itype(i,1)
3813 ! c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
3814 ! c(j,i)=(c(j,i-1)+c(j,i+1))/2
3818 if (itype(i,k).eq.ntyp1_molec(k)) then
3819 if (itype(i+1,k).eq.ntyp1_molec(k)) then
3820 if (itype(i+2,k).eq.0) then
3821 ! print *,"masz sieczke"
3823 if (itype(i+2,j).ne.0) then
3825 itype(i+1,j)=ntyp1_molec(j)
3826 nres_molec(k)=nres_molec(k)-1
3827 nres_molec(j)=nres_molec(j)+1
3833 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
3834 ! first is connected prevous chain (itype(i+1,1).eq.ntyp1)=true
3835 ! second dummy atom is conected to next chain itype(i+1,1).eq.ntyp1=false
3836 ! if (unres_pdb) then
3837 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
3838 ! print *,i,'tu dochodze'
3839 ! call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
3847 ! c(j,i)=c(j,i-1)-1.9d0*e2(j)
3851 dcj=(c(j,i-2)-c(j,i-3))/2.0
3852 if (dcj.eq.0) dcj=1.23591524223
3857 else !itype(i+1,1).eq.ntyp1
3858 ! if (unres_pdb) then
3859 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3860 ! call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
3867 ! c(j,i)=c(j,i+1)-1.9d0*e2(j)
3868 c(j,i)=c(j,i-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0)
3872 dcj=(c(j,i+3)-c(j,i+2))/2.0
3873 if (dcj.eq.0) dcj=1.23591524223
3878 endif !itype(i+1,1).eq.ntyp1
3879 endif !itype.eq.ntyp1
3883 ! Calculate the CM of the last side chain.
3887 dc(j,ires)=sccor(j,iii)
3890 call sccenter(ires,iii,sccor)
3896 ! print *,"molecule",molecule
3897 if ((itype(nres,1).ne.10)) then
3899 if (molecule.eq.5) molecule=molecprev
3900 itype(nres,molecule)=ntyp1_molec(molecule)
3901 nres_molec(molecule)=nres_molec(molecule)+1
3902 ! if (unres_pdb) then
3903 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
3904 ! call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
3911 ! c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
3915 dcj=(c(j,nres-2)-c(j,nres-3))/2.0
3916 c(j,nres)=c(j,nres-1)+dcj
3917 c(j,2*nres)=c(j,nres)
3921 ! print *,'I have',nres, nres_molec(:)
3923 !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
3925 if (nres.ne.nres0) then
3926 write (iout,*) "Error: wrong parameter value: NRES=",nres,&
3928 stop "Error nres value in WHAM input"
3931 !---------------------------------
3932 !el reallocate tables
3935 ! hfrag_alloc(j,i)=hfrag(j,i)
3938 ! bfrag_alloc(j,i)=bfrag(j,i)
3944 ! allocate(hfrag(2,nres/3)) !(2,maxres/3)
3945 !el allocate(hfrag(2,nhfrag)) !(2,maxres/3)
3946 !el allocate(bfrag(4,nbfrag)) !(4,maxres/3)
3947 ! allocate(bfrag(4,nres/3)) !(4,maxres/3)
3951 ! hfrag(j,i)=hfrag_alloc(j,i)
3956 ! bfrag(j,i)=bfrag_alloc(j,i)
3959 !el end reallocate tables
3960 !---------------------------------
3968 c(j,2*nres)=c(j,nres)
3971 if (itype(1,1).eq.ntyp1) then
3975 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3976 call refsys(2,3,4,e1,e2,e3,fail)
3983 ! c(j,1)=c(j,2)-1.9d0*e2(j)
3984 c(j,1)=c(j,2)+1.9d0*(e1(j)-e2(j))/sqrt(2.0d0)
3988 dcj=(c(j,4)-c(j,3))/2.0
3994 ! First lets assign correct dummy to correct type of chain
3996 if (itype(1,1).eq.ntyp1) then
3997 if (itype(2,1).eq.0) then
3999 if (itype(2,j).ne.0) then
4001 itype(1,j)=ntyp1_molec(j)
4002 nres_molec(1)=nres_molec(1)-1
4003 nres_molec(j)=nres_molec(j)+1
4010 print *,'I have',nres, nres_molec(:)
4012 ! Copy the coordinates to reference coordinates
4018 ! Calculate internal coordinates.
4020 write (iout,'(/a)') &
4021 "Cartesian coordinates of the reference structure"
4022 write (iout,'(a,16x,3(3x,a5),5x,3(3x,a5))') &
4023 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
4025 write (iout,'(5(a3,1x),i5,3f8.3,5x,3f8.3)') &
4026 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
4027 (c(j,ires+nres),j=1,3)
4030 ! znamy już nres wiec mozna alokowac tablice
4031 ! Calculate internal coordinates.
4032 if(me.eq.king.or..not.out1file)then
4033 write (iout,'(a)') &
4034 "Backbone and SC coordinates as read from the PDB"
4036 write (iout,'(i5,i3,2x,a,3f8.3,5x,3f8.3)') &
4037 ires,itype(ires,1),restyp(itype(ires,1),1),(c(j,ires),j=1,3),&
4038 (c(j,nres+ires),j=1,3)
4041 ! NOW LETS ROCK! SORTING
4042 allocate(c_temporary(3,2*nres))
4043 allocate(itype_temporary(nres,5))
4044 if (.not.allocated(molnum)) allocate(molnum(nres+1))
4045 if (.not.allocated(istype)) write(iout,*) &
4046 "SOMETHING WRONG WITH ISTYTPE"
4047 allocate(istype_temp(nres))
4048 itype_temporary(:,:)=0
4052 if (itype(i,k).ne.0) then
4054 c_temporary(j,seqalingbegin)=c(j,i)
4055 c_temporary(j,seqalingbegin+nres)=c(j,i+nres)
4058 itype_temporary(seqalingbegin,k)=itype(i,k)
4059 print *,i,k,itype(i,k),itype_temporary(seqalingbegin,k),seqalingbegin
4060 istype_temp(seqalingbegin)=istype(i)
4061 molnum(seqalingbegin)=k
4062 seqalingbegin=seqalingbegin+1
4068 c(j,i)=c_temporary(j,i)
4073 itype(i,k)=itype_temporary(i,k)
4074 istype(i)=istype_temp(i)
4077 if ((itype(1,1).eq.ntyp1).and.itype(2,5).ne.0) then
4078 ! I have only ions now dummy atoms in the system
4080 itype(1,5)=itype(2,5)
4083 itype(i,5)=itype(i+1,5)
4087 nres_molec(1)=nres_molec(1)-1
4089 ! if (itype(1,1).eq.ntyp1) then
4092 ! if (unres_pdb) then
4093 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
4094 ! call refsys(2,3,4,e1,e2,e3,fail)
4101 ! c(j,1)=c(j,2)-1.9d0*e2(j)
4105 ! dcj=(c(j,4)-c(j,3))/2.0
4107 ! c(j,nres+1)=c(j,1)
4113 write (iout,'(/a)') &
4114 "Cartesian coordinates of the reference structure after sorting"
4115 write (iout,'(a,16x,3(3x,a5),5x,3(3x,a5))') &
4116 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
4118 write (iout,'(5(a3,1x),i5,3f8.3,5x,3f8.3)') &
4119 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
4120 (c(j,ires+nres),j=1,3)
4124 ! print *,seqalingbegin,nres
4125 if(.not.allocated(vbld)) then
4126 allocate(vbld(2*nres))
4131 if(.not.allocated(vbld_inv)) then
4132 allocate(vbld_inv(2*nres))
4138 if(.not.allocated(theta)) then
4139 allocate(theta(nres+2))
4143 if(.not.allocated(phi)) allocate(phi(nres+2))
4144 if(.not.allocated(alph)) allocate(alph(nres+2))
4145 if(.not.allocated(omeg)) allocate(omeg(nres+2))
4146 if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
4147 if(.not.allocated(phiref)) allocate(phiref(nres+2))
4148 if(.not.allocated(costtab)) allocate(costtab(nres))
4149 if(.not.allocated(sinttab)) allocate(sinttab(nres))
4150 if(.not.allocated(cost2tab)) allocate(cost2tab(nres))
4151 if(.not.allocated(sint2tab)) allocate(sint2tab(nres))
4152 if(.not.allocated(xxref)) allocate(xxref(nres))
4153 if(.not.allocated(yyref)) allocate(yyref(nres))
4154 if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
4155 if(.not.allocated(dc_norm)) then
4156 ! if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres+2))
4157 allocate(dc_norm(3,0:2*nres+2))
4161 call int_from_cart(.true.,.false.)
4162 call sc_loc_geom(.false.)
4164 thetaref(i)=theta(i)
4174 dc(j,i)=c(j,i+1)-c(j,i)
4175 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
4180 dc(j,i+nres)=c(j,i+nres)-c(j,i)
4181 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
4183 ! write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),&
4187 ! Copy the coordinates to reference coordinates
4188 ! Splits to single chain if occurs
4192 ! cref(j,i,cou)=c(j,i)
4196 if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
4197 if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
4198 if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
4199 !-----------------------------
4203 write (iout,*) "symetr", symetr
4206 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
4208 if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then
4211 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
4216 cref(j,i,cou)=c(j,i)
4217 cref(j,i+nres,cou)=c(j,i+nres)
4219 chain_rep(j,lll,kkk)=c(j,i)
4220 chain_rep(j,lll+nres,kkk)=c(j,i+nres)
4224 write (iout,*) chain_length
4225 if (chain_length.eq.0) chain_length=nres
4227 chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
4228 chain_rep(j,chain_length+nres,symetr) &
4229 =chain_rep(j,chain_length+nres,1)
4232 ! write (iout,*) "spraw lancuchy",chain_length,symetr
4234 ! do kkk=1,chain_length
4235 ! write (iout,*) itype(kkk,1),(chain_rep(j,kkk,i), j=1,3)
4239 ! makes copy of chains
4240 write (iout,*) "symetr", symetr
4245 if (symetr.gt.1) then
4252 write(iout,*) (tabperm(i,kkk),kkk=1,4)
4258 write (iout,*) i,icha
4259 do lll=1,chain_length
4261 if (cou.le.nres) then
4263 kupa=mod(lll,chain_length)
4264 iprzes=(kkk-1)*chain_length+lll
4265 if (kupa.eq.0) kupa=chain_length
4266 write (iout,*) "kupa", kupa
4267 cref(j,iprzes,i)=chain_rep(j,kupa,icha)
4268 cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
4275 !-koniec robienia kopii
4278 write (iout,*) "nowa struktura", nperm
4280 write (iout,110) restyp(itype(i,1),1),i,cref(1,i,kkk),&
4282 cref(3,i,kkk),cref(1,nres+i,kkk),&
4283 cref(2,nres+i,kkk),cref(3,nres+i,kkk)
4285 100 format (//' alpha-carbon coordinates ',&
4286 ' centroid coordinates'/ &
4287 ' ', 6X,'X',11X,'Y',11X,'Z', &
4288 10X,'X',11X,'Y',11X,'Z')
4289 110 format (a,'(',i5,')',6f12.5)
4295 bfrag(i,j)=bfrag(i,j)-ishift
4301 hfrag(i,j)=hfrag(i,j)-ishift
4307 end subroutine readpdb
4308 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
4309 !-----------------------------------------------------------------------------
4311 !-----------------------------------------------------------------------------
4312 subroutine read_control
4326 use random, only: random_init
4327 ! implicit real*8 (a-h,o-z)
4328 ! include 'DIMENSIONS'
4330 use prng, only:prng_restart
4332 logical :: OKRandom!, prng_restart
4335 ! include 'COMMON.IOUNITS'
4336 ! include 'COMMON.TIME1'
4337 ! include 'COMMON.THREAD'
4338 ! include 'COMMON.SBRIDGE'
4339 ! include 'COMMON.CONTROL'
4340 ! include 'COMMON.MCM'
4341 ! include 'COMMON.MAP'
4342 ! include 'COMMON.HEADER'
4343 ! include 'COMMON.CSA'
4344 ! include 'COMMON.CHAIN'
4345 ! include 'COMMON.MUCA'
4346 ! include 'COMMON.MD'
4347 ! include 'COMMON.FFIELD'
4348 ! include 'COMMON.INTERACT'
4349 ! include 'COMMON.SETUP'
4350 !el integer :: KDIAG,ICORFL,IXDR
4351 !el COMMON /MACHSW/ KDIAG,ICORFL,IXDR
4352 character(len=8),dimension(0:3) :: diagmeth = reshape((/'Library ',&
4353 'EVVRSP ','Givens ','Jacobi '/),shape(diagmeth))
4354 ! character(len=80) :: ucase
4355 character(len=640) :: controlcard
4357 real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
4363 read (INP,'(a)') titel
4364 call card_concat(controlcard,.true.)
4365 ! out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
4366 ! print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
4367 call reada(controlcard,'SEED',seed,0.0D0)
4368 call random_init(seed)
4369 ! Set up the time limit (caution! The time must be input in minutes!)
4370 read_cart=index(controlcard,'READ_CART').gt.0
4371 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
4372 call readi(controlcard,'SYM',symetr,1)
4373 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
4374 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
4375 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
4376 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
4377 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
4378 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
4379 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
4380 call reada(controlcard,'DRMS',drms,0.1D0)
4381 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4382 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
4383 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
4384 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
4385 write (iout,'(a,f10.1)')'DRMS = ',drms
4386 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
4387 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
4389 call readi(controlcard,'NZ_START',nz_start,0)
4390 call readi(controlcard,'NZ_END',nz_end,0)
4391 call readi(controlcard,'IZ_SC',iz_sc,0)
4392 timlim=60.0D0*timlim
4393 safety = 60.0d0*safety
4396 call reada(controlcard,"T_BATH",t_bath,300.0d0)
4397 !C SHIELD keyword sets if the shielding effect of side-chains is used
4398 !C 0 denots no shielding is used all peptide are equally despite the
4399 !C solvent accesible area
4400 !C 1 the newly introduced function
4401 !C 2 reseved for further possible developement
4402 call readi(controlcard,'SHIELD',shield_mode,0)
4403 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4404 write(iout,*) "shield_mode",shield_mode
4405 call readi(controlcard,'TORMODE',tor_mode,0)
4406 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4407 write(iout,*) "torsional and valence angle mode",tor_mode
4409 !C Varibles set size of box
4410 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
4411 protein=index(controlcard,"PROTEIN").gt.0
4412 ions=index(controlcard,"IONS").gt.0
4413 call readi(controlcard,'OLDION',oldion,1)
4414 nucleic=index(controlcard,"NUCLEIC").gt.0
4415 write (iout,*) "with_theta_constr ",with_theta_constr
4416 AFMlog=(index(controlcard,'AFM'))
4417 selfguide=(index(controlcard,'SELFGUIDE'))
4418 print *,'AFMlog',AFMlog,selfguide,"KUPA"
4419 call readi(controlcard,'GENCONSTR',genconstr,0)
4420 call reada(controlcard,'BOXX',boxxsize,100.0d0)
4421 call reada(controlcard,'BOXY',boxysize,100.0d0)
4422 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
4423 call readi(controlcard,'TUBEMOD',tubemode,0)
4424 print *,"SCELE",scelemode
4425 call readi(controlcard,"SCELEMODE",scelemode,0)
4426 print *,"SCELE",scelemode
4428 ! elemode = 0 is orignal UNRES electrostatics
4429 ! elemode = 1 is "Momo" potentials in progress
4430 ! elemode = 2 is in development EVALD
4433 write (iout,*) TUBEmode,"TUBEMODE"
4434 if (TUBEmode.gt.0) then
4435 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
4436 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
4437 call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
4438 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
4439 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
4440 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
4441 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
4442 buftubebot=bordtubebot+tubebufthick
4443 buftubetop=bordtubetop-tubebufthick
4446 ! CUTOFFF ON ELECTROSTATICS
4447 call reada(controlcard,"R_CUT_ELE",r_cut_ele,25.0d0)
4448 call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
4449 write(iout,*) "R_CUT_ELE=",r_cut_ele
4450 ! Lipidic parameters
4451 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
4452 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
4453 if (lipthick.gt.0.0d0) then
4454 bordliptop=(boxzsize+lipthick)/2.0
4455 bordlipbot=bordliptop-lipthick
4456 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
4457 write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
4458 buflipbot=bordlipbot+lipbufthick
4459 bufliptop=bordliptop-lipbufthick
4460 if ((lipbufthick*2.0d0).gt.lipthick) &
4461 write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
4462 endif !lipthick.gt.0
4463 write(iout,*) "bordliptop=",bordliptop
4464 write(iout,*) "bordlipbot=",bordlipbot
4465 write(iout,*) "bufliptop=",bufliptop
4466 write(iout,*) "buflipbot=",buflipbot
4467 write (iout,*) "SHIELD MODE",shield_mode
4469 !C-------------------------
4470 minim=(index(controlcard,'MINIMIZE').gt.0)
4471 dccart=(index(controlcard,'CART').gt.0)
4472 overlapsc=(index(controlcard,'OVERLAP').gt.0)
4473 overlapsc=.not.overlapsc
4474 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
4475 searchsc=.not.searchsc
4476 sideadd=(index(controlcard,'SIDEADD').gt.0)
4477 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
4478 outpdb=(index(controlcard,'PDBOUT').gt.0)
4479 outmol2=(index(controlcard,'MOL2OUT').gt.0)
4480 pdbref=(index(controlcard,'PDBREF').gt.0)
4481 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
4482 indpdb=index(controlcard,'PDBSTART')
4483 extconf=(index(controlcard,'EXTCONF').gt.0)
4484 call readi(controlcard,'IPRINT',iprint,0)
4485 call readi(controlcard,'MAXGEN',maxgen,10000)
4486 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
4487 call readi(controlcard,"KDIAG",kdiag,0)
4488 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
4489 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) &
4490 write (iout,*) "RESCALE_MODE",rescale_mode
4491 split_ene=index(controlcard,'SPLIT_ENE').gt.0
4492 if (index(controlcard,'REGULAR').gt.0.0D0) then
4493 call reada(controlcard,'WEIDIS',weidis,0.1D0)
4497 if (index(controlcard,'CHECKGRAD').gt.0) then
4499 if (index(controlcard,'CART').gt.0) then
4501 elseif (index(controlcard,'CARINT').gt.0) then
4506 elseif (index(controlcard,'THREAD').gt.0) then
4508 call readi(controlcard,'THREAD',nthread,0)
4509 if (nthread.gt.0) then
4510 call reada(controlcard,'WEIDIS',weidis,0.1D0)
4513 write (iout,'(a)')'A number has to follow the THREAD keyword.'
4514 stop 'Error termination in Read_Control.'
4516 else if (index(controlcard,'MCMA').gt.0) then
4518 else if (index(controlcard,'MCEE').gt.0) then
4520 else if (index(controlcard,'MULTCONF').gt.0) then
4522 else if (index(controlcard,'MAP').gt.0) then
4524 call readi(controlcard,'MAP',nmap,0)
4525 else if (index(controlcard,'CSA').gt.0) then
4527 !rc else if (index(controlcard,'ZSCORE').gt.0) then
4529 !rc ZSCORE is rm from UNRES, modecalc=9 is available
4532 !fcm else if (index(controlcard,'MCMF').gt.0) then
4534 else if (index(controlcard,'SOFTREG').gt.0) then
4536 else if (index(controlcard,'CHECK_BOND').gt.0) then
4538 else if (index(controlcard,'TEST').gt.0) then
4540 else if (index(controlcard,'MD').gt.0) then
4542 else if (index(controlcard,'RE ').gt.0) then
4546 lmuca=index(controlcard,'MUCA').gt.0
4547 call readi(controlcard,'MUCADYN',mucadyn,0)
4548 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
4549 if (lmuca .and. (me.eq.king .or. .not.out1file )) &
4551 write (iout,*) 'MUCADYN=',mucadyn
4552 write (iout,*) 'MUCASMOOTH=',muca_smooth
4555 iscode=index(controlcard,'ONE_LETTER')
4556 indphi=index(controlcard,'PHI')
4557 indback=index(controlcard,'BACK')
4558 iranconf=index(controlcard,'RAND_CONF')
4559 i2ndstr=index(controlcard,'USE_SEC_PRED')
4560 gradout=index(controlcard,'GRADOUT').gt.0
4561 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
4562 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
4563 if (me.eq.king .or. .not.out1file ) &
4564 write (iout,*) "DISTCHAINMAX",distchainmax
4566 if(me.eq.king.or..not.out1file) &
4567 write (iout,'(2a)') diagmeth(kdiag),&
4568 ' routine used to diagonalize matrices.'
4569 if (shield_mode.gt.0) then
4570 pi=4.0D0*datan(1.0D0)
4571 !C VSolvSphere the volume of solving sphere
4573 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
4574 !C there will be no distinction between proline peptide group and normal peptide
4575 !C group in case of shielding parameters
4576 VSolvSphere=4.0/3.0*pi*(4.50d0)**3
4577 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3
4578 write (iout,*) VSolvSphere,VSolvSphere_div
4579 !C long axis of side chain
4581 ! long_r_sidechain(i)=vbldsc0(1,i)
4582 ! short_r_sidechain(i)=sigma0(i)
4583 ! write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
4589 end subroutine read_control
4590 !-----------------------------------------------------------------------------
4591 subroutine read_REMDpar
4593 ! Read REMD settings
4600 use control_data, only:out1file
4601 ! implicit real*8 (a-h,o-z)
4602 ! include 'DIMENSIONS'
4603 ! include 'COMMON.IOUNITS'
4604 ! include 'COMMON.TIME1'
4605 ! include 'COMMON.MD'
4608 !el include 'COMMON.LANGEVIN'
4610 !el include 'COMMON.LANGEVIN.lang0'
4612 ! include 'COMMON.INTERACT'
4613 ! include 'COMMON.NAMES'
4614 ! include 'COMMON.GEO'
4615 ! include 'COMMON.REMD'
4616 ! include 'COMMON.CONTROL'
4617 ! include 'COMMON.SETUP'
4618 ! character(len=80) :: ucase
4619 character(len=320) :: controlcard
4620 character(len=3200) :: controlcard1
4621 integer :: iremd_m_total
4624 ! real(kind=8) :: var,ene
4626 if(me.eq.king.or..not.out1file) &
4627 write (iout,*) "REMD setup"
4629 call card_concat(controlcard,.true.)
4630 call readi(controlcard,"NREP",nrep,3)
4631 call readi(controlcard,"NSTEX",nstex,1000)
4632 call reada(controlcard,"RETMIN",retmin,10.0d0)
4633 call reada(controlcard,"RETMAX",retmax,1000.0d0)
4634 mremdsync=(index(controlcard,'SYNC').gt.0)
4635 call readi(controlcard,"NSYN",i_sync_step,100)
4636 restart1file=(index(controlcard,'REST1FILE').gt.0)
4637 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
4638 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
4639 if(max_cache_traj_use.gt.max_cache_traj) &
4640 max_cache_traj_use=max_cache_traj
4641 if(me.eq.king.or..not.out1file) then
4642 !d if (traj1file) then
4643 !rc caching is in testing - NTWX is not ignored
4644 !d write (iout,*) "NTWX value is ignored"
4645 !d write (iout,*) " trajectory is stored to one file by master"
4646 !d write (iout,*) " before exchange at NSTEX intervals"
4648 write (iout,*) "NREP= ",nrep
4649 write (iout,*) "NSTEX= ",nstex
4650 write (iout,*) "SYNC= ",mremdsync
4651 write (iout,*) "NSYN= ",i_sync_step
4652 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
4655 allocate(remd_t(nrep),remd_m(nrep)) !(maxprocs)
4656 if (index(controlcard,'TLIST').gt.0) then
4658 call card_concat(controlcard1,.true.)
4659 read(controlcard1,*) (remd_t(i),i=1,nrep)
4660 if(me.eq.king.or..not.out1file) &
4661 write (iout,*)'tlist',(remd_t(i),i=1,nrep)
4664 if (index(controlcard,'MLIST').gt.0) then
4666 call card_concat(controlcard1,.true.)
4667 read(controlcard1,*) (remd_m(i),i=1,nrep)
4668 if(me.eq.king.or..not.out1file) then
4669 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
4672 iremd_m_total=iremd_m_total+remd_m(i)
4674 write (iout,*) 'Total number of replicas ',iremd_m_total
4677 if(me.eq.king.or..not.out1file) &
4678 write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
4680 end subroutine read_REMDpar
4681 !-----------------------------------------------------------------------------
4682 subroutine read_MDpar
4686 use control_data, only: r_cut,rlamb,out1file
4688 use geometry_data, only: pi
4690 ! implicit real*8 (a-h,o-z)
4691 ! include 'DIMENSIONS'
4692 ! include 'COMMON.IOUNITS'
4693 ! include 'COMMON.TIME1'
4694 ! include 'COMMON.MD'
4697 !el include 'COMMON.LANGEVIN'
4699 !el include 'COMMON.LANGEVIN.lang0'
4701 ! include 'COMMON.INTERACT'
4702 ! include 'COMMON.NAMES'
4703 ! include 'COMMON.GEO'
4704 ! include 'COMMON.SETUP'
4705 ! include 'COMMON.CONTROL'
4706 ! include 'COMMON.SPLITELE'
4707 ! character(len=80) :: ucase
4708 character(len=320) :: controlcard
4713 call card_concat(controlcard,.true.)
4714 call readi(controlcard,"NSTEP",n_timestep,1000000)
4715 call readi(controlcard,"NTWE",ntwe,100)
4716 call readi(controlcard,"NTWX",ntwx,1000)
4717 call reada(controlcard,"DT",d_time,1.0d-1)
4718 call reada(controlcard,"DVMAX",dvmax,2.0d1)
4719 call reada(controlcard,"DAMAX",damax,1.0d1)
4720 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
4721 call readi(controlcard,"LANG",lang,0)
4722 RESPA = index(controlcard,"RESPA") .gt. 0
4723 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
4724 ntime_split0=ntime_split
4725 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
4726 ntime_split0=ntime_split
4727 call reada(controlcard,"R_CUT",r_cut,2.0d0)
4728 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
4729 rest = index(controlcard,"REST").gt.0
4730 tbf = index(controlcard,"TBF").gt.0
4731 usampl = index(controlcard,"USAMPL").gt.0
4732 mdpdb = index(controlcard,"MDPDB").gt.0
4733 call reada(controlcard,"T_BATH",t_bath,300.0d0)
4734 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
4735 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
4736 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
4737 if (count_reset_moment.eq.0) count_reset_moment=1000000000
4738 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
4739 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
4740 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
4741 if (count_reset_vel.eq.0) count_reset_vel=1000000000
4742 large = index(controlcard,"LARGE").gt.0
4743 print_compon = index(controlcard,"PRINT_COMPON").gt.0
4744 rattle = index(controlcard,"RATTLE").gt.0
4745 preminim=(index(controlcard,'PREMINIM').gt.0)
4746 write (iout,*) "PREMINIM ",preminim
4747 dccart=(index(controlcard,'CART').gt.0)
4748 if (preminim) call read_minim
4749 ! if performing umbrella sampling, fragments constrained are read from the fragment file
4755 if(me.eq.king.or..not.out1file) then
4757 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
4759 write (iout,'(a)') "The units are:"
4760 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
4761 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",&
4762 " acceleration: angstrom/(48.9 fs)**2"
4763 write (iout,'(a)') "energy: kcal/mol, temperature: K"
4765 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
4766 write (iout,'(a60,f10.5,a)') &
4767 "Initial time step of numerical integration:",d_time,&
4769 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
4771 write (iout,'(2a,i4,a)') &
4772 "A-MTS algorithm used; initial time step for fast-varying",&
4773 " short-range forces split into",ntime_split," steps."
4774 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",&
4775 r_cut," lambda",rlamb
4777 write (iout,'(2a,f10.5)') &
4778 "Maximum acceleration threshold to reduce the time step",&
4779 "/increase split number:",damax
4780 write (iout,'(2a,f10.5)') &
4781 "Maximum predicted energy drift to reduce the timestep",&
4782 "/increase split number:",edriftmax
4783 write (iout,'(a60,f10.5)') &
4784 "Maximum velocity threshold to reduce velocities:",dvmax
4785 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
4786 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
4787 if (rattle) write (iout,'(a60)') &
4788 "Rattle algorithm used to constrain the virtual bonds"
4792 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
4793 call reada(controlcard,"RWAT",rwat,1.4d0)
4794 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
4795 surfarea=index(controlcard,"SURFAREA").gt.0
4796 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
4797 if(me.eq.king.or..not.out1file)then
4798 write (iout,'(/a,$)') "Langevin dynamics calculation"
4800 write (iout,'(a/)') &
4801 " with direct integration of Langevin equations"
4802 else if (lang.eq.2) then
4803 write (iout,'(a/)') " with TINKER stochasic MD integrator"
4804 else if (lang.eq.3) then
4805 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
4806 else if (lang.eq.4) then
4807 write (iout,'(a/)') " in overdamped mode"
4809 write (iout,'(//a,i5)') &
4810 "=========== ERROR: Unknown Langevin dynamics mode:",lang
4813 write (iout,'(a60,f10.5)') "Temperature:",t_bath
4814 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
4815 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
4816 write (iout,'(a60,f10.5)') &
4817 "Scaling factor of the friction forces:",scal_fric
4818 if (surfarea) write (iout,'(2a,i10,a)') &
4819 "Friction coefficients will be scaled by solvent-accessible",&
4820 " surface area every",reset_fricmat," steps."
4822 ! Calculate friction coefficients and bounds of stochastic forces
4823 eta=6*pi*cPoise*etawat
4824 if(me.eq.king.or..not.out1file) &
4825 write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
4828 do j=1,5 !types of molecules
4829 gamp(j)=scal_fric*(pstok(j)+rwat)*eta
4830 stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
4832 allocate(gamsc(ntyp1,5),stdfsc(ntyp1,5)) !(ntyp1)
4833 do j=1,5 !types of molecules
4835 gamsc(i,j)=scal_fric*(restok(i,j)+rwat)*eta
4836 stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
4840 if(me.eq.king.or..not.out1file)then
4841 write (iout,'(/2a/)') &
4842 "Radii of site types and friction coefficients and std's of",&
4843 " stochastic forces of fully exposed sites"
4844 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp(1),stdfp*dsqrt(gamp(1))
4846 write (iout,'(a5,f5.2,2f10.5)') restyp(i,1),restok(i,1),&
4847 gamsc(i,1),stdfsc(i,1)*dsqrt(gamsc(i,1))
4851 if(me.eq.king.or..not.out1file)then
4852 write (iout,'(a)') "Berendsen bath calculation"
4853 write (iout,'(a60,f10.5)') "Temperature:",t_bath
4854 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
4856 write (iout,'(a,i10,a)') "Momenta will be reset at zero every",&
4857 count_reset_moment," steps"
4859 write (iout,'(a,i10,a)') &
4860 "Velocities will be reset at random every",count_reset_vel,&
4864 if(me.eq.king.or..not.out1file) &
4865 write (iout,'(a31)') "Microcanonical mode calculation"
4867 if(me.eq.king.or..not.out1file)then
4868 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
4870 write(iout,*) "MD running with constraints."
4871 write(iout,*) "Equilibration time ", eq_time, " mtus."
4872 write(iout,*) "Constraining ", nfrag," fragments."
4873 write(iout,*) "Length of each fragment, weight and q0:"
4875 write (iout,*) "Set of restraints #",iset
4877 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),&
4878 ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
4880 write(iout,*) "constraints between ", npair, "fragments."
4881 write(iout,*) "constraint pairs, weights and q0:"
4883 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),&
4884 ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
4886 write(iout,*) "angle constraints within ", nfrag_back,&
4887 "backbone fragments."
4888 write(iout,*) "fragment, weights:"
4890 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),&
4891 ifrag_back(2,i,iset),wfrag_back(1,i,iset),&
4892 wfrag_back(2,i,iset),wfrag_back(3,i,iset)
4895 iset=mod(kolor,nset)+1
4898 if(me.eq.king.or..not.out1file) &
4899 write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
4901 end subroutine read_MDpar
4902 !-----------------------------------------------------------------------------
4906 ! implicit real*8 (a-h,o-z)
4907 ! include 'DIMENSIONS'
4908 ! include 'COMMON.MAP'
4909 ! include 'COMMON.IOUNITS'
4910 character(len=3) :: angid(4) = (/'THE','PHI','ALP','OME'/)
4911 character(len=80) :: mapcard !,ucase
4914 ! real(kind=8) :: var,ene
4917 read (inp,'(a)') mapcard
4918 mapcard=ucase(mapcard)
4919 if (index(mapcard,'PHI').gt.0) then
4921 else if (index(mapcard,'THE').gt.0) then
4923 else if (index(mapcard,'ALP').gt.0) then
4925 else if (index(mapcard,'OME').gt.0) then
4928 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
4929 stop 'Error - illegal variable spec in MAP card.'
4931 call readi (mapcard,'RES1',res1(imap),0)
4932 call readi (mapcard,'RES2',res2(imap),0)
4933 if (res1(imap).eq.0) then
4934 res1(imap)=res2(imap)
4935 else if (res2(imap).eq.0) then
4936 res2(imap)=res1(imap)
4938 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
4939 write (iout,'(a)') &
4940 'Error - illegal definition of variable group in MAP.'
4941 stop 'Error - illegal definition of variable group in MAP.'
4943 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
4944 call reada(mapcard,'TO',ang_to(imap),0.0D0)
4945 call readi(mapcard,'NSTEP',nstep(imap),0)
4946 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
4947 write (iout,'(a)') &
4948 'Illegal boundary and/or step size specification in MAP.'
4949 stop 'Illegal boundary and/or step size specification in MAP.'
4953 end subroutine map_read
4954 !-----------------------------------------------------------------------------
4957 use control_data, only: vdisulf
4959 ! implicit real*8 (a-h,o-z)
4960 ! include 'DIMENSIONS'
4961 ! include 'COMMON.IOUNITS'
4962 ! include 'COMMON.GEO'
4963 ! include 'COMMON.CSA'
4964 ! include 'COMMON.BANK'
4965 ! include 'COMMON.CONTROL'
4966 ! character(len=80) :: ucase
4967 character(len=620) :: mcmcard
4969 ! integer :: ntf,ik,iw_pdb
4970 ! real(kind=8) :: var,ene
4972 call card_concat(mcmcard,.true.)
4974 call readi(mcmcard,'NCONF',nconf,50)
4975 call readi(mcmcard,'NADD',nadd,0)
4976 call readi(mcmcard,'JSTART',jstart,1)
4977 call readi(mcmcard,'JEND',jend,1)
4978 call readi(mcmcard,'NSTMAX',nstmax,500000)
4979 call readi(mcmcard,'N0',n0,1)
4980 call readi(mcmcard,'N1',n1,6)
4981 call readi(mcmcard,'N2',n2,4)
4982 call readi(mcmcard,'N3',n3,0)
4983 call readi(mcmcard,'N4',n4,0)
4984 call readi(mcmcard,'N5',n5,0)
4985 call readi(mcmcard,'N6',n6,10)
4986 call readi(mcmcard,'N7',n7,0)
4987 call readi(mcmcard,'N8',n8,0)
4988 call readi(mcmcard,'N9',n9,0)
4989 call readi(mcmcard,'N14',n14,0)
4990 call readi(mcmcard,'N15',n15,0)
4991 call readi(mcmcard,'N16',n16,0)
4992 call readi(mcmcard,'N17',n17,0)
4993 call readi(mcmcard,'N18',n18,0)
4995 vdisulf=(index(mcmcard,'DYNSS').gt.0)
4997 call readi(mcmcard,'NDIFF',ndiff,2)
4998 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
4999 call readi(mcmcard,'IS1',is1,1)
5000 call readi(mcmcard,'IS2',is2,8)
5001 call readi(mcmcard,'NRAN0',nran0,4)
5002 call readi(mcmcard,'NRAN1',nran1,2)
5003 call readi(mcmcard,'IRR',irr,1)
5004 call readi(mcmcard,'NSEED',nseed,20)
5005 call readi(mcmcard,'NTOTAL',ntotal,10000)
5006 call reada(mcmcard,'CUT1',cut1,2.0d0)
5007 call reada(mcmcard,'CUT2',cut2,5.0d0)
5008 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
5009 call readi(mcmcard,'ICMAX',icmax,3)
5010 call readi(mcmcard,'IRESTART',irestart,0)
5011 !!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
5014 call reada(mcmcard,'DELE',dele,20.0d0)
5015 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
5016 call readi(mcmcard,'IREF',iref,0)
5017 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
5018 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
5019 call readi(mcmcard,'NCONF_IN',nconf_in,0)
5020 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
5021 write (iout,*) "NCONF_IN",nconf_in
5023 end subroutine csaread
5024 !-----------------------------------------------------------------------------
5028 use control_data, only: MaxMoveType
5031 ! implicit real*8 (a-h,o-z)
5032 ! include 'DIMENSIONS'
5033 ! include 'COMMON.MCM'
5034 ! include 'COMMON.MCE'
5035 ! include 'COMMON.IOUNITS'
5036 ! character(len=80) :: ucase
5037 character(len=320) :: mcmcard
5040 ! real(kind=8) :: var,ene
5042 call card_concat(mcmcard,.true.)
5043 call readi(mcmcard,'MAXACC',maxacc,100)
5044 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
5045 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
5046 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
5047 call readi(mcmcard,'MAXREPM',maxrepm,200)
5048 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
5049 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
5050 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
5051 call reada(mcmcard,'E_UP',e_up,5.0D0)
5052 call reada(mcmcard,'DELTE',delte,0.1D0)
5053 call readi(mcmcard,'NSWEEP',nsweep,5)
5054 call readi(mcmcard,'NSTEPH',nsteph,0)
5055 call readi(mcmcard,'NSTEPC',nstepc,0)
5056 call reada(mcmcard,'TMIN',tmin,298.0D0)
5057 call reada(mcmcard,'TMAX',tmax,298.0D0)
5058 call readi(mcmcard,'NWINDOW',nwindow,0)
5059 call readi(mcmcard,'PRINT_MC',print_mc,0)
5060 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
5061 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
5062 ent_read=(index(mcmcard,'ENT_READ').gt.0)
5063 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
5064 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
5065 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
5066 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
5067 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
5068 if (nwindow.gt.0) then
5069 allocate(winstart(nwindow)) !!el (maxres)
5070 allocate(winend(nwindow)) !!el
5071 allocate(winlen(nwindow)) !!el
5072 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
5074 winlen(i)=winend(i)-winstart(i)+1
5077 if (tmax.lt.tmin) tmax=tmin
5078 if (tmax.eq.tmin) then
5082 if (nstepc.gt.0 .and. nsteph.gt.0) then
5083 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
5084 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
5086 allocate(sumpro_type(0:MaxMoveType)) !(0:MaxMoveType)
5087 ! Probabilities of different move types
5088 sumpro_type(0)=0.0D0
5089 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
5090 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
5091 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
5092 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
5093 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
5094 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
5095 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
5097 print *,'i',i,' sumprotype',sumpro_type(i)
5098 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
5099 print *,'i',i,' sumprotype',sumpro_type(i)
5102 end subroutine mcmread
5103 !-----------------------------------------------------------------------------
5104 subroutine read_minim
5107 ! implicit real*8 (a-h,o-z)
5108 ! include 'DIMENSIONS'
5109 ! include 'COMMON.MINIM'
5110 ! include 'COMMON.IOUNITS'
5111 ! character(len=80) :: ucase
5112 character(len=320) :: minimcard
5114 ! integer :: ntf,ik,iw_pdb
5115 ! real(kind=8) :: var,ene
5117 call card_concat(minimcard,.true.)
5118 call readi(minimcard,'MAXMIN',maxmin,2000)
5119 call readi(minimcard,'MAXFUN',maxfun,5000)
5120 call readi(minimcard,'MINMIN',minmin,maxmin)
5121 call readi(minimcard,'MINFUN',minfun,maxmin)
5122 call reada(minimcard,'TOLF',tolf,1.0D-2)
5123 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
5124 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
5125 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
5126 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
5127 write (iout,'(/80(1h*)/20x,a/80(1h*))') &
5128 'Options in energy minimization:'
5129 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') &
5130 'MaxMin:',MaxMin,' MaxFun:',MaxFun,&
5131 'MinMin:',MinMin,' MinFun:',MinFun,&
5132 ' TolF:',TolF,' RTolF:',RTolF
5134 end subroutine read_minim
5135 !-----------------------------------------------------------------------------
5136 subroutine openunits
5138 use MD_data, only: usampl
5141 use control_data, only:out1file
5142 use control, only: getenv_loc
5143 ! implicit real*8 (a-h,o-z)
5144 ! include 'DIMENSIONS'
5147 character(len=16) :: form,nodename
5148 integer :: nodelen,ierror,npos
5150 ! include 'COMMON.SETUP'
5151 ! include 'COMMON.IOUNITS'
5152 ! include 'COMMON.MD'
5153 ! include 'COMMON.CONTROL'
5154 integer :: lenpre,lenpot,lentmp !,ilen
5156 character(len=3) :: out1file_text !,ucase
5157 character(len=3) :: ll
5160 ! integer :: ntf,ik,iw_pdb
5161 ! real(kind=8) :: var,ene
5163 ! print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
5164 call getenv_loc("PREFIX",prefix)
5166 call getenv_loc("POT",pot)
5167 call getenv_loc("DIRTMP",tmpdir)
5168 call getenv_loc("CURDIR",curdir)
5169 call getenv_loc("OUT1FILE",out1file_text)
5170 ! print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
5171 out1file_text=ucase(out1file_text)
5172 if (out1file_text(1:1).eq."Y") then
5175 out1file=fg_rank.gt.0
5180 if (lentmp.gt.0) then
5181 write (*,'(80(1h!))')
5182 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
5183 write (*,'(80(1h!))')
5184 write (*,*)"All output files will be on node /tmp directory."
5186 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
5187 if (me.eq.king) then
5188 write (*,*) "The master node is ",nodename
5189 else if (fg_rank.eq.0) then
5190 write (*,*) "I am the CG slave node ",nodename
5192 write (*,*) "I am the FG slave node ",nodename
5195 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
5196 lenpre = lentmp+lenpre+1
5198 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
5199 ! Get the names and open the input files
5200 #if defined(WINIFL) || defined(WINPGI)
5201 open(1,file=pref_orig(:ilen(pref_orig))// &
5202 '.inp',status='old',readonly,shared)
5203 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5204 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5205 ! Get parameter filenames and open the parameter files.
5206 call getenv_loc('BONDPAR',bondname)
5207 open (ibond,file=bondname,status='old',readonly,shared)
5208 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5209 open (ibond_nucl,file=bondname_nucl,status='old',readonly,shared)
5210 call getenv_loc('THETPAR',thetname)
5211 open (ithep,file=thetname,status='old',readonly,shared)
5212 call getenv_loc('ROTPAR',rotname)
5213 open (irotam,file=rotname,status='old',readonly,shared)
5214 call getenv_loc('TORPAR',torname)
5215 open (itorp,file=torname,status='old',readonly,shared)
5216 call getenv_loc('TORDPAR',tordname)
5217 open (itordp,file=tordname,status='old',readonly,shared)
5218 call getenv_loc('FOURIER',fouriername)
5219 open (ifourier,file=fouriername,status='old',readonly,shared)
5220 call getenv_loc('ELEPAR',elename)
5221 open (ielep,file=elename,status='old',readonly,shared)
5222 call getenv_loc('SIDEPAR',sidename)
5223 open (isidep,file=sidename,status='old',readonly,shared)
5225 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5226 open (ithep_nucl,file=thetname_nucl,status='old',readonly,shared)
5227 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5228 open (irotam_nucl,file=rotname_nucl,status='old',readonly,shared)
5229 call getenv_loc('TORPAR_NUCL',torname_nucl)
5230 open (itorp_nucl,file=torname_nucl,status='old',readonly,shared)
5231 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5232 open (itordp_nucl,file=tordname_nucl,status='old',readonly,shared)
5233 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5234 open (isidep_nucl,file=sidename_nucl,status='old',readonly,shared)
5237 #elif (defined CRAY) || (defined AIX)
5238 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
5240 ! print *,"Processor",myrank," opened file 1"
5241 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5242 ! print *,"Processor",myrank," opened file 9"
5243 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5244 ! Get parameter filenames and open the parameter files.
5245 call getenv_loc('BONDPAR',bondname)
5246 open (ibond,file=bondname,status='old',action='read')
5247 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5248 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
5250 ! print *,"Processor",myrank," opened file IBOND"
5251 call getenv_loc('THETPAR',thetname)
5252 open (ithep,file=thetname,status='old',action='read')
5253 ! print *,"Processor",myrank," opened file ITHEP"
5254 call getenv_loc('ROTPAR',rotname)
5255 open (irotam,file=rotname,status='old',action='read')
5256 ! print *,"Processor",myrank," opened file IROTAM"
5257 call getenv_loc('TORPAR',torname)
5258 open (itorp,file=torname,status='old',action='read')
5259 ! print *,"Processor",myrank," opened file ITORP"
5260 call getenv_loc('TORDPAR',tordname)
5261 open (itordp,file=tordname,status='old',action='read')
5262 ! print *,"Processor",myrank," opened file ITORDP"
5263 call getenv_loc('SCCORPAR',sccorname)
5264 open (isccor,file=sccorname,status='old',action='read')
5265 ! print *,"Processor",myrank," opened file ISCCOR"
5266 call getenv_loc('FOURIER',fouriername)
5267 open (ifourier,file=fouriername,status='old',action='read')
5268 ! print *,"Processor",myrank," opened file IFOURIER"
5269 call getenv_loc('ELEPAR',elename)
5270 open (ielep,file=elename,status='old',action='read')
5271 ! print *,"Processor",myrank," opened file IELEP"
5272 call getenv_loc('SIDEPAR',sidename)
5273 open (isidep,file=sidename,status='old',action='read')
5275 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5276 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
5277 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5278 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
5279 call getenv_loc('TORPAR_NUCL',torname_nucl)
5280 open (itorp_nucl,file=torname_nucl,status='old',action='read')
5281 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5282 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
5283 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5284 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
5286 call getenv_loc('LIPTRANPAR',liptranname)
5287 open (iliptranpar,file=liptranname,status='old',action='read')
5288 call getenv_loc('TUBEPAR',tubename)
5289 open (itube,file=tubename,status='old',action='read')
5290 call getenv_loc('IONPAR',ionname)
5291 open (iion,file=ionname,status='old',action='read')
5293 ! print *,"Processor",myrank," opened file ISIDEP"
5294 ! print *,"Processor",myrank," opened parameter files"
5296 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
5297 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5298 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5299 ! Get parameter filenames and open the parameter files.
5300 call getenv_loc('BONDPAR',bondname)
5301 open (ibond,file=bondname,status='old')
5302 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5303 open (ibond_nucl,file=bondname_nucl,status='old')
5305 call getenv_loc('THETPAR',thetname)
5306 open (ithep,file=thetname,status='old')
5307 call getenv_loc('ROTPAR',rotname)
5308 open (irotam,file=rotname,status='old')
5309 call getenv_loc('TORPAR',torname)
5310 open (itorp,file=torname,status='old')
5311 call getenv_loc('TORDPAR',tordname)
5312 open (itordp,file=tordname,status='old')
5313 call getenv_loc('SCCORPAR',sccorname)
5314 open (isccor,file=sccorname,status='old')
5315 call getenv_loc('FOURIER',fouriername)
5316 open (ifourier,file=fouriername,status='old')
5317 call getenv_loc('ELEPAR',elename)
5318 open (ielep,file=elename,status='old')
5319 call getenv_loc('SIDEPAR',sidename)
5320 open (isidep,file=sidename,status='old')
5322 open (ithep_nucl,file=thetname_nucl,status='old')
5323 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5324 open (irotam_nucl,file=rotname_nucl,status='old')
5325 call getenv_loc('TORPAR_NUCL',torname_nucl)
5326 open (itorp_nucl,file=torname_nucl,status='old')
5327 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5328 open (itordp_nucl,file=tordname_nucl,status='old')
5329 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5330 open (isidep_nucl,file=sidename_nucl,status='old')
5332 call getenv_loc('LIPTRANPAR',liptranname)
5333 open (iliptranpar,file=liptranname,status='old')
5334 call getenv_loc('TUBEPAR',tubename)
5335 open (itube,file=tubename,status='old')
5336 call getenv_loc('IONPAR',ionname)
5337 open (iion,file=ionname,status='old')
5339 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
5341 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5342 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5343 ! Get parameter filenames and open the parameter files.
5344 call getenv_loc('BONDPAR',bondname)
5345 open (ibond,file=bondname,status='old',action='read')
5346 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5347 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
5348 call getenv_loc('THETPAR',thetname)
5349 open (ithep,file=thetname,status='old',action='read')
5350 call getenv_loc('ROTPAR',rotname)
5351 open (irotam,file=rotname,status='old',action='read')
5352 call getenv_loc('TORPAR',torname)
5353 open (itorp,file=torname,status='old',action='read')
5354 call getenv_loc('TORDPAR',tordname)
5355 open (itordp,file=tordname,status='old',action='read')
5356 call getenv_loc('SCCORPAR',sccorname)
5357 open (isccor,file=sccorname,status='old',action='read')
5359 call getenv_loc('THETPARPDB',thetname_pdb)
5360 print *,"thetname_pdb ",thetname_pdb
5361 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
5362 print *,ithep_pdb," opened"
5364 call getenv_loc('FOURIER',fouriername)
5365 open (ifourier,file=fouriername,status='old',readonly)
5366 call getenv_loc('ELEPAR',elename)
5367 open (ielep,file=elename,status='old',readonly)
5368 call getenv_loc('SIDEPAR',sidename)
5369 open (isidep,file=sidename,status='old',readonly)
5371 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5372 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
5373 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5374 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
5375 call getenv_loc('TORPAR_NUCL',torname_nucl)
5376 open (itorp_nucl,file=torname_nucl,status='old',action='read')
5377 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5378 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
5379 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5380 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
5381 call getenv_loc('SIDEPAR_SCBASE',sidename_scbase)
5382 open (isidep_scbase,file=sidename_scbase,status='old',action='read')
5383 call getenv_loc('PEPPAR_PEPBASE',pepname_pepbase)
5384 open (isidep_pepbase,file=pepname_pepbase,status='old',action='read')
5385 call getenv_loc('SCPAR_PHOSPH',pepname_scpho)
5386 open (isidep_scpho,file=pepname_scpho,status='old',action='read')
5387 call getenv_loc('PEPPAR_PHOSPH',pepname_peppho)
5388 open (isidep_peppho,file=pepname_peppho,status='old',action='read')
5391 call getenv_loc('LIPTRANPAR',liptranname)
5392 open (iliptranpar,file=liptranname,status='old',action='read')
5393 call getenv_loc('TUBEPAR',tubename)
5394 open (itube,file=tubename,status='old',action='read')
5395 call getenv_loc('IONPAR',ionname)
5396 open (iion,file=ionname,status='old',action='read')
5399 call getenv_loc('ROTPARPDB',rotname_pdb)
5400 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
5403 call getenv_loc('SCPPAR_NUCL',scpname_nucl)
5404 #if defined(WINIFL) || defined(WINPGI)
5405 open (iscpp_nucl,file=scpname_nucl,status='old',readonly,shared)
5406 #elif (defined CRAY) || (defined AIX)
5407 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
5409 open (iscpp_nucl,file=scpname_nucl,status='old')
5411 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
5416 ! 8/9/01 In the newest version SCp interaction constants are read from a file
5417 ! Use -DOLDSCP to use hard-coded constants instead.
5419 call getenv_loc('SCPPAR',scpname)
5420 #if defined(WINIFL) || defined(WINPGI)
5421 open (iscpp,file=scpname,status='old',readonly,shared)
5422 #elif (defined CRAY) || (defined AIX)
5423 open (iscpp,file=scpname,status='old',action='read')
5425 open (iscpp,file=scpname,status='old')
5427 open (iscpp,file=scpname,status='old',action='read')
5430 call getenv_loc('PATTERN',patname)
5431 #if defined(WINIFL) || defined(WINPGI)
5432 open (icbase,file=patname,status='old',readonly,shared)
5433 #elif (defined CRAY) || (defined AIX)
5434 open (icbase,file=patname,status='old',action='read')
5436 open (icbase,file=patname,status='old')
5438 open (icbase,file=patname,status='old',action='read')
5441 ! Open output file only for CG processes
5442 ! print *,"Processor",myrank," fg_rank",fg_rank
5443 if (fg_rank.eq.0) then
5445 if (nodes.eq.1) then
5448 npos = dlog10(dfloat(nodes-1))+1
5450 if (npos.lt.3) npos=3
5451 write (liczba,'(i1)') npos
5452 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba)) &
5454 write (liczba,form) me
5455 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)// &
5456 liczba(:ilen(liczba))
5457 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
5459 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
5461 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)// &
5462 liczba(:ilen(liczba))//'.mol2'
5463 statname=prefix(:lenpre)//'_'//pot(:lenpot)// &
5464 liczba(:ilen(liczba))//'.stat'
5466 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) &
5467 //liczba(:ilen(liczba))//'.stat')
5468 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba)) &
5471 qname=prefix(:lenpre)//'_'//pot(:lenpot)// &
5472 liczba(:ilen(liczba))//'.const'
5477 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
5478 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
5479 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
5480 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
5481 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
5483 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)// &
5485 rest2name=prefix(:ilen(prefix))//'.rst'
5487 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
5490 #if defined(AIX) || defined(PGI)
5491 if (me.eq.king .or. .not. out1file) &
5492 open(iout,file=outname,status='unknown')
5494 if (fg_rank.gt.0) then
5495 write (liczba,'(i3.3)') myrank/nfgtasks
5496 write (ll,'(bz,i3.3)') fg_rank
5497 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
5502 open(igeom,file=intname,status='unknown',position='append')
5503 open(ipdb,file=pdbname,status='unknown')
5504 open(imol2,file=mol2name,status='unknown')
5505 open(istat,file=statname,status='unknown',position='append')
5507 !1out open(iout,file=outname,status='unknown')
5510 if (me.eq.king .or. .not.out1file) &
5511 open(iout,file=outname,status='unknown')
5513 if (fg_rank.gt.0) then
5514 write (liczba,'(i3.3)') myrank/nfgtasks
5515 write (ll,'(bz,i3.3)') fg_rank
5516 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
5521 open(igeom,file=intname,status='unknown',access='append')
5522 open(ipdb,file=pdbname,status='unknown')
5523 open(imol2,file=mol2name,status='unknown')
5524 open(istat,file=statname,status='unknown',access='append')
5526 !1out open(iout,file=outname,status='unknown')
5529 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
5530 csa_seed=prefix(:lenpre)//'.CSA.seed'
5531 csa_history=prefix(:lenpre)//'.CSA.history'
5532 csa_bank=prefix(:lenpre)//'.CSA.bank'
5533 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
5534 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
5535 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
5536 !!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
5537 csa_int=prefix(:lenpre)//'.int'
5538 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
5539 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
5540 csa_in=prefix(:lenpre)//'.CSA.in'
5541 ! print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
5544 write (iout,'(80(1h-))')
5545 write (iout,'(30x,a)') "FILE ASSIGNMENT"
5546 write (iout,'(80(1h-))')
5547 write (iout,*) "Input file : ",&
5548 pref_orig(:ilen(pref_orig))//'.inp'
5549 write (iout,*) "Output file : ",&
5550 outname(:ilen(outname))
5552 write (iout,*) "Sidechain potential file : ",&
5553 sidename(:ilen(sidename))
5555 write (iout,*) "SCp potential file : ",&
5556 scpname(:ilen(scpname))
5558 write (iout,*) "Electrostatic potential file : ",&
5559 elename(:ilen(elename))
5560 write (iout,*) "Cumulant coefficient file : ",&
5561 fouriername(:ilen(fouriername))
5562 write (iout,*) "Torsional parameter file : ",&
5563 torname(:ilen(torname))
5564 write (iout,*) "Double torsional parameter file : ",&
5565 tordname(:ilen(tordname))
5566 write (iout,*) "SCCOR parameter file : ",&
5567 sccorname(:ilen(sccorname))
5568 write (iout,*) "Bond & inertia constant file : ",&
5569 bondname(:ilen(bondname))
5570 write (iout,*) "Bending parameter file : ",&
5571 thetname(:ilen(thetname))
5572 write (iout,*) "Rotamer parameter file : ",&
5573 rotname(:ilen(rotname))
5576 write (iout,*) "Thetpdb parameter file : ",&
5577 thetname_pdb(:ilen(thetname_pdb))
5580 write (iout,*) "Threading database : ",&
5581 patname(:ilen(patname))
5583 write (iout,*)" DIRTMP : ",&
5585 write (iout,'(80(1h-))')
5588 end subroutine openunits
5589 !-----------------------------------------------------------------------------
5592 use geometry_data, only: nres,dc
5594 ! implicit real*8 (a-h,o-z)
5595 ! include 'DIMENSIONS'
5596 ! include 'COMMON.CHAIN'
5597 ! include 'COMMON.IOUNITS'
5598 ! include 'COMMON.MD'
5601 ! real(kind=8) :: var,ene
5603 open(irest2,file=rest2name,status='unknown')
5604 read(irest2,*) totT,EK,potE,totE,t_bath
5607 ! AL 4/17/17: Now reading d_t(0,:) too
5609 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
5612 ! AL 4/17/17: Now reading d_c(0,:) too
5614 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
5617 read (irest2,*) iset
5621 end subroutine readrst
5622 !-----------------------------------------------------------------------------
5623 subroutine read_fragments
5627 use control_data, only:out1file
5630 ! implicit real*8 (a-h,o-z)
5631 ! include 'DIMENSIONS'
5635 ! include 'COMMON.SETUP'
5636 ! include 'COMMON.CHAIN'
5637 ! include 'COMMON.IOUNITS'
5638 ! include 'COMMON.MD'
5639 ! include 'COMMON.CONTROL'
5642 ! real(kind=8) :: var,ene
5644 read(inp,*) nset,nfrag,npair,nfrag_back
5646 !el from module energy
5647 ! if(.not.allocated(mset)) allocate(mset(nset)) !(maxprocs/20)
5648 if(.not.allocated(wfrag_back)) then
5649 allocate(wfrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
5650 allocate(ifrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
5652 allocate(qinfrag(nfrag,nset),wfrag(nfrag,nset)) !(50,maxprocs/20)
5653 allocate(ifrag(2,nfrag,nset)) !(2,50,maxprocs/20)
5655 allocate(qinpair(npair,nset),wpair(npair,nset)) !(100,maxprocs/20)
5656 allocate(ipair(2,npair,nset)) !(2,100,maxprocs/20)
5659 if(me.eq.king.or..not.out1file) &
5660 write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,&
5661 " nfrag_back",nfrag_back
5663 read(inp,*) mset(iset)
5665 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),&
5667 if(me.eq.king.or..not.out1file) &
5668 write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),&
5669 ifrag(2,i,iset), qinfrag(i,iset)
5672 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),&
5674 if(me.eq.king.or..not.out1file) &
5675 write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),&
5676 ipair(2,i,iset), qinpair(i,iset)
5679 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
5680 wfrag_back(3,i,iset),&
5681 ifrag_back(1,i,iset),ifrag_back(2,i,iset)
5682 if(me.eq.king.or..not.out1file) &
5683 write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
5684 wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
5688 end subroutine read_fragments
5689 !-----------------------------------------------------------------------------
5691 !-----------------------------------------------------------------------------
5695 ! implicit real*8 (a-h,o-z)
5696 ! include 'DIMENSIONS'
5697 ! include 'COMMON.CSA'
5698 ! include 'COMMON.BANK'
5699 ! include 'COMMON.IOUNITS'
5701 ! integer :: ntf,ik,iw_pdb
5702 ! real(kind=8) :: var,ene
5704 open(icsa_in,file=csa_in,status="old",err=100)
5705 read(icsa_in,*) nconf
5706 read(icsa_in,*) jstart,jend
5707 read(icsa_in,*) nstmax
5708 read(icsa_in,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
5709 read(icsa_in,*) nran0,nran1,irr
5710 read(icsa_in,*) nseed
5711 read(icsa_in,*) ntotal,cut1,cut2
5712 read(icsa_in,*) estop
5713 read(icsa_in,*) icmax,irestart
5714 read(icsa_in,*) ntbankm,dele,difcut
5715 read(icsa_in,*) iref,rmscut,pnccut
5716 read(icsa_in,*) ndiff
5723 end subroutine csa_read
5724 !-----------------------------------------------------------------------------
5725 subroutine initial_write
5728 ! implicit real*8 (a-h,o-z)
5729 ! include 'DIMENSIONS'
5730 ! include 'COMMON.CSA'
5731 ! include 'COMMON.BANK'
5732 ! include 'COMMON.IOUNITS'
5734 ! integer :: ntf,ik,iw_pdb
5735 ! real(kind=8) :: var,ene
5737 open(icsa_seed,file=csa_seed,status="unknown")
5738 write(icsa_seed,*) "seed"
5740 #if defined(AIX) || defined(PGI)
5741 open(icsa_history,file=csa_history,status="unknown",&
5744 open(icsa_history,file=csa_history,status="unknown",&
5747 write(icsa_history,*) nconf
5748 write(icsa_history,*) jstart,jend
5749 write(icsa_history,*) nstmax
5750 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
5751 write(icsa_history,*) nran0,nran1,irr
5752 write(icsa_history,*) nseed
5753 write(icsa_history,*) ntotal,cut1,cut2
5754 write(icsa_history,*) estop
5755 write(icsa_history,*) icmax,irestart
5756 write(icsa_history,*) ntbankm,dele,difcut
5757 write(icsa_history,*) iref,rmscut,pnccut
5758 write(icsa_history,*) ndiff
5760 write(icsa_history,*)
5763 open(icsa_bank1,file=csa_bank1,status="unknown")
5764 write(icsa_bank1,*) 0
5768 end subroutine initial_write
5769 !-----------------------------------------------------------------------------
5770 subroutine restart_write
5773 ! implicit real*8 (a-h,o-z)
5774 ! include 'DIMENSIONS'
5775 ! include 'COMMON.IOUNITS'
5776 ! include 'COMMON.CSA'
5777 ! include 'COMMON.BANK'
5779 ! integer :: ntf,ik,iw_pdb
5780 ! real(kind=8) :: var,ene
5782 #if defined(AIX) || defined(PGI)
5783 open(icsa_history,file=csa_history,position="append")
5785 open(icsa_history,file=csa_history,access="append")
5787 write(icsa_history,*)
5788 write(icsa_history,*) "This is restart"
5789 write(icsa_history,*)
5790 write(icsa_history,*) nconf
5791 write(icsa_history,*) jstart,jend
5792 write(icsa_history,*) nstmax
5793 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
5794 write(icsa_history,*) nran0,nran1,irr
5795 write(icsa_history,*) nseed
5796 write(icsa_history,*) ntotal,cut1,cut2
5797 write(icsa_history,*) estop
5798 write(icsa_history,*) icmax,irestart
5799 write(icsa_history,*) ntbankm,dele,difcut
5800 write(icsa_history,*) iref,rmscut,pnccut
5801 write(icsa_history,*) ndiff
5802 write(icsa_history,*)
5803 write(icsa_history,*) "irestart is: ", irestart
5805 write(icsa_history,*)
5809 end subroutine restart_write
5810 !-----------------------------------------------------------------------------
5812 !-----------------------------------------------------------------------------
5813 subroutine write_pdb(npdb,titelloc,ee)
5815 ! implicit real*8 (a-h,o-z)
5816 ! include 'DIMENSIONS'
5817 ! include 'COMMON.IOUNITS'
5818 character(len=50) :: titelloc1
5819 character*(*) :: titelloc
5820 character(len=3) :: zahl
5821 character(len=5) :: liczba5
5823 integer :: npdb !,ilen
5827 ! real(kind=8) :: var,ene
5831 if (npdb.lt.1000) then
5832 call numstr(npdb,zahl)
5833 open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
5835 if (npdb.lt.10000) then
5836 write(liczba5,'(i1,i4)') 0,npdb
5838 write(liczba5,'(i5)') npdb
5840 open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
5842 call pdbout(ee,titelloc1,ipdb)
5845 end subroutine write_pdb
5846 !-----------------------------------------------------------------------------
5848 !-----------------------------------------------------------------------------
5849 subroutine write_thread_summary
5850 ! Thread the sequence through a database of known structures
5851 use control_data, only: refstr
5853 use energy_data, only: n_ene_comp
5855 ! implicit real*8 (a-h,o-z)
5856 ! include 'DIMENSIONS'
5858 use MPI_data !include 'COMMON.INFO'
5861 ! include 'COMMON.CONTROL'
5862 ! include 'COMMON.CHAIN'
5863 ! include 'COMMON.DBASE'
5864 ! include 'COMMON.INTERACT'
5865 ! include 'COMMON.VAR'
5866 ! include 'COMMON.THREAD'
5867 ! include 'COMMON.FFIELD'
5868 ! include 'COMMON.SBRIDGE'
5869 ! include 'COMMON.HEADER'
5870 ! include 'COMMON.NAMES'
5871 ! include 'COMMON.IOUNITS'
5872 ! include 'COMMON.TIME1'
5874 integer,dimension(maxthread) :: ip
5875 real(kind=8),dimension(0:n_ene) :: energia
5877 integer :: i,j,ii,jj,ipj,ik,kk,ist
5878 real(kind=8) :: enet,etot,rmsnat,rms,frac,frac_nn
5880 write (iout,'(30x,a/)') &
5881 ' *********** Summary threading statistics ************'
5882 write (iout,'(a)') 'Initial energies:'
5883 write (iout,'(a4,2x,a12,14a14,3a8)') &
5884 'No','seq',(ename(print_order(i)),i=1,nprint_ene),'ETOT',&
5885 'RMSnat','NatCONT','NNCONT','RMS'
5886 ! Energy sort patterns
5891 enet=ener(n_ene-1,ip(i))
5894 if (ener(n_ene-1,ip(j)).lt.enet) then
5896 enet=ener(n_ene-1,ip(j))
5908 ist=nres_base(2,ii)+ipatt(2,i)
5910 energia(i)=ener0(kk,i)
5912 etot=ener0(n_ene_comp+1,i)
5913 rmsnat=ener0(n_ene_comp+2,i)
5914 rms=ener0(n_ene_comp+3,i)
5915 frac=ener0(n_ene_comp+4,i)
5916 frac_nn=ener0(n_ene_comp+5,i)
5919 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5920 i,str_nam(ii),ist+1,&
5921 (energia(print_order(kk)),kk=1,nprint_ene),&
5922 etot,rmsnat,frac,frac_nn,rms
5924 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3)') &
5925 i,str_nam(ii),ist+1,&
5926 (energia(print_order(kk)),kk=1,nprint_ene),etot
5929 write (iout,'(//a)') 'Final energies:'
5930 write (iout,'(a4,2x,a12,17a14,3a8)') &
5931 'No','seq',(ename(print_order(kk)),kk=1,nprint_ene),'ETOT',&
5932 'RMSnat','NatCONT','NNCONT','RMS'
5936 ist=nres_base(2,ii)+ipatt(2,i)
5938 energia(kk)=ener(kk,ik)
5940 etot=ener(n_ene_comp+1,i)
5941 rmsnat=ener(n_ene_comp+2,i)
5942 rms=ener(n_ene_comp+3,i)
5943 frac=ener(n_ene_comp+4,i)
5944 frac_nn=ener(n_ene_comp+5,i)
5945 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5946 i,str_nam(ii),ist+1,&
5947 (energia(print_order(kk)),kk=1,nprint_ene),&
5948 etot,rmsnat,frac,frac_nn,rms
5950 write (iout,'(/a/)') 'IEXAM array:'
5951 write (iout,'(i5)') nexcl
5953 write (iout,'(2i5)') iexam(1,i),iexam(2,i)
5955 write (iout,'(/a,1pe14.4/a,1pe14.4/)') &
5956 'Max. time for threading step ',max_time_for_thread,&
5957 'Average time for threading step: ',ave_time_for_thread
5959 end subroutine write_thread_summary
5960 !-----------------------------------------------------------------------------
5961 subroutine write_stat_thread(ithread,ipattern,ist)
5963 use energy_data, only: n_ene_comp
5965 ! implicit real*8 (a-h,o-z)
5966 ! include "DIMENSIONS"
5967 ! include "COMMON.CONTROL"
5968 ! include "COMMON.IOUNITS"
5969 ! include "COMMON.THREAD"
5970 ! include "COMMON.FFIELD"
5971 ! include "COMMON.DBASE"
5972 ! include "COMMON.NAMES"
5973 real(kind=8),dimension(0:n_ene) :: energia
5975 integer :: ithread,ipattern,ist,i
5976 real(kind=8) :: etot,rmsnat,rms,frac,frac_nn
5978 #if defined(AIX) || defined(PGI)
5979 open(istat,file=statname,position='append')
5981 open(istat,file=statname,access='append')
5984 energia(i)=ener(i,ithread)
5986 etot=ener(n_ene_comp+1,ithread)
5987 rmsnat=ener(n_ene_comp+2,ithread)
5988 rms=ener(n_ene_comp+3,ithread)
5989 frac=ener(n_ene_comp+4,ithread)
5990 frac_nn=ener(n_ene_comp+5,ithread)
5991 write (istat,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5992 ithread,str_nam(ipattern),ist+1,&
5993 (energia(print_order(i)),i=1,nprint_ene),&
5994 etot,rmsnat,frac,frac_nn,rms
5997 end subroutine write_stat_thread
5998 !-----------------------------------------------------------------------------
6000 !-----------------------------------------------------------------------------
6001 end module io_config