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)
871 read(iion,*) msc(i,5),restok(i,5)
872 print *,msc(i,5),restok(i,5)
876 read (iion,*) ncatprotparm
877 allocate(catprm(ncatprotparm,4))
879 read (iion,*) (catprm(i,k),i=1,ncatprotparm)
882 ! read (iion,*) (vcatprm(k),k=1,ncatprotpram)
883 !----------------------------------------------------
884 allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
885 allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp)) !(-ntyp:ntyp)
886 allocate(athet(2,-ntyp:ntyp,-1:1,-1:1))
887 allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
888 allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
889 allocate(gthet(3,-ntyp:ntyp)) !(3,-ntyp:ntyp)
899 allocate(liptranene(ntyp))
900 !C reading lipid parameters
901 write (iout,*) "iliptranpar",iliptranpar
903 read(iliptranpar,*) pepliptran
906 read(iliptranpar,*) liptranene(i)
907 print *,liptranene(i)
913 ! Read the parameters of the probability distribution/energy expression
914 ! of the virtual-bond valence angles theta
917 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),&
918 (bthet(j,i,1,1),j=1,2)
919 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
920 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
921 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
925 athet(1,i,1,-1)=athet(1,i,1,1)
926 athet(2,i,1,-1)=athet(2,i,1,1)
927 bthet(1,i,1,-1)=-bthet(1,i,1,1)
928 bthet(2,i,1,-1)=-bthet(2,i,1,1)
929 athet(1,i,-1,1)=-athet(1,i,1,1)
930 athet(2,i,-1,1)=-athet(2,i,1,1)
931 bthet(1,i,-1,1)=bthet(1,i,1,1)
932 bthet(2,i,-1,1)=bthet(2,i,1,1)
936 athet(1,i,-1,-1)=athet(1,-i,1,1)
937 athet(2,i,-1,-1)=-athet(2,-i,1,1)
938 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
939 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
940 athet(1,i,-1,1)=athet(1,-i,1,1)
941 athet(2,i,-1,1)=-athet(2,-i,1,1)
942 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
943 bthet(2,i,-1,1)=bthet(2,-i,1,1)
944 athet(1,i,1,-1)=-athet(1,-i,1,1)
945 athet(2,i,1,-1)=athet(2,-i,1,1)
946 bthet(1,i,1,-1)=bthet(1,-i,1,1)
947 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
952 polthet(j,i)=polthet(j,-i)
955 gthet(j,i)=gthet(j,-i)
963 'Parameters of the virtual-bond valence angles:'
964 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',&
965 ' ATHETA0 ',' A1 ',' A2 ',&
968 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
969 a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
971 write (iout,'(/a/9x,5a/79(1h-))') &
972 'Parameters of the expression for sigma(theta_c):',&
973 ' ALPH0 ',' ALPH1 ',' ALPH2 ',&
974 ' ALPH3 ',' SIGMA0C '
976 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
977 (polthet(j,i),j=0,3),sigc0(i)
979 write (iout,'(/a/9x,5a/79(1h-))') &
980 'Parameters of the second gaussian:',&
981 ' THETA0 ',' SIGMA0 ',' G1 ',&
984 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,theta0(i),&
985 sig0(i),(gthet(j,i),j=1,3)
989 'Parameters of the virtual-bond valence angles:'
990 write (iout,'(/a/9x,5a/79(1h-))') &
991 'Coefficients of expansion',&
992 ' theta0 ',' a1*10^2 ',' a2*10^2 ',&
993 ' b1*10^1 ',' b2*10^1 '
995 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),&
996 a0thet(i),(100*athet(j,i,1,1),j=1,2),&
997 (10*bthet(j,i,1,1),j=1,2)
999 write (iout,'(/a/9x,5a/79(1h-))') &
1000 'Parameters of the expression for sigma(theta_c):',&
1001 ' alpha0 ',' alph1 ',' alph2 ',&
1002 ' alhp3 ',' sigma0c '
1004 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i,1),&
1005 (polthet(j,i),j=0,3),sigc0(i)
1007 write (iout,'(/a/9x,5a/79(1h-))') &
1008 'Parameters of the second gaussian:',&
1009 ' theta0 ',' sigma0*10^2 ',' G1*10^-1',&
1012 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),theta0(i),&
1013 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
1019 ! Read the parameters of Utheta determined from ab initio surfaces
1020 ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
1022 IF (tor_mode.eq.0) THEN
1023 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,&
1024 ntheterm3,nsingle,ndouble
1025 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
1027 !----------------------------------------------------
1028 allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1029 allocate(aa0thet(-nthetyp-1:nthetyp+1,&
1030 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1031 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1032 allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,&
1033 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1034 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1035 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1036 allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1037 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1038 allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1039 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1040 allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1041 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1042 allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1043 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1044 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1045 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1046 allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1047 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1048 allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1049 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1050 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1051 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1053 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
1055 ithetyp(i)=-ithetyp(-i)
1058 aa0thet(:,:,:,:)=0.0d0
1059 aathet(:,:,:,:,:)=0.0d0
1060 bbthet(:,:,:,:,:,:)=0.0d0
1061 ccthet(:,:,:,:,:,:)=0.0d0
1062 ddthet(:,:,:,:,:,:)=0.0d0
1063 eethet(:,:,:,:,:,:)=0.0d0
1064 ffthet(:,:,:,:,:,:,:)=0.0d0
1065 ggthet(:,:,:,:,:,:,:)=0.0d0
1067 ! VAR:iblock means terminally blocking group 1=non-proline 2=proline
1069 ! VAR:ntethtyp is type of theta potentials type currently 0=glycine
1070 ! VAR:1=non-glicyne non-proline 2=proline
1071 ! VAR:negative values for D-aminoacid
1073 do j=-nthetyp,nthetyp
1074 do k=-nthetyp,nthetyp
1075 read (ithep,'(6a)',end=111,err=111) res1
1076 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
1077 ! VAR: aa0thet is variable describing the average value of Foureir
1078 ! VAR: expansion series
1079 ! VAR: aathet is foureir expansion in theta/2 angle for full formula
1080 ! VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
1081 !ondens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
1082 read (ithep,*,end=111,err=111) &
1083 (aathet(l,i,j,k,iblock),l=1,ntheterm)
1084 read (ithep,*,end=111,err=111) &
1085 ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1086 (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1087 (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1088 (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1090 read (ithep,*,end=111,err=111) &
1091 (((ffthet(llll,lll,ll,i,j,k,iblock),&
1092 ffthet(lll,llll,ll,i,j,k,iblock),&
1093 ggthet(llll,lll,ll,i,j,k,iblock),&
1094 ggthet(lll,llll,ll,i,j,k,iblock),&
1095 llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
1100 ! For dummy ends assign glycine-type coefficients of theta-only terms; the
1101 ! coefficients of theta-and-gamma-dependent terms are zero.
1102 ! IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
1103 ! RECOMENTDED AFTER VERSION 3.3)
1107 ! aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
1108 ! aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
1110 ! aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
1111 ! aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
1114 ! aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
1116 ! aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
1119 ! AND COMMENT THE LOOPS BELOW
1123 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
1124 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
1126 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
1127 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
1130 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1132 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1137 ! Substitution for D aminoacids from symmetry.
1140 do j=-nthetyp,nthetyp
1141 do k=-nthetyp,nthetyp
1142 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
1144 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
1148 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
1149 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
1150 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
1151 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
1157 ffthet(llll,lll,ll,i,j,k,iblock)= &
1158 ffthet(llll,lll,ll,-i,-j,-k,iblock)
1159 ffthet(lll,llll,ll,i,j,k,iblock)= &
1160 ffthet(lll,llll,ll,-i,-j,-k,iblock)
1161 ggthet(llll,lll,ll,i,j,k,iblock)= &
1162 -ggthet(llll,lll,ll,-i,-j,-k,iblock)
1163 ggthet(lll,llll,ll,i,j,k,iblock)= &
1164 -ggthet(lll,llll,ll,-i,-j,-k,iblock)
1173 ! Control printout of the coefficients of virtual-bond-angle potentials
1176 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
1181 write (iout,'(//4a)') &
1182 'Type ',onelett(i),onelett(j),onelett(k)
1183 write (iout,'(//a,10x,a)') " l","a[l]"
1184 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
1185 write (iout,'(i2,1pe15.5)') &
1186 (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
1188 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))') &
1189 "b",l,"c",l,"d",l,"e",l
1191 write (iout,'(i2,4(1pe15.5))') m,&
1192 bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
1193 ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
1197 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))') &
1198 "f+",l,"f-",l,"g+",l,"g-",l
1201 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
1202 ffthet(n,m,l,i,j,k,iblock),&
1203 ffthet(m,n,l,i,j,k,iblock),&
1204 ggthet(n,m,l,i,j,k,iblock),&
1205 ggthet(m,n,l,i,j,k,iblock)
1216 !C here will be the apropriate recalibrating for D-aminoacid
1217 read (ithep,*,end=121,err=121) nthetyp
1218 allocate(nbend_kcc_Tb(-nthetyp:nthetyp))
1219 allocate(v1bend_chyb(0:36,-nthetyp:nthetyp))
1220 do i=-nthetyp+1,nthetyp-1
1221 read (ithep,*,end=121,err=121) nbend_kcc_Tb(i)
1222 do j=0,nbend_kcc_Tb(i)
1223 read (ithep,*,end=121,err=121) ijunk,v1bend_chyb(j,i)
1227 write (iout,'(a)') &
1228 "Parameters of the valence-only potentials"
1229 do i=-nthetyp+1,nthetyp-1
1230 write (iout,'(2a)') "Type ",toronelet(i)
1231 do k=0,nbend_kcc_Tb(i)
1232 write(iout,'(i5,f15.5)') k,v1bend_chyb(k,i)
1238 write (2,*) "Start reading THETA_PDB",ithep_pdb
1240 ! write (2,*) 'i=',i
1241 read (ithep_pdb,*,err=111,end=111) &
1242 a0thet(i),(athet(j,i,1,1),j=1,2),&
1243 (bthet(j,i,1,1),j=1,2)
1244 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
1245 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
1246 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
1247 sigc0(i)=sigc0(i)**2
1250 athet(1,i,1,-1)=athet(1,i,1,1)
1251 athet(2,i,1,-1)=athet(2,i,1,1)
1252 bthet(1,i,1,-1)=-bthet(1,i,1,1)
1253 bthet(2,i,1,-1)=-bthet(2,i,1,1)
1254 athet(1,i,-1,1)=-athet(1,i,1,1)
1255 athet(2,i,-1,1)=-athet(2,i,1,1)
1256 bthet(1,i,-1,1)=bthet(1,i,1,1)
1257 bthet(2,i,-1,1)=bthet(2,i,1,1)
1260 a0thet(i)=a0thet(-i)
1261 athet(1,i,-1,-1)=athet(1,-i,1,1)
1262 athet(2,i,-1,-1)=-athet(2,-i,1,1)
1263 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
1264 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
1265 athet(1,i,-1,1)=athet(1,-i,1,1)
1266 athet(2,i,-1,1)=-athet(2,-i,1,1)
1267 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
1268 bthet(2,i,-1,1)=bthet(2,-i,1,1)
1269 athet(1,i,1,-1)=-athet(1,-i,1,1)
1270 athet(2,i,1,-1)=athet(2,-i,1,1)
1271 bthet(1,i,1,-1)=bthet(1,-i,1,1)
1272 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
1273 theta0(i)=theta0(-i)
1277 polthet(j,i)=polthet(j,-i)
1280 gthet(j,i)=gthet(j,-i)
1283 write (2,*) "End reading THETA_PDB"
1287 !--------------- Reading theta parameters for nucleic acid-------
1288 read (ithep_nucl,*,err=111,end=111) nthetyp_nucl,ntheterm_nucl,&
1289 ntheterm2_nucl,ntheterm3_nucl,nsingle_nucl,ndouble_nucl
1290 nntheterm_nucl=max0(ntheterm_nucl,ntheterm2_nucl,ntheterm3_nucl)
1291 allocate(ithetyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
1292 allocate(aa0thet_nucl(nthetyp_nucl+1,&
1293 nthetyp_nucl+1,nthetyp_nucl+1))
1294 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1295 allocate(aathet_nucl(ntheterm_nucl+1,nthetyp_nucl+1,&
1296 nthetyp_nucl+1,nthetyp_nucl+1))
1297 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1298 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1299 allocate(bbthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1300 nthetyp_nucl+1,nthetyp_nucl+1))
1301 allocate(ccthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1302 nthetyp_nucl+1,nthetyp_nucl+1))
1303 allocate(ddthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1304 nthetyp_nucl+1,nthetyp_nucl+1))
1305 allocate(eethet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1306 nthetyp_nucl+1,nthetyp_nucl+1))
1307 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1308 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1309 allocate(ffthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1310 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1311 allocate(ggthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1312 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1314 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1315 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1317 read (ithep_nucl,*,err=111,end=111) (ithetyp_nucl(i),i=1,ntyp1_molec(2))
1319 aa0thet_nucl(:,:,:)=0.0d0
1320 aathet_nucl(:,:,:,:)=0.0d0
1321 bbthet_nucl(:,:,:,:,:)=0.0d0
1322 ccthet_nucl(:,:,:,:,:)=0.0d0
1323 ddthet_nucl(:,:,:,:,:)=0.0d0
1324 eethet_nucl(:,:,:,:,:)=0.0d0
1325 ffthet_nucl(:,:,:,:,:,:)=0.0d0
1326 ggthet_nucl(:,:,:,:,:,:)=0.0d0
1331 read (ithep_nucl,'(3a)',end=111,err=111) t1,t2,t3
1332 read (ithep_nucl,*,end=111,err=111) aa0thet_nucl(i,j,k)
1333 read (ithep_nucl,*,end=111,err=111)(aathet_nucl(l,i,j,k),l=1,ntheterm_nucl)
1334 read (ithep_nucl,*,end=111,err=111) &
1335 (((bbthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1336 (ccthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1337 (ddthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1338 (eethet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl)),ll=1,ntheterm2_nucl)
1339 read (ithep_nucl,*,end=111,err=111) &
1340 (((ffthet_nucl(llll,lll,ll,i,j,k),ffthet_nucl(lll,llll,ll,i,j,k), &
1341 ggthet_nucl(llll,lll,ll,i,j,k),ggthet_nucl(lll,llll,ll,i,j,k), &
1342 llll=1,lll-1),lll=2,ndouble_nucl),ll=1,ntheterm3_nucl)
1347 !-------------------------------------------
1348 allocate(nlob(ntyp1)) !(ntyp1)
1349 allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
1350 allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
1351 allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
1358 gaussc(:,:,:,:)=0.0D0
1362 ! Read the parameters of the probability distribution/energy expression
1363 ! of the side chains.
1366 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1370 dsc_inv(i)=1.0D0/dsc(i)
1381 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1382 ((blower(k,l,1),l=1,k),k=1,3)
1383 censc(1,1,-i)=censc(1,1,i)
1384 censc(2,1,-i)=censc(2,1,i)
1385 censc(3,1,-i)=-censc(3,1,i)
1387 read (irotam,*,end=112,err=112) bsc(j,i)
1388 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1389 ((blower(k,l,j),l=1,k),k=1,3)
1390 censc(1,j,-i)=censc(1,j,i)
1391 censc(2,j,-i)=censc(2,j,i)
1392 censc(3,j,-i)=-censc(3,j,i)
1393 ! BSC is amplitude of Gaussian
1400 akl=akl+blower(k,m,j)*blower(l,m,j)
1404 if (((k.eq.3).and.(l.ne.3)) &
1405 .or.((l.eq.3).and.(k.ne.3))) then
1406 gaussc(k,l,j,-i)=-akl
1407 gaussc(l,k,j,-i)=-akl
1409 gaussc(k,l,j,-i)=akl
1410 gaussc(l,k,j,-i)=akl
1419 write (iout,'(/a)') 'Parameters of side-chain local geometry'
1422 if (nlobi.gt.0) then
1424 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i,1),&
1425 ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
1426 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
1427 'log h',(bsc(j,i),j=1,nlobi)
1428 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
1429 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
1431 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
1432 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
1435 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
1436 write (iout,'(a,f10.4,4(16x,f10.4))') &
1437 'Center ',(bsc(j,i),j=1,nlobi)
1438 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),&
1447 ! Read scrot parameters for potentials determined from all-atom AM1 calculations
1448 ! added by Urszula Kozlowska 07/11/2007
1450 !el Maximum number of SC local term fitting function coefficiants
1451 !el integer,parameter :: maxsccoef=65
1453 allocate(sc_parmin(65,ntyp)) !(maxsccoef,ntyp)
1456 read (irotam,*,end=112,err=112)
1458 read (irotam,*,end=112,err=112)
1461 read(irotam,*,end=112,err=112) sc_parmin(j,i)
1465 !---------reading nucleic acid parameters for rotamers-------------------
1466 allocate(sc_parmin_nucl(9,ntyp_molec(2))) !(maxsccoef,ntyp)
1467 do i=1,ntyp_molec(2)
1468 read (irotam_nucl,*,end=112,err=112)
1470 read(irotam_nucl,*,end=112,err=112) sc_parmin_nucl(j,i)
1476 write (iout,*) "Base rotamer parameters"
1477 do i=1,ntyp_molec(2)
1478 write (iout,'(a)') restyp(i,2)
1479 write (iout,'(i5,f10.5)') (i,sc_parmin_nucl(j,i),j=1,9)
1484 ! Read the parameters of the probability distribution/energy expression
1485 ! of the side chains.
1487 write (2,*) "Start reading ROTAM_PDB"
1489 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1493 dsc_inv(i)=1.0D0/dsc(i)
1504 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1505 ((blower(k,l,1),l=1,k),k=1,3)
1507 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
1508 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1509 ((blower(k,l,j),l=1,k),k=1,3)
1516 akl=akl+blower(k,m,j)*blower(l,m,j)
1526 write (2,*) "End reading ROTAM_PDB"
1532 !C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
1533 !C interaction energy of the Gly, Ala, and Pro prototypes.
1535 read (ifourier,*) nloctyp
1536 SPLIT_FOURIERTOR = nloctyp.lt.0
1537 nloctyp = iabs(nloctyp)
1538 !C allocate(b1(2,nres)) !(2,-maxtor:maxtor)
1539 !C allocate(b2(2,nres)) !(2,-maxtor:maxtor)
1540 !C allocate(b1tilde(2,nres)) !(2,-maxtor:maxtor)
1541 !C allocate(ctilde(2,2,nres))
1542 !C allocate(dtilde(2,2,nres)) !(2,2,-maxtor:maxtor)
1543 !C allocate(gtb1(2,nres))
1544 !C allocate(gtb2(2,nres))
1545 !C allocate(cc(2,2,nres))
1546 !C allocate(dd(2,2,nres))
1547 !C allocate(ee(2,2,nres))
1548 !C allocate(gtcc(2,2,nres))
1549 !C allocate(gtdd(2,2,nres))
1550 !C allocate(gtee(2,2,nres))
1553 allocate(itype2loc(-ntyp1:ntyp1))
1554 allocate(iloctyp(-nloctyp:nloctyp))
1555 allocate(bnew1(3,2,-nloctyp:nloctyp))
1556 allocate(bnew2(3,2,-nloctyp:nloctyp))
1557 allocate(ccnew(3,2,-nloctyp:nloctyp))
1558 allocate(ddnew(3,2,-nloctyp:nloctyp))
1559 allocate(e0new(3,-nloctyp:nloctyp))
1560 allocate(eenew(2,2,2,-nloctyp:nloctyp))
1561 allocate(bnew1tor(3,2,-nloctyp:nloctyp))
1562 allocate(bnew2tor(3,2,-nloctyp:nloctyp))
1563 allocate(ccnewtor(3,2,-nloctyp:nloctyp))
1564 allocate(ddnewtor(3,2,-nloctyp:nloctyp))
1565 allocate(e0newtor(3,-nloctyp:nloctyp))
1566 allocate(eenewtor(2,2,2,-nloctyp:nloctyp))
1568 read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp)
1569 read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1)
1570 itype2loc(ntyp1)=nloctyp
1571 iloctyp(nloctyp)=ntyp1
1573 itype2loc(-i)=-itype2loc(i)
1582 iloctyp(-i)=-iloctyp(i)
1584 !c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1585 !c write (iout,*) "nloctyp",nloctyp,
1586 !c & " iloctyp",(iloctyp(i),i=0,nloctyp)
1587 !c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1588 !c write (iout,*) "nloctyp",nloctyp,
1589 !c & " iloctyp",(iloctyp(i),i=0,nloctyp)
1592 !c write (iout,*) "NEWCORR",i
1593 read (ifourier,*,end=115,err=115)
1596 read (ifourier,*,end=115,err=115) bnew1(ii,j,i)
1599 !c write (iout,*) "NEWCORR BNEW1"
1600 !c write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2)
1603 read (ifourier,*,end=115,err=115) bnew2(ii,j,i)
1606 !c write (iout,*) "NEWCORR BNEW2"
1607 !c write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2)
1609 read (ifourier,*,end=115,err=115) ccnew(kk,1,i)
1610 read (ifourier,*,end=115,err=115) ccnew(kk,2,i)
1612 !c write (iout,*) "NEWCORR CCNEW"
1613 !c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
1615 read (ifourier,*,end=115,err=115) ddnew(kk,1,i)
1616 read (ifourier,*,end=115,err=115) ddnew(kk,2,i)
1618 !c write (iout,*) "NEWCORR DDNEW"
1619 !c write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2)
1623 read (ifourier,*,end=115,err=115) eenew(ii,jj,kk,i)
1627 !c write (iout,*) "NEWCORR EENEW1"
1628 !c write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
1630 read (ifourier,*,end=115,err=115) e0new(ii,i)
1632 !c write (iout,*) (e0new(ii,i),ii=1,3)
1634 !c write (iout,*) "NEWCORR EENEW"
1637 ccnew(ii,1,i)=ccnew(ii,1,i)/2
1638 ccnew(ii,2,i)=ccnew(ii,2,i)/2
1639 ddnew(ii,1,i)=ddnew(ii,1,i)/2
1640 ddnew(ii,2,i)=ddnew(ii,2,i)/2
1645 bnew1(ii,1,-i)= bnew1(ii,1,i)
1646 bnew1(ii,2,-i)=-bnew1(ii,2,i)
1647 bnew2(ii,1,-i)= bnew2(ii,1,i)
1648 bnew2(ii,2,-i)=-bnew2(ii,2,i)
1651 !c ccnew(ii,1,i)=ccnew(ii,1,i)/2
1652 !c ccnew(ii,2,i)=ccnew(ii,2,i)/2
1653 !c ddnew(ii,1,i)=ddnew(ii,1,i)/2
1654 !c ddnew(ii,2,i)=ddnew(ii,2,i)/2
1655 ccnew(ii,1,-i)=ccnew(ii,1,i)
1656 ccnew(ii,2,-i)=-ccnew(ii,2,i)
1657 ddnew(ii,1,-i)=ddnew(ii,1,i)
1658 ddnew(ii,2,-i)=-ddnew(ii,2,i)
1660 e0new(1,-i)= e0new(1,i)
1661 e0new(2,-i)=-e0new(2,i)
1662 e0new(3,-i)=-e0new(3,i)
1664 eenew(kk,1,1,-i)= eenew(kk,1,1,i)
1665 eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
1666 eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
1667 eenew(kk,2,2,-i)= eenew(kk,2,2,i)
1671 write (iout,'(a)') "Coefficients of the multibody terms"
1672 do i=-nloctyp+1,nloctyp-1
1673 write (iout,*) "Type: ",onelet(iloctyp(i))
1674 write (iout,*) "Coefficients of the expansion of B1"
1676 write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
1678 write (iout,*) "Coefficients of the expansion of B2"
1680 write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
1682 write (iout,*) "Coefficients of the expansion of C"
1683 write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
1684 write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
1685 write (iout,*) "Coefficients of the expansion of D"
1686 write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
1687 write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
1688 write (iout,*) "Coefficients of the expansion of E"
1689 write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
1692 write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
1697 IF (SPLIT_FOURIERTOR) THEN
1699 !c write (iout,*) "NEWCORR TOR",i
1700 read (ifourier,*,end=115,err=115)
1703 read (ifourier,*,end=115,err=115) bnew1tor(ii,j,i)
1706 !c write (iout,*) "NEWCORR BNEW1 TOR"
1707 !c write (iout,*) ((bnew1tor(ii,j,i),ii=1,3),j=1,2)
1710 read (ifourier,*,end=115,err=115) bnew2tor(ii,j,i)
1713 !c write (iout,*) "NEWCORR BNEW2 TOR"
1714 !c write (iout,*) ((bnew2tor(ii,j,i),ii=1,3),j=1,2)
1716 read (ifourier,*,end=115,err=115) ccnewtor(kk,1,i)
1717 read (ifourier,*,end=115,err=115) ccnewtor(kk,2,i)
1719 !c write (iout,*) "NEWCORR CCNEW TOR"
1720 !c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
1722 read (ifourier,*,end=115,err=115) ddnewtor(kk,1,i)
1723 read (ifourier,*,end=115,err=115) ddnewtor(kk,2,i)
1725 !c write (iout,*) "NEWCORR DDNEW TOR"
1726 !c write (iout,*) ((ddnewtor(ii,j,i),ii=1,3),j=1,2)
1730 read (ifourier,*,end=115,err=115) eenewtor(ii,jj,kk,i)
1734 !c write (iout,*) "NEWCORR EENEW1 TOR"
1735 !c write(iout,*)(((eenewtor(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
1737 read (ifourier,*,end=115,err=115) e0newtor(ii,i)
1739 !c write (iout,*) (e0newtor(ii,i),ii=1,3)
1741 !c write (iout,*) "NEWCORR EENEW TOR"
1744 ccnewtor(ii,1,i)=ccnewtor(ii,1,i)/2
1745 ccnewtor(ii,2,i)=ccnewtor(ii,2,i)/2
1746 ddnewtor(ii,1,i)=ddnewtor(ii,1,i)/2
1747 ddnewtor(ii,2,i)=ddnewtor(ii,2,i)/2
1752 bnew1tor(ii,1,-i)= bnew1tor(ii,1,i)
1753 bnew1tor(ii,2,-i)=-bnew1tor(ii,2,i)
1754 bnew2tor(ii,1,-i)= bnew2tor(ii,1,i)
1755 bnew2tor(ii,2,-i)=-bnew2tor(ii,2,i)
1758 ccnewtor(ii,1,-i)=ccnewtor(ii,1,i)
1759 ccnewtor(ii,2,-i)=-ccnewtor(ii,2,i)
1760 ddnewtor(ii,1,-i)=ddnewtor(ii,1,i)
1761 ddnewtor(ii,2,-i)=-ddnewtor(ii,2,i)
1763 e0newtor(1,-i)= e0newtor(1,i)
1764 e0newtor(2,-i)=-e0newtor(2,i)
1765 e0newtor(3,-i)=-e0newtor(3,i)
1767 eenewtor(kk,1,1,-i)= eenewtor(kk,1,1,i)
1768 eenewtor(kk,1,2,-i)=-eenewtor(kk,1,2,i)
1769 eenewtor(kk,2,1,-i)=-eenewtor(kk,2,1,i)
1770 eenewtor(kk,2,2,-i)= eenewtor(kk,2,2,i)
1774 write (iout,'(a)') &
1775 "Single-body coefficients of the torsional potentials"
1776 do i=-nloctyp+1,nloctyp-1
1777 write (iout,*) "Type: ",onelet(iloctyp(i))
1778 write (iout,*) "Coefficients of the expansion of B1tor"
1780 write (iout,'(3hB1(,i1,1h),3f10.5)') &
1781 j,(bnew1tor(k,j,i),k=1,3)
1783 write (iout,*) "Coefficients of the expansion of B2tor"
1785 write (iout,'(3hB2(,i1,1h),3f10.5)') &
1786 j,(bnew2tor(k,j,i),k=1,3)
1788 write (iout,*) "Coefficients of the expansion of Ctor"
1789 write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3)
1790 write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3)
1791 write (iout,*) "Coefficients of the expansion of Dtor"
1792 write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3)
1793 write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3)
1794 write (iout,*) "Coefficients of the expansion of Etor"
1795 write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3)
1798 write (iout,'(1hE,2i1,2f10.5)') &
1799 j,k,(eenewtor(l,j,k,i),l=1,2)
1805 do i=-nloctyp+1,nloctyp-1
1808 bnew1tor(ii,j,i)=bnew1(ii,j,i)
1813 bnew2tor(ii,j,i)=bnew2(ii,j,i)
1817 ccnewtor(ii,1,i)=ccnew(ii,1,i)
1818 ccnewtor(ii,2,i)=ccnew(ii,2,i)
1819 ddnewtor(ii,1,i)=ddnew(ii,1,i)
1820 ddnewtor(ii,2,i)=ddnew(ii,2,i)
1823 ENDIF !SPLIT_FOURIER_TOR
1825 allocate(ccold(2,2,-nloctyp-1:nloctyp+1))
1826 allocate(ddold(2,2,-nloctyp-1:nloctyp+1))
1827 allocate(eeold(2,2,-nloctyp-1:nloctyp+1))
1830 write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)"
1832 read (ifourier,*,end=115,err=115)
1833 read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
1835 write (iout,*) 'Type ',onelet(iloctyp(i))
1836 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
1850 c B1tilde(1,i) = b(3)
1851 c B1tilde(2,i) =-b(5)
1852 c B1tilde(1,-i) =-b(3)
1853 c B1tilde(2,-i) =b(5)
1854 c b1tilde(1,i)=0.0d0
1855 c b1tilde(2,i)=0.0d0
1860 cc B1tilde(1,i) = b(3,i)
1861 cc B1tilde(2,i) =-b(5,i)
1862 C B1tilde(1,-i) =-b(3,i)
1863 C B1tilde(2,-i) =b(5,i)
1864 cc b1tilde(1,i)=0.0d0
1865 cc b1tilde(2,i)=0.0d0
1873 CCold(1,1,i)= b(7,i)
1874 CCold(2,2,i)=-b(7,i)
1875 CCold(2,1,i)= b(9,i)
1876 CCold(1,2,i)= b(9,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)
1885 c Ctilde(1,1,i)= CCold(1,1,i)
1886 c Ctilde(1,2,i)= CCold(1,2,i)
1887 c Ctilde(2,1,i)=-CCold(2,1,i)
1888 c Ctilde(2,2,i)=-CCold(2,2,i)
1893 c Ctilde(1,1,i)= CCold(1,1,i)
1894 c Ctilde(1,2,i)= CCold(1,2,i)
1895 c Ctilde(2,1,i)=-CCold(2,1,i)
1896 c Ctilde(2,2,i)=-CCold(2,2,i)
1898 c Ctilde(1,1,i)=0.0d0
1899 c Ctilde(1,2,i)=0.0d0
1900 c Ctilde(2,1,i)=0.0d0
1901 c Ctilde(2,2,i)=0.0d0
1902 DDold(1,1,i)= b(6,i)
1903 DDold(2,2,i)=-b(6,i)
1904 DDold(2,1,i)= b(8,i)
1905 DDold(1,2,i)= b(8,i)
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)
1914 c Dtilde(1,1,i)= DD(1,1,i)
1915 c Dtilde(1,2,i)= DD(1,2,i)
1916 c Dtilde(2,1,i)=-DD(2,1,i)
1917 c Dtilde(2,2,i)=-DD(2,2,i)
1919 c Dtilde(1,1,i)=0.0d0
1920 c Dtilde(1,2,i)=0.0d0
1921 c Dtilde(2,1,i)=0.0d0
1922 c Dtilde(2,2,i)=0.0d0
1923 EEold(1,1,i)= b(10,i)+b(11,i)
1924 EEold(2,2,i)=-b(10,i)+b(11,i)
1925 EEold(2,1,i)= b(12,i)-b(13,i)
1926 EEold(1,2,i)= b(12,i)+b(13,i)
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 write(iout,*) "TU DOCHODZE"
1937 c ee(2,1,i)=ee(1,2,i)
1942 "Coefficients of the cumulants (independent of valence angles)"
1943 do i=-nloctyp+1,nloctyp-1
1944 write (iout,*) 'Type ',onelet(iloctyp(i))
1946 write(iout,'(2f10.5)') B(3,i),B(5,i)
1948 write(iout,'(2f10.5)') B(2,i),B(4,i)
1951 write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
1955 write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
1959 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1968 ! Read torsional parameters in old format
1970 allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
1972 read (itorp,*,end=113,err=113) ntortyp,nterm_old
1973 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1974 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1976 !el from energy module--------
1977 allocate(v1(nterm_old,ntortyp,ntortyp))
1978 allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
1979 !el---------------------------
1984 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
1990 write (iout,'(/a/)') 'Torsional constants:'
1993 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1994 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
2000 ! Read torsional parameters
2002 IF (TOR_MODE.eq.0) THEN
2003 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
2004 read (itorp,*,end=113,err=113) ntortyp
2005 !el from energy module---------
2006 allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2007 allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2009 allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
2010 allocate(vlor2(maxlor,ntortyp,ntortyp))
2011 allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
2012 allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2014 allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
2015 allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
2016 !el---------------------------
2019 !el---------------------------
2021 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
2023 itortyp(i)=-itortyp(-i)
2025 itortyp(ntyp1)=ntortyp
2026 itortyp(-ntyp1)=-ntortyp
2028 write (iout,*) 'ntortyp',ntortyp
2030 do j=-ntortyp+1,ntortyp-1
2031 read (itorp,*,end=113,err=113) nterm(i,j,iblock),&
2033 nterm(-i,-j,iblock)=nterm(i,j,iblock)
2034 nlor(-i,-j,iblock)=nlor(i,j,iblock)
2037 do k=1,nterm(i,j,iblock)
2038 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),&
2040 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
2041 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
2042 v0ij=v0ij+si*v1(k,i,j,iblock)
2044 ! write(iout,*) i,j,k,iblock,nterm(i,j,iblock) !
2045 ! write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),&!
2046 ! v2(k,-i,-j,iblock),v2(k,i,j,iblock)!
2048 do k=1,nlor(i,j,iblock)
2049 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),&
2050 vlor2(k,i,j),vlor3(k,i,j)
2051 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
2054 v0(-i,-j,iblock)=v0ij
2060 write (iout,'(/a/)') 'Torsional constants:'
2062 do i=-ntortyp,ntortyp
2063 do j=-ntortyp,ntortyp
2064 write (iout,*) 'ityp',i,' jtyp',j
2065 write (iout,*) 'Fourier constants'
2066 do k=1,nterm(i,j,iblock)
2067 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
2070 write (iout,*) 'Lorenz constants'
2071 do k=1,nlor(i,j,iblock)
2072 write (iout,'(3(1pe15.5))') &
2073 vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
2079 !elwrite (iout,'(/a/)') 'Torsional constants:',vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
2081 ! 6/23/01 Read parameters for double torsionals
2083 !el from energy module------------
2084 allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2085 allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2086 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
2087 allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2088 allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2089 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
2090 allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2091 allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2092 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
2093 !---------------------------------
2097 do j=-ntortyp+1,ntortyp-1
2098 do k=-ntortyp+1,ntortyp-1
2099 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
2100 ! write (iout,*) "OK onelett",
2103 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
2104 .or. t3.ne.toronelet(k)) then
2105 write (iout,*) "Error in double torsional parameter file",&
2108 call MPI_Finalize(Ierror)
2110 stop "Error in double torsional parameter file"
2112 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),&
2113 ntermd_2(i,j,k,iblock)
2114 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
2115 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
2116 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,&
2117 ntermd_1(i,j,k,iblock))
2118 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,&
2119 ntermd_1(i,j,k,iblock))
2120 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,&
2121 ntermd_1(i,j,k,iblock))
2122 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,&
2123 ntermd_1(i,j,k,iblock))
2124 ! Martix of D parameters for one dimesional foureir series
2125 do l=1,ntermd_1(i,j,k,iblock)
2126 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
2127 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
2128 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
2129 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
2130 ! write(iout,*) "whcodze" ,
2131 ! & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
2133 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),&
2134 v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
2135 v2s(m,l,i,j,k,iblock),&
2136 m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
2137 ! Martix of D parameters for two dimesional fourier series
2138 do l=1,ntermd_2(i,j,k,iblock)
2140 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
2141 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
2142 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
2143 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
2152 write (iout,*) 'Constants for double torsionals'
2155 do j=-ntortyp+1,ntortyp-1
2156 do k=-ntortyp+1,ntortyp-1
2157 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
2158 ' nsingle',ntermd_1(i,j,k,iblock),&
2159 ' ndouble',ntermd_2(i,j,k,iblock)
2161 write (iout,*) 'Single angles:'
2162 do l=1,ntermd_1(i,j,k,iblock)
2163 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
2164 v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
2165 v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
2166 v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
2169 write (iout,*) 'Pairs of angles:'
2170 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
2171 do l=1,ntermd_2(i,j,k,iblock)
2172 write (iout,'(i5,20f10.5)') &
2173 l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
2176 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
2177 do l=1,ntermd_2(i,j,k,iblock)
2178 write (iout,'(i5,20f10.5)') &
2179 l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
2180 (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
2188 ELSE IF (TOR_MODE.eq.1) THEN
2190 !C read valence-torsional parameters
2191 read (itorp,*,end=121,err=121) ntortyp
2193 write (iout,*) "Valence-torsional parameters read in ntortyp",&
2195 read (itorp,*,end=121,err=121) (itortyp(i),i=1,ntyp)
2196 write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp)
2199 itype2loc(i)=itortyp(i)
2203 itortyp(i)=-itortyp(-i)
2205 do i=-ntortyp+1,ntortyp-1
2206 do j=-ntortyp+1,ntortyp-1
2207 !C first we read the cos and sin gamma parameters
2208 read (itorp,'(13x,a)',end=121,err=121) string
2209 write (iout,*) i,j,string
2210 read (itorp,*,end=121,err=121) &
2211 nterm_kcc(j,i),nterm_kcc_Tb(j,i)
2212 !C read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
2213 do k=1,nterm_kcc(j,i)
2214 do l=1,nterm_kcc_Tb(j,i)
2215 do ll=1,nterm_kcc_Tb(j,i)
2216 read (itorp,*,end=121,err=121) ii,jj,kk, &
2217 v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
2225 !c AL 4/8/16: Calculate coefficients from one-body parameters
2227 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
2228 allocate(nterm_kcc(-ntyp1:ntyp1,-ntyp1:ntyp1))
2229 allocate(nterm_kcc_Tb(-ntyp1:ntyp1,-ntyp1:ntyp1))
2230 allocate(v1_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1))
2231 allocate(v2_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1))
2234 print *,i,itortyp(i)
2235 itortyp(i)=itype2loc(i)
2238 "Val-tor parameters calculated from cumulant coefficients ntortyp"&
2240 do i=-ntortyp+1,ntortyp-1
2241 do j=-ntortyp+1,ntortyp-1
2244 do k=1,nterm_kcc_Tb(j,i)
2245 do l=1,nterm_kcc_Tb(j,i)
2246 v1_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,1,j)&
2247 +bnew1tor(k,2,i)*bnew2tor(l,2,j)
2248 v2_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,2,j)&
2249 +bnew1tor(k,2,i)*bnew2tor(l,1,j)
2252 do k=1,nterm_kcc_Tb(j,i)
2253 do l=1,nterm_kcc_Tb(j,i)
2255 v1_kcc(k,l,2,i,j)=-(ccnewtor(k,1,i)*ddnewtor(l,1,j) &
2256 -ccnewtor(k,2,i)*ddnewtor(l,2,j))
2257 v2_kcc(k,l,2,i,j)=-(ccnewtor(k,2,i)*ddnewtor(l,1,j) &
2258 +ccnewtor(k,1,i)*ddnewtor(l,2,j))
2260 v1_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,1,i)*ddnewtor(l,1,j) &
2261 -ccnewtor(k,2,i)*ddnewtor(l,2,j))
2262 v2_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,2,i)*ddnewtor(l,1,j) &
2263 +ccnewtor(k,1,i)*ddnewtor(l,2,j))
2267 !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)
2271 write (iout,*) "TOR_MODE>1 only with NEWCORR"
2276 if (tor_mode.gt.0 .and. lprint) then
2277 !c Print valence-torsional parameters
2278 write (iout,'(a)') &
2279 "Parameters of the valence-torsional potentials"
2280 do i=-ntortyp+1,ntortyp-1
2281 do j=-ntortyp+1,ntortyp-1
2282 write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
2283 write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
2284 do k=1,nterm_kcc(j,i)
2285 do l=1,nterm_kcc_Tb(j,i)
2286 do ll=1,nterm_kcc_Tb(j,i)
2287 write (iout,'(3i5,2f15.4)')&
2288 k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
2296 allocate(itortyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
2297 read (itorp_nucl,*,end=113,err=113) ntortyp_nucl
2298 ! print *,"ntortyp_nucl",ntortyp_nucl,ntyp_molec(2)
2299 !el from energy module---------
2300 allocate(nterm_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2301 allocate(nlor_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2303 allocate(vlor1_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
2304 allocate(vlor2_nucl(maxlor,ntortyp_nucl,ntortyp_nucl))
2305 allocate(vlor3_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,maxtor,maxtor)
2306 allocate(v0_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2308 allocate(v1_nucl(maxterm,ntortyp_nucl,ntortyp_nucl))
2309 allocate(v2_nucl(maxterm,ntortyp_nucl,ntortyp_nucl)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
2310 !el---------------------------
2313 !el--------------------
2314 read (itorp_nucl,*,end=113,err=113) &
2315 (itortyp_nucl(i),i=1,ntyp_molec(2))
2316 ! print *,itortyp_nucl(:)
2317 !c write (iout,*) 'ntortyp',ntortyp
2320 read (itorp_nucl,*,end=113,err=113) nterm_nucl(i,j),nlor_nucl(i,j)
2321 ! print *,nterm_nucl(i,j),nlor_nucl(i,j)
2324 do k=1,nterm_nucl(i,j)
2325 read (itorp_nucl,*,end=113,err=113) kk,v1_nucl(k,i,j),v2_nucl(k,i,j)
2326 v0ij=v0ij+si*v1_nucl(k,i,j)
2329 do k=1,nlor_nucl(i,j)
2330 read (itorp,*,end=113,err=113) kk,vlor1_nucl(k,i,j),&
2331 vlor2_nucl(k,i,j),vlor3_nucl(k,i,j)
2332 v0ij=v0ij+vlor1_nucl(k,i,j)/(1+vlor3_nucl(k,i,j)**2)
2338 ! Read of Side-chain backbone correlation parameters
2339 ! Modified 11 May 2012 by Adasko
2342 read (isccor,*,end=119,err=119) nsccortyp
2344 !el from module energy-------------
2345 allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2346 allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
2347 allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
2348 allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp)) !(maxterm_sccor,20,20)
2349 !-----------------------------------
2351 !el from module energy-------------
2352 allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
2354 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
2356 isccortyp(i)=-isccortyp(-i)
2358 iscprol=isccortyp(20)
2359 ! write (iout,*) 'ntortyp',ntortyp
2361 !c maxinter is maximum interaction sites
2362 !el from module energy---------
2363 allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2364 allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
2365 -nsccortyp:nsccortyp))
2366 allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
2367 -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
2368 allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
2369 -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
2370 !-----------------------------------
2374 read (isccor,*,end=119,err=119) &
2375 nterm_sccor(i,j),nlor_sccor(i,j)
2381 nterm_sccor(-i,j)=nterm_sccor(i,j)
2382 nterm_sccor(-i,-j)=nterm_sccor(i,j)
2383 nterm_sccor(i,-j)=nterm_sccor(i,j)
2384 do k=1,nterm_sccor(i,j)
2385 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
2387 if (j.eq.iscprol) then
2388 if (i.eq.isccortyp(10)) then
2389 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
2390 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2392 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
2393 +v2sccor(k,l,i,j)*dsqrt(0.75d0)
2394 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
2395 +v1sccor(k,l,i,j)*dsqrt(0.75d0)
2396 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
2397 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
2398 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
2399 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
2402 if (i.eq.isccortyp(10)) then
2403 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
2404 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2406 if (j.eq.isccortyp(10)) then
2407 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
2408 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
2410 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
2411 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2412 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
2413 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
2414 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
2415 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
2419 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
2420 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
2421 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
2422 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
2425 do k=1,nlor_sccor(i,j)
2426 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
2427 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2428 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
2429 (1+vlor3sccor(k,i,j)**2)
2431 v0sccor(l,i,j)=v0ijsccor
2432 v0sccor(l,-i,j)=v0ijsccor1
2433 v0sccor(l,i,-j)=v0ijsccor2
2434 v0sccor(l,-i,-j)=v0ijsccor3
2440 !el from module energy-------------
2441 allocate(isccortyp(ntyp)) !(-ntyp:ntyp)
2443 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
2444 ! write (iout,*) 'ntortyp',ntortyp
2446 !c maxinter is maximum interaction sites
2447 !el from module energy---------
2448 allocate(nterm_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2449 allocate(v1sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp))
2450 allocate(v2sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
2451 allocate(v0sccor(maxinter,nsccortyp,nsccortyp)) !???(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
2452 !-----------------------------------
2456 read (isccor,*,end=119,err=119) &
2457 nterm_sccor(i,j),nlor_sccor(i,j)
2461 do k=1,nterm_sccor(i,j)
2462 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
2464 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
2467 do k=1,nlor_sccor(i,j)
2468 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
2469 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2470 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
2471 (1+vlor3sccor(k,i,j)**2)
2473 v0sccor(l,i,j)=v0ijsccor !el ,iblock
2482 write (iout,'(/a/)') 'Torsional constants:'
2485 write (iout,*) 'ityp',i,' jtyp',j
2486 write (iout,*) 'Fourier constants'
2487 do k=1,nterm_sccor(i,j)
2488 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
2490 write (iout,*) 'Lorenz constants'
2491 do k=1,nlor_sccor(i,j)
2492 write (iout,'(3(1pe15.5))') &
2493 vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2500 ! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
2501 ! interaction energy of the Gly, Ala, and Pro prototypes.
2504 ! Read electrostatic-interaction parameters
2509 write (iout,'(/a)') 'Electrostatic interaction constants:'
2510 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
2511 'IT','JT','APP','BPP','AEL6','AEL3'
2513 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
2514 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
2515 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
2516 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
2521 app (i,j)=epp(i,j)*rri*rri
2522 bpp (i,j)=-2.0D0*epp(i,j)*rri
2523 ael6(i,j)=elpp6(i,j)*4.2D0**6
2524 ael3(i,j)=elpp3(i,j)*4.2D0**3
2526 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
2532 ! Read side-chain interaction parameters.
2534 !el from module energy - COMMON.INTERACT-------
2535 allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
2536 allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
2537 allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
2539 allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
2540 allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
2541 allocate(epslip(ntyp,ntyp))
2549 !--------------------------------
2551 read (isidep,*,end=117,err=117) ipot,expon
2552 if (ipot.lt.1 .or. ipot.gt.5) then
2553 write (iout,'(2a)') 'Error while reading SC interaction',&
2554 'potential file - unknown potential type.'
2556 call MPI_Finalize(Ierror)
2562 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
2563 ', exponents are ',expon,2*expon
2564 ! goto (10,20,30,30,40) ipot
2566 !----------------------- LJ potential ---------------------------------
2568 ! 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2569 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2570 (sigma0(i),i=1,ntyp)
2572 write (iout,'(/a/)') 'Parameters of the LJ potential:'
2573 write (iout,'(a/)') 'The epsilon array:'
2574 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2575 write (iout,'(/a)') 'One-body parameters:'
2576 write (iout,'(a,4x,a)') 'residue','sigma'
2577 write (iout,'(a3,6x,f10.5)') (restyp(i,1),sigma0(i),i=1,ntyp)
2580 !----------------------- LJK potential --------------------------------
2582 ! 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2583 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2584 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
2586 write (iout,'(/a/)') 'Parameters of the LJK potential:'
2587 write (iout,'(a/)') 'The epsilon array:'
2588 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2589 write (iout,'(/a)') 'One-body parameters:'
2590 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
2591 write (iout,'(a3,6x,2f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2595 !---------------------- GB or BP potential -----------------------------
2598 ! print *,"I AM in SCELE",scelemode
2599 if (scelemode.eq.0) then
2601 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
2603 read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
2604 read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
2605 read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
2606 read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
2608 read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
2611 ! For the GB potential convert sigma'**2 into chi'
2614 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
2618 write (iout,'(/a/)') 'Parameters of the BP potential:'
2619 write (iout,'(a/)') 'The epsilon array:'
2620 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2621 write (iout,'(/a)') 'One-body parameters:'
2622 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',&
2624 write (iout,'(a3,6x,4f10.5)') (restyp(i,1),sigma0(i),sigii(i),&
2625 chip(i),alp(i),i=1,ntyp)
2628 ! print *,ntyp,"NTYP"
2629 allocate(icharge(ntyp1))
2630 ! print *,ntyp,icharge(i)
2632 read (isidep,*) (icharge(i),i=1,ntyp)
2633 print *,ntyp,icharge(i)
2634 ! if(.not.allocated(eps)) allocate(eps(-ntyp
2635 !c write (2,*) "icharge",(icharge(i),i=1,ntyp)
2636 allocate(alphapol(ntyp,ntyp),epshead(ntyp,ntyp),sig0head(ntyp,ntyp))
2637 allocate(sigiso1(ntyp,ntyp),rborn(ntyp,ntyp),sigmap1(ntyp,ntyp))
2638 allocate(sigmap2(ntyp,ntyp),sigiso2(ntyp,ntyp))
2639 allocate(chis(ntyp,ntyp),wquad(ntyp,ntyp),chipp(ntyp,ntyp))
2640 allocate(epsintab(ntyp,ntyp))
2641 allocate(dtail(2,ntyp,ntyp))
2642 allocate(alphasur(4,ntyp,ntyp),alphiso(4,ntyp,ntyp))
2643 allocate(wqdip(2,ntyp,ntyp))
2644 allocate(wstate(4,ntyp,ntyp))
2645 allocate(dhead(2,2,ntyp,ntyp))
2646 allocate(nstate(ntyp,ntyp))
2647 if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2648 if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2651 ! write (*,*) "Im in ALAB", i, " ", j
2653 eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), &
2654 (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), &
2655 chis(i,j),chis(j,i), &
2656 nstate(i,j),(wstate(k,i,j),k=1,4), &
2657 dhead(1,1,i,j),dhead(1,2,i,j),dhead(2,1,i,j),dhead(2,2,i,j),&
2658 dtail(1,i,j),dtail(2,i,j), &
2659 epshead(i,j),sig0head(i,j), &
2660 rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j), &
2661 alphapol(i,j),alphapol(j,i), &
2662 (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j)
2663 ! print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i)
2669 sigma(i,j) = sigma(j,i)
2670 sigmap1(i,j)=sigmap1(j,i)
2671 sigmap2(i,j)=sigmap2(j,i)
2672 sigiso1(i,j)=sigiso1(j,i)
2673 sigiso2(i,j)=sigiso2(j,i)
2674 ! print *,"ATU",sigma(j,i),sigma(i,j),i,j
2675 nstate(i,j) = nstate(j,i)
2676 dtail(1,i,j) = dtail(1,j,i)
2677 dtail(2,i,j) = dtail(2,j,i)
2679 alphasur(k,i,j) = alphasur(k,j,i)
2680 wstate(k,i,j) = wstate(k,j,i)
2681 alphiso(k,i,j) = alphiso(k,j,i)
2684 dhead(2,1,i,j) = dhead(1,1,j,i)
2685 dhead(2,2,i,j) = dhead(1,2,j,i)
2686 dhead(1,1,i,j) = dhead(2,1,j,i)
2687 dhead(1,2,i,j) = dhead(2,2,j,i)
2689 epshead(i,j) = epshead(j,i)
2690 sig0head(i,j) = sig0head(j,i)
2693 wqdip(k,i,j) = wqdip(k,j,i)
2696 wquad(i,j) = wquad(j,i)
2697 epsintab(i,j) = epsintab(j,i)
2698 ! if (epsintab(i,j).ne.1.0) print *,"WHAT?",i,j,epsintab(i,j)
2703 !--------------------- GBV potential -----------------------------------
2705 ! 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2706 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2707 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
2708 (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
2710 write (iout,'(/a/)') 'Parameters of the GBV potential:'
2711 write (iout,'(a/)') 'The epsilon array:'
2712 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2713 write (iout,'(/a)') 'One-body parameters:'
2714 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',&
2715 's||/s_|_^2',' chip ',' alph '
2716 write (iout,'(a3,6x,5f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2717 sigii(i),chip(i),alp(i),i=1,ntyp)
2720 write(iout,*)"Wrong ipot"
2726 !-----------------------------------------------------------------------
2727 ! Calculate the "working" parameters of SC interactions.
2729 !el from module energy - COMMON.INTERACT-------
2730 allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1))
2731 if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2732 allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)
2733 if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2734 allocate(r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
2735 allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),&
2737 allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),&
2743 if (scelemode.eq.0) then
2752 sc_aa_tube_par(:)=0.0d0
2753 sc_bb_tube_par(:)=0.0d0
2755 !--------------------------------
2760 epslip(i,j)=epslip(j,i)
2763 if (scelemode.eq.0) then
2766 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
2767 sigma(j,i)=sigma(i,j)
2768 rs0(i,j)=dwa16*sigma(i,j)
2773 if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
2774 'Working parameters of the SC interactions:',&
2775 ' a ',' b ',' augm ',' sigma ',' r0 ',&
2780 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2782 ! print *,"SIGMA ZLA?",sigma(i,j)
2790 sigeps=dsign(1.0D0,epsij)
2792 aa_aq(i,j)=epsij*rrij*rrij
2793 print *,"ADASKO",epsij,rrij,expon
2794 bb_aq(i,j)=-sigeps*epsij*rrij
2795 aa_aq(j,i)=aa_aq(i,j)
2796 bb_aq(j,i)=bb_aq(i,j)
2797 epsijlip=epslip(i,j)
2798 sigeps=dsign(1.0D0,epsijlip)
2799 epsijlip=dabs(epsijlip)
2800 aa_lip(i,j)=epsijlip*rrij*rrij
2801 bb_lip(i,j)=-sigeps*epsijlip*rrij
2802 aa_lip(j,i)=aa_lip(i,j)
2803 bb_lip(j,i)=bb_lip(i,j)
2804 !C write(iout,*) aa_lip
2805 if ((ipot.gt.2).and. (scelemode.eq.0)) then
2806 sigt1sq=sigma0(i)**2
2807 sigt2sq=sigma0(j)**2
2810 ratsig1=sigt2sq/sigt1sq
2811 ratsig2=1.0D0/ratsig1
2812 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
2813 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
2814 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
2818 ! if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2819 sigmaii(i,j)=rsum_max
2820 sigmaii(j,i)=rsum_max
2822 ! sigmaii(i,j)=r0(i,j)
2823 ! sigmaii(j,i)=r0(i,j)
2825 !d write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
2826 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
2827 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
2828 augm(i,j)=epsij*r_augm**(2*expon)
2829 ! augm(i,j)=0.5D0**(2*expon)*aa(i,j)
2836 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
2837 restyp(i,1),restyp(j,1),aa_aq(i,j),bb_aq(i,j),augm(i,j),&
2838 sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
2843 allocate(eps_nucl(ntyp_molec(2),ntyp_molec(2)))
2844 allocate(sigma_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2845 allocate(elpp6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2846 allocate(elpp3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2847 allocate(elpp63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2848 allocate(elpp32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2849 allocate(chi_nucl(ntyp_molec(2),ntyp_molec(2)),chip_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2850 allocate(ael3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2851 allocate(ael6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2852 allocate(ael32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2853 allocate(ael63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2854 allocate(aa_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2855 allocate(bb_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2856 allocate(r0_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2857 allocate(sigmaii_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2858 allocate(eps_scp_nucl(ntyp_molec(2)),rscp_nucl(ntyp_molec(2))) !(ntyp,2)
2867 read (isidep_nucl,*) ipot_nucl
2868 ! print *,"TU?!",ipot_nucl
2869 if (ipot_nucl.eq.1) then
2870 do i=1,ntyp_molec(2)
2871 do j=i,ntyp_molec(2)
2872 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),elpp6_nucl(i,j),&
2873 elpp3_nucl(i,j), elpp63_nucl(i,j),elpp32_nucl(i,j)
2877 do i=1,ntyp_molec(2)
2878 do j=i,ntyp_molec(2)
2879 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),chi_nucl(i,j),&
2880 chi_nucl(j,i),chip_nucl(i,j),chip_nucl(j,i),&
2881 elpp6_nucl(i,j),elpp3_nucl(i,j),elpp63_nucl(i,j),elpp32_nucl(i,j)
2885 ! rpp(1,1)=2**(1.0/6.0)*5.16158
2886 do i=1,ntyp_molec(2)
2887 do j=i,ntyp_molec(2)
2888 rrij=sigma_nucl(i,j)
2892 epsij=4*eps_nucl(i,j)
2893 sigeps=dsign(1.0D0,epsij)
2895 aa_nucl(i,j)=epsij*rrij*rrij
2896 bb_nucl(i,j)=-sigeps*epsij*rrij
2897 ael3_nucl(i,j)=elpp3_nucl(i,j)*dsqrt(rrij)
2898 ael6_nucl(i,j)=elpp6_nucl(i,j)*rrij
2899 ael63_nucl(i,j)=elpp63_nucl(i,j)*rrij
2900 ael32_nucl(i,j)=elpp32_nucl(i,j)*rrij
2901 sigmaii_nucl(i,j)=sigma_nucl(i,j)/sqrt(1-(chi_nucl(i,j)+chi_nucl(j,i)- &
2902 2*chi_nucl(i,j)*chi_nucl(j,i))/(1-chi_nucl(i,j)*chi_nucl(j,i)))
2905 aa_nucl(i,j)=aa_nucl(j,i)
2906 bb_nucl(i,j)=bb_nucl(j,i)
2907 ael3_nucl(i,j)=ael3_nucl(j,i)
2908 ael6_nucl(i,j)=ael6_nucl(j,i)
2909 ael63_nucl(i,j)=ael63_nucl(j,i)
2910 ael32_nucl(i,j)=ael32_nucl(j,i)
2911 elpp3_nucl(i,j)=elpp3_nucl(j,i)
2912 elpp6_nucl(i,j)=elpp6_nucl(j,i)
2913 elpp63_nucl(i,j)=elpp63_nucl(j,i)
2914 elpp32_nucl(i,j)=elpp32_nucl(j,i)
2915 eps_nucl(i,j)=eps_nucl(j,i)
2916 sigma_nucl(i,j)=sigma_nucl(j,i)
2917 sigmaii_nucl(i,j)=sigmaii_nucl(j,i)
2921 write(iout,*) "tube param"
2922 read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep, &
2923 ccavtubpep,dcavtubpep,tubetranenepep
2924 sigmapeptube=sigmapeptube**6
2925 sigeps=dsign(1.0D0,epspeptube)
2926 epspeptube=dabs(epspeptube)
2927 pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
2928 pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
2929 write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
2931 read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i), &
2932 ccavtub(i),dcavtub(i),tubetranene(i)
2933 sigmasctube=sigmasctube**6
2934 sigeps=dsign(1.0D0,epssctube)
2935 epssctube=dabs(epssctube)
2936 sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
2937 sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
2938 write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
2940 !-----------------READING SC BASE POTENTIALS-----------------------------
2941 allocate(eps_scbase(ntyp_molec(1),ntyp_molec(2)))
2942 allocate(sigma_scbase(ntyp_molec(1),ntyp_molec(2)))
2943 allocate(chi_scbase(ntyp_molec(1),ntyp_molec(2),2))
2944 allocate(chipp_scbase(ntyp_molec(1),ntyp_molec(2),2))
2945 allocate(alphasur_scbase(4,ntyp_molec(1),ntyp_molec(2)))
2946 allocate(sigmap1_scbase(ntyp_molec(1),ntyp_molec(2)))
2947 allocate(sigmap2_scbase(ntyp_molec(1),ntyp_molec(2)))
2948 allocate(chis_scbase(ntyp_molec(1),ntyp_molec(2),2))
2949 allocate(dhead_scbasei(ntyp_molec(1),ntyp_molec(2)))
2950 allocate(dhead_scbasej(ntyp_molec(1),ntyp_molec(2)))
2951 allocate(rborn_scbasei(ntyp_molec(1),ntyp_molec(2)))
2952 allocate(rborn_scbasej(ntyp_molec(1),ntyp_molec(2)))
2953 allocate(wdipdip_scbase(3,ntyp_molec(1),ntyp_molec(2)))
2954 allocate(wqdip_scbase(2,ntyp_molec(1),ntyp_molec(2)))
2955 allocate(alphapol_scbase(ntyp_molec(1),ntyp_molec(2)))
2956 allocate(epsintab_scbase(ntyp_molec(1),ntyp_molec(2)))
2959 do i=1,ntyp_molec(1)
2960 do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
2961 write (*,*) "Im in ", i, " ", j
2962 read(isidep_scbase,*) &
2963 eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
2964 chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
2965 write(*,*) "eps",eps_scbase(i,j)
2966 read(isidep_scbase,*) &
2967 (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
2968 chis_scbase(i,j,1),chis_scbase(i,j,2)
2969 read(isidep_scbase,*) &
2970 dhead_scbasei(i,j), &
2971 dhead_scbasej(i,j), &
2972 rborn_scbasei(i,j),rborn_scbasej(i,j)
2973 read(isidep_scbase,*) &
2974 (wdipdip_scbase(k,i,j),k=1,3), &
2975 (wqdip_scbase(k,i,j),k=1,2)
2976 read(isidep_scbase,*) &
2977 alphapol_scbase(i,j), &
2978 epsintab_scbase(i,j)
2981 allocate(aa_scbase(ntyp_molec(1),ntyp_molec(2)))
2982 allocate(bb_scbase(ntyp_molec(1),ntyp_molec(2)))
2984 do i=1,ntyp_molec(1)
2985 do j=1,ntyp_molec(2)-1
2986 epsij=eps_scbase(i,j)
2987 rrij=sigma_scbase(i,j)
2992 sigeps=dsign(1.0D0,epsij)
2994 aa_scbase(i,j)=epsij*rrij*rrij
2995 bb_scbase(i,j)=-sigeps*epsij*rrij
2998 !-----------------READING PEP BASE POTENTIALS-------------------
2999 allocate(eps_pepbase(ntyp_molec(2)))
3000 allocate(sigma_pepbase(ntyp_molec(2)))
3001 allocate(chi_pepbase(ntyp_molec(2),2))
3002 allocate(chipp_pepbase(ntyp_molec(2),2))
3003 allocate(alphasur_pepbase(4,ntyp_molec(2)))
3004 allocate(sigmap1_pepbase(ntyp_molec(2)))
3005 allocate(sigmap2_pepbase(ntyp_molec(2)))
3006 allocate(chis_pepbase(ntyp_molec(2),2))
3007 allocate(wdipdip_pepbase(3,ntyp_molec(2)))
3010 do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
3011 write (*,*) "Im in ", i, " ", j
3012 read(isidep_pepbase,*) &
3013 eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
3014 chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
3015 write(*,*) "eps",eps_pepbase(j)
3016 read(isidep_pepbase,*) &
3017 (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
3018 chis_pepbase(j,1),chis_pepbase(j,2)
3019 read(isidep_pepbase,*) &
3020 (wdipdip_pepbase(k,j),k=1,3)
3022 allocate(aa_pepbase(ntyp_molec(2)))
3023 allocate(bb_pepbase(ntyp_molec(2)))
3025 do j=1,ntyp_molec(2)-1
3026 epsij=eps_pepbase(j)
3027 rrij=sigma_pepbase(j)
3032 sigeps=dsign(1.0D0,epsij)
3034 aa_pepbase(j)=epsij*rrij*rrij
3035 bb_pepbase(j)=-sigeps*epsij*rrij
3037 !--------------READING SC PHOSPHATE-------------------------------------
3038 allocate(eps_scpho(ntyp_molec(1)))
3039 allocate(sigma_scpho(ntyp_molec(1)))
3040 allocate(chi_scpho(ntyp_molec(1),2))
3041 allocate(chipp_scpho(ntyp_molec(1),2))
3042 allocate(alphasur_scpho(4,ntyp_molec(1)))
3043 allocate(sigmap1_scpho(ntyp_molec(1)))
3044 allocate(sigmap2_scpho(ntyp_molec(1)))
3045 allocate(chis_scpho(ntyp_molec(1),2))
3046 allocate(wqq_scpho(ntyp_molec(1)))
3047 allocate(wqdip_scpho(2,ntyp_molec(1)))
3048 allocate(alphapol_scpho(ntyp_molec(1)))
3049 allocate(epsintab_scpho(ntyp_molec(1)))
3050 allocate(dhead_scphoi(ntyp_molec(1)))
3051 allocate(rborn_scphoi(ntyp_molec(1)))
3052 allocate(rborn_scphoj(ntyp_molec(1)))
3053 allocate(alphi_scpho(ntyp_molec(1)))
3057 do j=1,ntyp_molec(1) ! without U then we will take T for U
3058 write (*,*) "Im in scpho ", i, " ", j
3059 read(isidep_scpho,*) &
3060 eps_scpho(j),sigma_scpho(j),chi_scpho(j,1),&
3061 chi_scpho(j,2),chipp_scpho(j,1),chipp_scpho(j,2)
3062 write(*,*) "eps",eps_scpho(j)
3063 read(isidep_scpho,*) &
3064 (alphasur_scpho(k,j),k=1,4),sigmap1_scpho(j),sigmap2_scpho(j), &
3065 chis_scpho(j,1),chis_scpho(j,2)
3066 read(isidep_scpho,*) &
3067 (wqdip_scpho(k,j),k=1,2),wqq_scpho(j),dhead_scphoi(j)
3068 read(isidep_scpho,*) &
3069 epsintab_scpho(j),alphapol_scpho(j),rborn_scphoi(j),rborn_scphoj(j), &
3073 allocate(aa_scpho(ntyp_molec(1)))
3074 allocate(bb_scpho(ntyp_molec(1)))
3076 do j=1,ntyp_molec(1)
3083 sigeps=dsign(1.0D0,epsij)
3085 aa_scpho(j)=epsij*rrij*rrij
3086 bb_scpho(j)=-sigeps*epsij*rrij
3090 read(isidep_peppho,*) &
3091 eps_peppho,sigma_peppho
3092 read(isidep_peppho,*) &
3093 (alphasur_peppho(k),k=1,4),sigmap1_peppho,sigmap2_peppho
3094 read(isidep_peppho,*) &
3095 (wqdip_peppho(k),k=1,2)
3103 sigeps=dsign(1.0D0,epsij)
3105 aa_peppho=epsij*rrij*rrij
3106 bb_peppho=-sigeps*epsij*rrij
3109 allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
3114 ! Define the SC-p interaction constants (hard-coded; old style)
3117 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
3119 ! aad(i,1)=0.3D0*4.0D0**12
3120 ! Following line for constants currently implemented
3121 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
3122 aad(i,1)=1.5D0*4.0D0**12
3123 ! aad(i,1)=0.17D0*5.6D0**12
3125 ! "Soft" SC-p repulsion
3127 ! Following line for constants currently implemented
3128 ! aad(i,1)=0.3D0*4.0D0**6
3129 ! "Hard" SC-p repulsion
3130 bad(i,1)=3.0D0*4.0D0**6
3131 ! bad(i,1)=-2.0D0*0.17D0*5.6D0**6
3140 ! 8/9/01 Read the SC-p interaction constants from file
3143 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
3146 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
3147 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
3148 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
3149 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
3153 write (iout,*) "Parameters of SC-p interactions:"
3155 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
3156 eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
3161 allocate(aad_nucl(ntyp_molec(2)),bad_nucl(ntyp_molec(2))) !(ntyp,2)
3163 do i=1,ntyp_molec(2)
3164 read (iscpp_nucl,*,end=118,err=118) eps_scp_nucl(i),rscp_nucl(i)
3166 do i=1,ntyp_molec(2)
3167 aad_nucl(i)=dabs(eps_scp_nucl(i))*rscp_nucl(i)**12
3168 bad_nucl(i)=-2*eps_scp_nucl(i)*rscp_nucl(i)**6
3170 r0pp=1.12246204830937298142*5.16158
3176 ! Define the constants of the disulfide bridge
3180 ! Old arbitrary potential - commented out.
3185 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
3186 ! energy surface of diethyl disulfide.
3187 ! A. Liwo and U. Kozlowska, 11/24/03
3204 write (iout,'(/a)') "Disulfide bridge parameters:"
3205 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
3206 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
3207 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
3208 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
3211 if (shield_mode.gt.0) then
3212 pi=4.0D0*datan(1.0D0)
3213 !C VSolvSphere the volume of solving sphere
3215 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
3216 !C there will be no distinction between proline peptide group and normal peptide
3217 !C group in case of shielding parameters
3218 VSolvSphere=4.0/3.0*pi*(4.50d0)**3
3219 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3
3220 write (iout,*) VSolvSphere,VSolvSphere_div
3221 !C long axis of side chain
3223 long_r_sidechain(i)=vbldsc0(1,i)
3224 ! if (scelemode.eq.0) then
3225 short_r_sidechain(i)=sigma(i,i)/sqrt(2.0)
3226 if (short_r_sidechain(i).eq.0.0) short_r_sidechain(i)=0.2
3228 ! short_r_sidechain(i)=sigma(i,i)
3230 write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
3237 111 write (iout,*) "Error reading bending energy parameters."
3239 112 write (iout,*) "Error reading rotamer energy parameters."
3241 113 write (iout,*) "Error reading torsional energy parameters."
3243 114 write (iout,*) "Error reading double torsional energy parameters."
3245 115 write (iout,*) &
3246 "Error reading cumulant (multibody energy) parameters."
3248 116 write (iout,*) "Error reading electrostatic energy parameters."
3250 117 write (iout,*) "Error reading side chain interaction parameters."
3252 118 write (iout,*) "Error reading SCp interaction parameters."
3254 119 write (iout,*) "Error reading SCCOR parameters"
3256 121 write (iout,*) "Error in Czybyshev parameters"
3259 call MPI_Finalize(Ierror)
3263 end subroutine parmread
3265 !-----------------------------------------------------------------------------
3267 !-----------------------------------------------------------------------------
3268 subroutine printmat(ldim,m,n,iout,key,a)
3271 character(len=3),dimension(n) :: key
3272 real(kind=8),dimension(ldim,n) :: a
3274 integer :: i,j,k,m,iout,nlim
3278 write (iout,1000) (key(k),k=i,nlim)
3280 1000 format (/5x,8(6x,a3))
3281 1020 format (/80(1h-)/)
3283 write (iout,1010) key(j),(a(j,k),k=i,nlim)
3286 1010 format (a3,2x,8(f9.4))
3288 end subroutine printmat
3289 !-----------------------------------------------------------------------------
3291 !-----------------------------------------------------------------------------
3293 ! Read the PDB file and convert the peptide geometry into virtual-chain
3296 use energy_data, only: itype,istype
3300 ! use control, only: rescode,sugarcode
3301 ! implicit real*8 (a-h,o-z)
3302 ! include 'DIMENSIONS'
3303 ! include 'COMMON.LOCAL'
3304 ! include 'COMMON.VAR'
3305 ! include 'COMMON.CHAIN'
3306 ! include 'COMMON.INTERACT'
3307 ! include 'COMMON.IOUNITS'
3308 ! include 'COMMON.GEO'
3309 ! include 'COMMON.NAMES'
3310 ! include 'COMMON.CONTROL'
3311 ! include 'COMMON.DISTFIT'
3312 ! include 'COMMON.SETUP'
3313 integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift,k!,ity!,&
3315 logical :: lprn=.true.,fail
3316 real(kind=8),dimension(3) :: e1,e2,e3
3317 real(kind=8) :: dcj,efree_temp
3318 character(len=3) :: seq,res,res2
3319 character(len=5) :: atom
3320 character(len=80) :: card
3321 real(kind=8),dimension(3,20) :: sccor
3322 integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin !rescode,
3323 integer :: isugar,molecprev,firstion
3324 character*1 :: sugar
3326 real(kind=8),dimension(3) :: ccc
3328 integer,dimension(2,maxres/3) :: hfrag_alloc
3329 integer,dimension(4,maxres/3) :: bfrag_alloc
3330 real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm)
3331 real(kind=8),dimension(:,:), allocatable :: c_temporary
3332 integer,dimension(:,:) , allocatable :: itype_temporary
3333 integer,dimension(:),allocatable :: istype_temp
3340 ! write (2,*) "UNRES_PDB",unres_pdb
3360 !-----------------------------
3361 allocate(hfrag(2,maxres/3)) !(2,maxres/3)
3362 allocate(bfrag(4,maxres/3)) !(4,maxres/3)
3363 if(.not. allocated(istype)) allocate(istype(maxres))
3365 read (ipdbin,'(a80)',end=10) card
3366 write (iout,'(a)') card
3367 if (card(:5).eq.'HELIX') then
3370 read(card(22:25),*) hfrag(1,nhfrag)
3371 read(card(34:37),*) hfrag(2,nhfrag)
3373 if (card(:5).eq.'SHEET') then
3376 read(card(24:26),*) bfrag(1,nbfrag)
3377 read(card(35:37),*) bfrag(2,nbfrag)
3378 !rc----------------------------------------
3379 !rc to be corrected !!!
3380 bfrag(3,nbfrag)=bfrag(1,nbfrag)
3381 bfrag(4,nbfrag)=bfrag(2,nbfrag)
3382 !rc----------------------------------------
3384 if (card(:3).eq.'END') then
3386 else if (card(:3).eq.'TER') then
3391 itype(ires_old,molecule)=ntyp1_molec(molecule)
3392 itype(ires_old-1,molecule)=ntyp1_molec(molecule)
3393 nres_molec(molecule)=nres_molec(molecule)+2
3395 ! write (iout,*) "Chain ended",ires,ishift,ires_old
3398 dc(j,ires)=sccor(j,iii)
3401 call sccenter(ires,iii,sccor)
3407 if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
3408 ! Fish out the ATOM cards.
3409 ! write(iout,*) 'card',card(1:20)
3410 if (index(card(1:4),'ATOM').gt.0) then
3411 read (card(12:16),*) atom
3412 ! write (iout,*) "! ",atom," !",ires
3413 ! if (atom.eq.'CA' .or. atom.eq.'CH3') then
3414 read (card(23:26),*) ires
3415 read (card(18:20),'(a3)') res
3416 ! write (iout,*) "ires",ires,ires-ishift+ishift1,
3417 ! & " ires_old",ires_old
3418 ! write (iout,*) "ishift",ishift," ishift1",ishift1
3419 ! write (iout,*) "IRES",ires-ishift+ishift1,ires_old
3420 if (ires-ishift+ishift1.ne.ires_old) then
3421 ! Calculate the CM of the preceding residue.
3422 ! if (ibeg.eq.0) call sccenter(ires,iii,sccor)
3424 ! write (iout,*) "Calculating sidechain center iii",iii
3427 dc(j,ires+ishift1-ishift-1)=sccor(j,iii)
3430 call sccenter(ires_old,iii,sccor)
3434 ! Start new residue.
3435 if (res.eq.'Cl-' .or. res.eq.'Na+') then
3438 else if (ibeg.eq.1) then
3439 write (iout,*) "BEG ires",ires
3441 if (res.ne.'GLY' .and. res.ne. 'ACE') then
3444 nres_molec(molecule)=nres_molec(molecule)+1
3446 ires=ires-ishift+ishift1
3448 ! write (iout,*) "ishift",ishift," ires",ires,&
3449 ! " ires_old",ires_old
3451 else if (ibeg.eq.2) then
3453 ishift=-ires_old+ires-1 !!!!!
3454 ishift1=ishift1-1 !!!!!
3455 ! write (iout,*) "New chain started",ires,ishift,ishift1,"!"
3456 ires=ires-ishift+ishift1
3457 ! print *,ires,ishift,ishift1
3461 ishift=ishift-(ires-ishift+ishift1-ires_old-1)
3462 ires=ires-ishift+ishift1
3465 ! print *,'atom',ires,atom
3466 if (res.eq.'ACE' .or. res.eq.'NHE') then
3469 if (atom.eq.'CA '.or.atom.eq.'N ') then
3471 itype(ires,molecule)=rescode(ires,res,0,molecule)
3473 ! nres_molec(molecule)=nres_molec(molecule)+1
3477 itype(ires,molecule)=rescode(ires,res2,0,molecule)
3478 ! nres_molec(molecule)=nres_molec(molecule)+1
3479 read (card(19:19),'(a1)') sugar
3480 isugar=sugarcode(sugar,ires)
3481 ! if (ibeg.eq.1) then
3485 ! print *,"ires=",ires,istype(ires)
3491 ires=ires-ishift+ishift1
3493 ! write (iout,*) "ires_old",ires_old," ires",ires
3494 if (card(27:27).eq."A" .or. card(27:27).eq."B") then
3497 ! write (2,*) "ires",ires," res ",res!," ity"!,ity
3498 if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
3499 res.eq.'NHE'.and.atom(:2).eq.'HN') then
3500 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3501 ! print *,ires,ishift,ishift1
3502 ! write (iout,*) "backbone ",atom
3504 write (iout,'(2i3,2x,a,3f8.3)') &
3505 ires,itype(ires,1),res,(c(j,ires),j=1,3)
3508 nres_molec(molecule)=nres_molec(molecule)+1
3510 sccor(j,iii)=c(j,ires)
3512 else if (.not.unres_pdb .and. (atom.eq."C1'" .or. &
3513 atom.eq."C2'" .or. atom.eq."C3'" &
3514 .or. atom.eq."C4'" .or. atom.eq."O4'")) then
3515 read(card(31:54),'(3f8.3)') (ccc(j),j=1,3)
3516 !c write (2,'(i5,3f10.5)') ires,(ccc(j),j=1,3)
3517 ! print *,ires,ishift,ishift1
3521 ! sccor(j,iii)=c(j,ires)
3524 c(j,ires)=c(j,ires)+ccc(j)/5.0
3526 print *,counter,molecule
3527 if (counter.eq.5) then
3529 nres_molec(molecule)=nres_molec(molecule)+1
3532 ! sccor(j,iii)=c(j,ires)
3536 ! print *, "ATOM",atom(1:3)
3537 else if (atom.eq."C5'") then
3538 read (card(19:19),'(a1)') sugar
3539 isugar=sugarcode(sugar,ires)
3544 ! print *,ires,istype(ires)
3548 ! print *,"nres_molec(molecule)",nres_molec(molecule),ires
3549 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3550 nres_molec(molecule)=nres_molec(molecule)+1
3551 print *,"nres_molec(molecule)",nres_molec(molecule),ires
3555 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
3557 else if ((atom.eq."C1'").and.unres_pdb) then
3559 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
3560 ! write (*,*) card(23:27),ires,itype(ires,1)
3561 else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. &
3562 atom.ne.'N' .and. atom.ne.'C' .and. &
3563 atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. &
3564 atom.ne.'OXT' .and. atom(:2).ne.'3H' &
3565 .and. atom.ne.'P '.and. &
3566 atom(1:1).ne.'H' .and. &
3567 atom.ne.'OP1' .and. atom.ne.'OP2 '.and. atom.ne.'OP3'&
3569 ! write (iout,*) "sidechain ",atom
3570 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
3571 if ((molecule.ne.2).or.(atom(3:3).ne."'")) then
3572 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
3574 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
3577 else if ((ions).and.(card(1:6).eq.'HETATM')) then
3578 if (firstion.eq.0) then
3582 dc(j,ires)=sccor(j,iii)
3585 call sccenter(ires,iii,sccor)
3588 read (card(12:16),*) atom
3589 print *,"HETATOM", atom
3590 read (card(18:20),'(a3)') res
3591 if ((atom(1:2).eq.'NA').or.(atom(1:2).eq.'CL').or.&
3592 (atom(1:2).eq.'CA').or.(atom(1:2).eq.'MG') &
3593 .or.(atom(1:2).eq.'K ')) &
3596 if (molecule.ne.5) molecprev=molecule
3598 nres_molec(molecule)=nres_molec(molecule)+1
3599 print *,"HERE",nres_molec(molecule)
3601 itype(ires,molecule)=rescode(ires,res,0,molecule)
3602 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3606 10 write (iout,'(a,i5)') ' Number of residues found: ',ires
3607 if (ires.eq.0) return
3608 ! Calculate dummy residue coordinates inside the "chain" of a multichain
3611 if (((ires_old.ne.ires).and.(molecule.ne.5)) &
3613 nres_molec(molecule)=nres_molec(molecule)-2
3614 print *,'I have',nres, nres_molec(:)
3616 do k=1,4 ! ions are without dummy
3617 if (nres_molec(k).eq.0) cycle
3619 ! write (iout,*) i,itype(i,1)
3620 ! if (itype(i,1).eq.ntyp1) then
3621 ! write (iout,*) "dummy",i,itype(i,1)
3623 ! c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
3624 ! c(j,i)=(c(j,i-1)+c(j,i+1))/2
3628 if (itype(i,k).eq.ntyp1_molec(k)) then
3629 if (itype(i+1,k).eq.ntyp1_molec(k)) then
3630 if (itype(i+2,k).eq.0) then
3631 ! print *,"masz sieczke"
3633 if (itype(i+2,j).ne.0) then
3635 itype(i+1,j)=ntyp1_molec(j)
3636 nres_molec(k)=nres_molec(k)-1
3637 nres_molec(j)=nres_molec(j)+1
3643 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
3644 ! first is connected prevous chain (itype(i+1,1).eq.ntyp1)=true
3645 ! second dummy atom is conected to next chain itype(i+1,1).eq.ntyp1=false
3646 ! if (unres_pdb) then
3647 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
3648 ! print *,i,'tu dochodze'
3649 ! call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
3657 ! c(j,i)=c(j,i-1)-1.9d0*e2(j)
3661 dcj=(c(j,i-2)-c(j,i-3))/2.0
3662 if (dcj.eq.0) dcj=1.23591524223
3667 else !itype(i+1,1).eq.ntyp1
3668 ! if (unres_pdb) then
3669 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3670 ! call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
3677 c(j,i)=c(j,i+1)-1.9d0*e2(j)
3681 dcj=(c(j,i+3)-c(j,i+2))/2.0
3682 if (dcj.eq.0) dcj=1.23591524223
3687 endif !itype(i+1,1).eq.ntyp1
3688 endif !itype.eq.ntyp1
3692 ! Calculate the CM of the last side chain.
3696 dc(j,ires)=sccor(j,iii)
3699 call sccenter(ires,iii,sccor)
3705 ! print *,"molecule",molecule
3706 if ((itype(nres,1).ne.10)) then
3708 if (molecule.eq.5) molecule=molecprev
3709 itype(nres,molecule)=ntyp1_molec(molecule)
3710 nres_molec(molecule)=nres_molec(molecule)+1
3711 ! if (unres_pdb) then
3712 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
3713 ! call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
3720 ! c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
3724 dcj=(c(j,nres-2)-c(j,nres-3))/2.0
3725 c(j,nres)=c(j,nres-1)+dcj
3726 c(j,2*nres)=c(j,nres)
3730 ! print *,'I have',nres, nres_molec(:)
3732 !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
3734 if (nres.ne.nres0) then
3735 write (iout,*) "Error: wrong parameter value: NRES=",nres,&
3737 stop "Error nres value in WHAM input"
3740 !---------------------------------
3741 !el reallocate tables
3744 ! hfrag_alloc(j,i)=hfrag(j,i)
3747 ! bfrag_alloc(j,i)=bfrag(j,i)
3753 ! allocate(hfrag(2,nres/3)) !(2,maxres/3)
3754 !el allocate(hfrag(2,nhfrag)) !(2,maxres/3)
3755 !el allocate(bfrag(4,nbfrag)) !(4,maxres/3)
3756 ! allocate(bfrag(4,nres/3)) !(4,maxres/3)
3760 ! hfrag(j,i)=hfrag_alloc(j,i)
3765 ! bfrag(j,i)=bfrag_alloc(j,i)
3768 !el end reallocate tables
3769 !---------------------------------
3777 c(j,2*nres)=c(j,nres)
3780 if (itype(1,1).eq.ntyp1) then
3784 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3785 call refsys(2,3,4,e1,e2,e3,fail)
3792 c(j,1)=c(j,2)-1.9d0*e2(j)
3796 dcj=(c(j,4)-c(j,3))/2.0
3802 ! First lets assign correct dummy to correct type of chain
3804 if (itype(1,1).eq.ntyp1) then
3805 if (itype(2,1).eq.0) then
3807 if (itype(2,j).ne.0) then
3809 itype(1,j)=ntyp1_molec(j)
3810 nres_molec(1)=nres_molec(1)-1
3811 nres_molec(j)=nres_molec(j)+1
3818 print *,'I have',nres, nres_molec(:)
3820 ! Copy the coordinates to reference coordinates
3826 ! Calculate internal coordinates.
3828 write (iout,'(/a)') &
3829 "Cartesian coordinates of the reference structure"
3830 write (iout,'(a,16x,3(3x,a5),5x,3(3x,a5))') &
3831 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
3833 write (iout,'(5(a3,1x),i5,3f8.3,5x,3f8.3)') &
3834 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
3835 (c(j,ires+nres),j=1,3)
3838 ! znamy już nres wiec mozna alokowac tablice
3839 ! Calculate internal coordinates.
3840 if(me.eq.king.or..not.out1file)then
3841 write (iout,'(a)') &
3842 "Backbone and SC coordinates as read from the PDB"
3844 write (iout,'(i5,i3,2x,a,3f8.3,5x,3f8.3)') &
3845 ires,itype(ires,1),restyp(itype(ires,1),1),(c(j,ires),j=1,3),&
3846 (c(j,nres+ires),j=1,3)
3849 ! NOW LETS ROCK! SORTING
3850 allocate(c_temporary(3,2*nres))
3851 allocate(itype_temporary(nres,5))
3852 if (.not.allocated(molnum)) allocate(molnum(nres+1))
3853 if (.not.allocated(istype)) write(iout,*) &
3854 "SOMETHING WRONG WITH ISTYTPE"
3855 allocate(istype_temp(nres))
3856 itype_temporary(:,:)=0
3860 if (itype(i,k).ne.0) then
3862 c_temporary(j,seqalingbegin)=c(j,i)
3863 c_temporary(j,seqalingbegin+nres)=c(j,i+nres)
3866 itype_temporary(seqalingbegin,k)=itype(i,k)
3867 print *,i,k,itype(i,k),itype_temporary(seqalingbegin,k),seqalingbegin
3868 istype_temp(seqalingbegin)=istype(i)
3869 molnum(seqalingbegin)=k
3870 seqalingbegin=seqalingbegin+1
3876 c(j,i)=c_temporary(j,i)
3881 itype(i,k)=itype_temporary(i,k)
3882 istype(i)=istype_temp(i)
3885 ! if (itype(1,1).eq.ntyp1) then
3888 ! if (unres_pdb) then
3889 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3890 ! call refsys(2,3,4,e1,e2,e3,fail)
3897 ! c(j,1)=c(j,2)-1.9d0*e2(j)
3901 ! dcj=(c(j,4)-c(j,3))/2.0
3903 ! c(j,nres+1)=c(j,1)
3909 write (iout,'(/a)') &
3910 "Cartesian coordinates of the reference structure after sorting"
3911 write (iout,'(a,16x,3(3x,a5),5x,3(3x,a5))') &
3912 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
3914 write (iout,'(5(a3,1x),i5,3f8.3,5x,3f8.3)') &
3915 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
3916 (c(j,ires+nres),j=1,3)
3920 ! print *,seqalingbegin,nres
3921 if(.not.allocated(vbld)) then
3922 allocate(vbld(2*nres))
3927 if(.not.allocated(vbld_inv)) then
3928 allocate(vbld_inv(2*nres))
3934 if(.not.allocated(theta)) then
3935 allocate(theta(nres+2))
3939 if(.not.allocated(phi)) allocate(phi(nres+2))
3940 if(.not.allocated(alph)) allocate(alph(nres+2))
3941 if(.not.allocated(omeg)) allocate(omeg(nres+2))
3942 if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
3943 if(.not.allocated(phiref)) allocate(phiref(nres+2))
3944 if(.not.allocated(costtab)) allocate(costtab(nres))
3945 if(.not.allocated(sinttab)) allocate(sinttab(nres))
3946 if(.not.allocated(cost2tab)) allocate(cost2tab(nres))
3947 if(.not.allocated(sint2tab)) allocate(sint2tab(nres))
3948 if(.not.allocated(xxref)) allocate(xxref(nres))
3949 if(.not.allocated(yyref)) allocate(yyref(nres))
3950 if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
3951 if(.not.allocated(dc_norm)) then
3952 ! if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres+2))
3953 allocate(dc_norm(3,0:2*nres+2))
3957 call int_from_cart(.true.,.false.)
3958 call sc_loc_geom(.false.)
3960 thetaref(i)=theta(i)
3970 dc(j,i)=c(j,i+1)-c(j,i)
3971 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
3976 dc(j,i+nres)=c(j,i+nres)-c(j,i)
3977 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
3979 ! write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),&
3983 ! Copy the coordinates to reference coordinates
3984 ! Splits to single chain if occurs
3988 ! cref(j,i,cou)=c(j,i)
3992 if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
3993 if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
3994 if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
3995 !-----------------------------
3999 write (iout,*) "symetr", symetr
4002 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
4004 if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then
4007 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
4012 cref(j,i,cou)=c(j,i)
4013 cref(j,i+nres,cou)=c(j,i+nres)
4015 chain_rep(j,lll,kkk)=c(j,i)
4016 chain_rep(j,lll+nres,kkk)=c(j,i+nres)
4020 write (iout,*) chain_length
4021 if (chain_length.eq.0) chain_length=nres
4023 chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
4024 chain_rep(j,chain_length+nres,symetr) &
4025 =chain_rep(j,chain_length+nres,1)
4028 ! write (iout,*) "spraw lancuchy",chain_length,symetr
4030 ! do kkk=1,chain_length
4031 ! write (iout,*) itype(kkk,1),(chain_rep(j,kkk,i), j=1,3)
4035 ! makes copy of chains
4036 write (iout,*) "symetr", symetr
4041 if (symetr.gt.1) then
4048 write(iout,*) (tabperm(i,kkk),kkk=1,4)
4054 write (iout,*) i,icha
4055 do lll=1,chain_length
4057 if (cou.le.nres) then
4059 kupa=mod(lll,chain_length)
4060 iprzes=(kkk-1)*chain_length+lll
4061 if (kupa.eq.0) kupa=chain_length
4062 write (iout,*) "kupa", kupa
4063 cref(j,iprzes,i)=chain_rep(j,kupa,icha)
4064 cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
4071 !-koniec robienia kopii
4074 write (iout,*) "nowa struktura", nperm
4076 write (iout,110) restyp(itype(i,1),1),i,cref(1,i,kkk),&
4078 cref(3,i,kkk),cref(1,nres+i,kkk),&
4079 cref(2,nres+i,kkk),cref(3,nres+i,kkk)
4081 100 format (//' alpha-carbon coordinates ',&
4082 ' centroid coordinates'/ &
4083 ' ', 6X,'X',11X,'Y',11X,'Z', &
4084 10X,'X',11X,'Y',11X,'Z')
4085 110 format (a,'(',i5,')',6f12.5)
4091 bfrag(i,j)=bfrag(i,j)-ishift
4097 hfrag(i,j)=hfrag(i,j)-ishift
4103 end subroutine readpdb
4104 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
4105 !-----------------------------------------------------------------------------
4107 !-----------------------------------------------------------------------------
4108 subroutine read_control
4122 use random, only: random_init
4123 ! implicit real*8 (a-h,o-z)
4124 ! include 'DIMENSIONS'
4126 use prng, only:prng_restart
4128 logical :: OKRandom!, prng_restart
4131 ! include 'COMMON.IOUNITS'
4132 ! include 'COMMON.TIME1'
4133 ! include 'COMMON.THREAD'
4134 ! include 'COMMON.SBRIDGE'
4135 ! include 'COMMON.CONTROL'
4136 ! include 'COMMON.MCM'
4137 ! include 'COMMON.MAP'
4138 ! include 'COMMON.HEADER'
4139 ! include 'COMMON.CSA'
4140 ! include 'COMMON.CHAIN'
4141 ! include 'COMMON.MUCA'
4142 ! include 'COMMON.MD'
4143 ! include 'COMMON.FFIELD'
4144 ! include 'COMMON.INTERACT'
4145 ! include 'COMMON.SETUP'
4146 !el integer :: KDIAG,ICORFL,IXDR
4147 !el COMMON /MACHSW/ KDIAG,ICORFL,IXDR
4148 character(len=8),dimension(0:3) :: diagmeth = reshape((/'Library ',&
4149 'EVVRSP ','Givens ','Jacobi '/),shape(diagmeth))
4150 ! character(len=80) :: ucase
4151 character(len=640) :: controlcard
4153 real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
4159 read (INP,'(a)') titel
4160 call card_concat(controlcard,.true.)
4161 ! out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
4162 ! print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
4163 call reada(controlcard,'SEED',seed,0.0D0)
4164 call random_init(seed)
4165 ! Set up the time limit (caution! The time must be input in minutes!)
4166 read_cart=index(controlcard,'READ_CART').gt.0
4167 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
4168 call readi(controlcard,'SYM',symetr,1)
4169 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
4170 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
4171 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
4172 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
4173 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
4174 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
4175 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
4176 call reada(controlcard,'DRMS',drms,0.1D0)
4177 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4178 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
4179 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
4180 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
4181 write (iout,'(a,f10.1)')'DRMS = ',drms
4182 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
4183 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
4185 call readi(controlcard,'NZ_START',nz_start,0)
4186 call readi(controlcard,'NZ_END',nz_end,0)
4187 call readi(controlcard,'IZ_SC',iz_sc,0)
4188 timlim=60.0D0*timlim
4189 safety = 60.0d0*safety
4192 call reada(controlcard,"T_BATH",t_bath,300.0d0)
4193 !C SHIELD keyword sets if the shielding effect of side-chains is used
4194 !C 0 denots no shielding is used all peptide are equally despite the
4195 !C solvent accesible area
4196 !C 1 the newly introduced function
4197 !C 2 reseved for further possible developement
4198 call readi(controlcard,'SHIELD',shield_mode,0)
4199 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4200 write(iout,*) "shield_mode",shield_mode
4201 call readi(controlcard,'TORMODE',tor_mode,0)
4202 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4203 write(iout,*) "torsional and valence angle mode",tor_mode
4205 !C Varibles set size of box
4206 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
4207 protein=index(controlcard,"PROTEIN").gt.0
4208 ions=index(controlcard,"IONS").gt.0
4209 nucleic=index(controlcard,"NUCLEIC").gt.0
4210 write (iout,*) "with_theta_constr ",with_theta_constr
4211 AFMlog=(index(controlcard,'AFM'))
4212 selfguide=(index(controlcard,'SELFGUIDE'))
4213 print *,'AFMlog',AFMlog,selfguide,"KUPA"
4214 call readi(controlcard,'GENCONSTR',genconstr,0)
4215 call reada(controlcard,'BOXX',boxxsize,100.0d0)
4216 call reada(controlcard,'BOXY',boxysize,100.0d0)
4217 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
4218 call readi(controlcard,'TUBEMOD',tubemode,0)
4219 print *,"SCELE",scelemode
4220 call readi(controlcard,"SCELEMODE",scelemode,0)
4221 print *,"SCELE",scelemode
4223 ! elemode = 0 is orignal UNRES electrostatics
4224 ! elemode = 1 is "Momo" potentials in progress
4225 ! elemode = 2 is in development EVALD
4226 write (iout,*) TUBEmode,"TUBEMODE"
4227 if (TUBEmode.gt.0) then
4228 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
4229 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
4230 call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
4231 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
4232 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
4233 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
4234 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
4235 buftubebot=bordtubebot+tubebufthick
4236 buftubetop=bordtubetop-tubebufthick
4239 ! CUTOFFF ON ELECTROSTATICS
4240 call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
4241 call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
4242 write(iout,*) "R_CUT_ELE=",r_cut_ele
4243 ! Lipidic parameters
4244 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
4245 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
4246 if (lipthick.gt.0.0d0) then
4247 bordliptop=(boxzsize+lipthick)/2.0
4248 bordlipbot=bordliptop-lipthick
4249 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
4250 write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
4251 buflipbot=bordlipbot+lipbufthick
4252 bufliptop=bordliptop-lipbufthick
4253 if ((lipbufthick*2.0d0).gt.lipthick) &
4254 write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
4255 endif !lipthick.gt.0
4256 write(iout,*) "bordliptop=",bordliptop
4257 write(iout,*) "bordlipbot=",bordlipbot
4258 write(iout,*) "bufliptop=",bufliptop
4259 write(iout,*) "buflipbot=",buflipbot
4260 write (iout,*) "SHIELD MODE",shield_mode
4262 !C-------------------------
4263 minim=(index(controlcard,'MINIMIZE').gt.0)
4264 dccart=(index(controlcard,'CART').gt.0)
4265 overlapsc=(index(controlcard,'OVERLAP').gt.0)
4266 overlapsc=.not.overlapsc
4267 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
4268 searchsc=.not.searchsc
4269 sideadd=(index(controlcard,'SIDEADD').gt.0)
4270 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
4271 outpdb=(index(controlcard,'PDBOUT').gt.0)
4272 outmol2=(index(controlcard,'MOL2OUT').gt.0)
4273 pdbref=(index(controlcard,'PDBREF').gt.0)
4274 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
4275 indpdb=index(controlcard,'PDBSTART')
4276 extconf=(index(controlcard,'EXTCONF').gt.0)
4277 call readi(controlcard,'IPRINT',iprint,0)
4278 call readi(controlcard,'MAXGEN',maxgen,10000)
4279 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
4280 call readi(controlcard,"KDIAG",kdiag,0)
4281 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
4282 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) &
4283 write (iout,*) "RESCALE_MODE",rescale_mode
4284 split_ene=index(controlcard,'SPLIT_ENE').gt.0
4285 if (index(controlcard,'REGULAR').gt.0.0D0) then
4286 call reada(controlcard,'WEIDIS',weidis,0.1D0)
4290 if (index(controlcard,'CHECKGRAD').gt.0) then
4292 if (index(controlcard,'CART').gt.0) then
4294 elseif (index(controlcard,'CARINT').gt.0) then
4299 elseif (index(controlcard,'THREAD').gt.0) then
4301 call readi(controlcard,'THREAD',nthread,0)
4302 if (nthread.gt.0) then
4303 call reada(controlcard,'WEIDIS',weidis,0.1D0)
4306 write (iout,'(a)')'A number has to follow the THREAD keyword.'
4307 stop 'Error termination in Read_Control.'
4309 else if (index(controlcard,'MCMA').gt.0) then
4311 else if (index(controlcard,'MCEE').gt.0) then
4313 else if (index(controlcard,'MULTCONF').gt.0) then
4315 else if (index(controlcard,'MAP').gt.0) then
4317 call readi(controlcard,'MAP',nmap,0)
4318 else if (index(controlcard,'CSA').gt.0) then
4320 !rc else if (index(controlcard,'ZSCORE').gt.0) then
4322 !rc ZSCORE is rm from UNRES, modecalc=9 is available
4325 !fcm else if (index(controlcard,'MCMF').gt.0) then
4327 else if (index(controlcard,'SOFTREG').gt.0) then
4329 else if (index(controlcard,'CHECK_BOND').gt.0) then
4331 else if (index(controlcard,'TEST').gt.0) then
4333 else if (index(controlcard,'MD').gt.0) then
4335 else if (index(controlcard,'RE ').gt.0) then
4339 lmuca=index(controlcard,'MUCA').gt.0
4340 call readi(controlcard,'MUCADYN',mucadyn,0)
4341 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
4342 if (lmuca .and. (me.eq.king .or. .not.out1file )) &
4344 write (iout,*) 'MUCADYN=',mucadyn
4345 write (iout,*) 'MUCASMOOTH=',muca_smooth
4348 iscode=index(controlcard,'ONE_LETTER')
4349 indphi=index(controlcard,'PHI')
4350 indback=index(controlcard,'BACK')
4351 iranconf=index(controlcard,'RAND_CONF')
4352 i2ndstr=index(controlcard,'USE_SEC_PRED')
4353 gradout=index(controlcard,'GRADOUT').gt.0
4354 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
4355 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
4356 if (me.eq.king .or. .not.out1file ) &
4357 write (iout,*) "DISTCHAINMAX",distchainmax
4359 if(me.eq.king.or..not.out1file) &
4360 write (iout,'(2a)') diagmeth(kdiag),&
4361 ' routine used to diagonalize matrices.'
4362 if (shield_mode.gt.0) then
4363 pi=4.0D0*datan(1.0D0)
4364 !C VSolvSphere the volume of solving sphere
4366 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
4367 !C there will be no distinction between proline peptide group and normal peptide
4368 !C group in case of shielding parameters
4369 VSolvSphere=4.0/3.0*pi*(4.50d0)**3
4370 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3
4371 write (iout,*) VSolvSphere,VSolvSphere_div
4372 !C long axis of side chain
4374 ! long_r_sidechain(i)=vbldsc0(1,i)
4375 ! short_r_sidechain(i)=sigma0(i)
4376 ! write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
4382 end subroutine read_control
4383 !-----------------------------------------------------------------------------
4384 subroutine read_REMDpar
4386 ! Read REMD settings
4393 use control_data, only:out1file
4394 ! implicit real*8 (a-h,o-z)
4395 ! include 'DIMENSIONS'
4396 ! include 'COMMON.IOUNITS'
4397 ! include 'COMMON.TIME1'
4398 ! include 'COMMON.MD'
4401 !el include 'COMMON.LANGEVIN'
4403 !el include 'COMMON.LANGEVIN.lang0'
4405 ! include 'COMMON.INTERACT'
4406 ! include 'COMMON.NAMES'
4407 ! include 'COMMON.GEO'
4408 ! include 'COMMON.REMD'
4409 ! include 'COMMON.CONTROL'
4410 ! include 'COMMON.SETUP'
4411 ! character(len=80) :: ucase
4412 character(len=320) :: controlcard
4413 character(len=3200) :: controlcard1
4414 integer :: iremd_m_total
4417 ! real(kind=8) :: var,ene
4419 if(me.eq.king.or..not.out1file) &
4420 write (iout,*) "REMD setup"
4422 call card_concat(controlcard,.true.)
4423 call readi(controlcard,"NREP",nrep,3)
4424 call readi(controlcard,"NSTEX",nstex,1000)
4425 call reada(controlcard,"RETMIN",retmin,10.0d0)
4426 call reada(controlcard,"RETMAX",retmax,1000.0d0)
4427 mremdsync=(index(controlcard,'SYNC').gt.0)
4428 call readi(controlcard,"NSYN",i_sync_step,100)
4429 restart1file=(index(controlcard,'REST1FILE').gt.0)
4430 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
4431 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
4432 if(max_cache_traj_use.gt.max_cache_traj) &
4433 max_cache_traj_use=max_cache_traj
4434 if(me.eq.king.or..not.out1file) then
4435 !d if (traj1file) then
4436 !rc caching is in testing - NTWX is not ignored
4437 !d write (iout,*) "NTWX value is ignored"
4438 !d write (iout,*) " trajectory is stored to one file by master"
4439 !d write (iout,*) " before exchange at NSTEX intervals"
4441 write (iout,*) "NREP= ",nrep
4442 write (iout,*) "NSTEX= ",nstex
4443 write (iout,*) "SYNC= ",mremdsync
4444 write (iout,*) "NSYN= ",i_sync_step
4445 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
4448 allocate(remd_t(nrep),remd_m(nrep)) !(maxprocs)
4449 if (index(controlcard,'TLIST').gt.0) then
4451 call card_concat(controlcard1,.true.)
4452 read(controlcard1,*) (remd_t(i),i=1,nrep)
4453 if(me.eq.king.or..not.out1file) &
4454 write (iout,*)'tlist',(remd_t(i),i=1,nrep)
4457 if (index(controlcard,'MLIST').gt.0) then
4459 call card_concat(controlcard1,.true.)
4460 read(controlcard1,*) (remd_m(i),i=1,nrep)
4461 if(me.eq.king.or..not.out1file) then
4462 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
4465 iremd_m_total=iremd_m_total+remd_m(i)
4467 write (iout,*) 'Total number of replicas ',iremd_m_total
4470 if(me.eq.king.or..not.out1file) &
4471 write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
4473 end subroutine read_REMDpar
4474 !-----------------------------------------------------------------------------
4475 subroutine read_MDpar
4479 use control_data, only: r_cut,rlamb,out1file
4481 use geometry_data, only: pi
4483 ! implicit real*8 (a-h,o-z)
4484 ! include 'DIMENSIONS'
4485 ! include 'COMMON.IOUNITS'
4486 ! include 'COMMON.TIME1'
4487 ! include 'COMMON.MD'
4490 !el include 'COMMON.LANGEVIN'
4492 !el include 'COMMON.LANGEVIN.lang0'
4494 ! include 'COMMON.INTERACT'
4495 ! include 'COMMON.NAMES'
4496 ! include 'COMMON.GEO'
4497 ! include 'COMMON.SETUP'
4498 ! include 'COMMON.CONTROL'
4499 ! include 'COMMON.SPLITELE'
4500 ! character(len=80) :: ucase
4501 character(len=320) :: controlcard
4506 call card_concat(controlcard,.true.)
4507 call readi(controlcard,"NSTEP",n_timestep,1000000)
4508 call readi(controlcard,"NTWE",ntwe,100)
4509 call readi(controlcard,"NTWX",ntwx,1000)
4510 call reada(controlcard,"DT",d_time,1.0d-1)
4511 call reada(controlcard,"DVMAX",dvmax,2.0d1)
4512 call reada(controlcard,"DAMAX",damax,1.0d1)
4513 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
4514 call readi(controlcard,"LANG",lang,0)
4515 RESPA = index(controlcard,"RESPA") .gt. 0
4516 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
4517 ntime_split0=ntime_split
4518 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
4519 ntime_split0=ntime_split
4520 call reada(controlcard,"R_CUT",r_cut,2.0d0)
4521 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
4522 rest = index(controlcard,"REST").gt.0
4523 tbf = index(controlcard,"TBF").gt.0
4524 usampl = index(controlcard,"USAMPL").gt.0
4525 mdpdb = index(controlcard,"MDPDB").gt.0
4526 call reada(controlcard,"T_BATH",t_bath,300.0d0)
4527 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
4528 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
4529 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
4530 if (count_reset_moment.eq.0) count_reset_moment=1000000000
4531 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
4532 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
4533 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
4534 if (count_reset_vel.eq.0) count_reset_vel=1000000000
4535 large = index(controlcard,"LARGE").gt.0
4536 print_compon = index(controlcard,"PRINT_COMPON").gt.0
4537 rattle = index(controlcard,"RATTLE").gt.0
4538 preminim=(index(controlcard,'PREMINIM').gt.0)
4539 write (iout,*) "PREMINIM ",preminim
4540 dccart=(index(controlcard,'CART').gt.0)
4541 if (preminim) call read_minim
4542 ! if performing umbrella sampling, fragments constrained are read from the fragment file
4548 if(me.eq.king.or..not.out1file) then
4550 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
4552 write (iout,'(a)') "The units are:"
4553 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
4554 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",&
4555 " acceleration: angstrom/(48.9 fs)**2"
4556 write (iout,'(a)') "energy: kcal/mol, temperature: K"
4558 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
4559 write (iout,'(a60,f10.5,a)') &
4560 "Initial time step of numerical integration:",d_time,&
4562 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
4564 write (iout,'(2a,i4,a)') &
4565 "A-MTS algorithm used; initial time step for fast-varying",&
4566 " short-range forces split into",ntime_split," steps."
4567 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",&
4568 r_cut," lambda",rlamb
4570 write (iout,'(2a,f10.5)') &
4571 "Maximum acceleration threshold to reduce the time step",&
4572 "/increase split number:",damax
4573 write (iout,'(2a,f10.5)') &
4574 "Maximum predicted energy drift to reduce the timestep",&
4575 "/increase split number:",edriftmax
4576 write (iout,'(a60,f10.5)') &
4577 "Maximum velocity threshold to reduce velocities:",dvmax
4578 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
4579 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
4580 if (rattle) write (iout,'(a60)') &
4581 "Rattle algorithm used to constrain the virtual bonds"
4585 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
4586 call reada(controlcard,"RWAT",rwat,1.4d0)
4587 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
4588 surfarea=index(controlcard,"SURFAREA").gt.0
4589 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
4590 if(me.eq.king.or..not.out1file)then
4591 write (iout,'(/a,$)') "Langevin dynamics calculation"
4593 write (iout,'(a/)') &
4594 " with direct integration of Langevin equations"
4595 else if (lang.eq.2) then
4596 write (iout,'(a/)') " with TINKER stochasic MD integrator"
4597 else if (lang.eq.3) then
4598 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
4599 else if (lang.eq.4) then
4600 write (iout,'(a/)') " in overdamped mode"
4602 write (iout,'(//a,i5)') &
4603 "=========== ERROR: Unknown Langevin dynamics mode:",lang
4606 write (iout,'(a60,f10.5)') "Temperature:",t_bath
4607 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
4608 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
4609 write (iout,'(a60,f10.5)') &
4610 "Scaling factor of the friction forces:",scal_fric
4611 if (surfarea) write (iout,'(2a,i10,a)') &
4612 "Friction coefficients will be scaled by solvent-accessible",&
4613 " surface area every",reset_fricmat," steps."
4615 ! Calculate friction coefficients and bounds of stochastic forces
4616 eta=6*pi*cPoise*etawat
4617 if(me.eq.king.or..not.out1file) &
4618 write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
4621 do j=1,5 !types of molecules
4622 gamp(j)=scal_fric*(pstok(j)+rwat)*eta
4623 stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
4625 allocate(gamsc(ntyp1,5),stdfsc(ntyp1,5)) !(ntyp1)
4626 do j=1,5 !types of molecules
4628 gamsc(i,j)=scal_fric*(restok(i,j)+rwat)*eta
4629 stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
4633 if(me.eq.king.or..not.out1file)then
4634 write (iout,'(/2a/)') &
4635 "Radii of site types and friction coefficients and std's of",&
4636 " stochastic forces of fully exposed sites"
4637 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp(1),stdfp*dsqrt(gamp(1))
4639 write (iout,'(a5,f5.2,2f10.5)') restyp(i,1),restok(i,1),&
4640 gamsc(i,1),stdfsc(i,1)*dsqrt(gamsc(i,1))
4644 if(me.eq.king.or..not.out1file)then
4645 write (iout,'(a)') "Berendsen bath calculation"
4646 write (iout,'(a60,f10.5)') "Temperature:",t_bath
4647 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
4649 write (iout,'(a,i10,a)') "Momenta will be reset at zero every",&
4650 count_reset_moment," steps"
4652 write (iout,'(a,i10,a)') &
4653 "Velocities will be reset at random every",count_reset_vel,&
4657 if(me.eq.king.or..not.out1file) &
4658 write (iout,'(a31)') "Microcanonical mode calculation"
4660 if(me.eq.king.or..not.out1file)then
4661 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
4663 write(iout,*) "MD running with constraints."
4664 write(iout,*) "Equilibration time ", eq_time, " mtus."
4665 write(iout,*) "Constraining ", nfrag," fragments."
4666 write(iout,*) "Length of each fragment, weight and q0:"
4668 write (iout,*) "Set of restraints #",iset
4670 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),&
4671 ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
4673 write(iout,*) "constraints between ", npair, "fragments."
4674 write(iout,*) "constraint pairs, weights and q0:"
4676 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),&
4677 ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
4679 write(iout,*) "angle constraints within ", nfrag_back,&
4680 "backbone fragments."
4681 write(iout,*) "fragment, weights:"
4683 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),&
4684 ifrag_back(2,i,iset),wfrag_back(1,i,iset),&
4685 wfrag_back(2,i,iset),wfrag_back(3,i,iset)
4688 iset=mod(kolor,nset)+1
4691 if(me.eq.king.or..not.out1file) &
4692 write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
4694 end subroutine read_MDpar
4695 !-----------------------------------------------------------------------------
4699 ! implicit real*8 (a-h,o-z)
4700 ! include 'DIMENSIONS'
4701 ! include 'COMMON.MAP'
4702 ! include 'COMMON.IOUNITS'
4703 character(len=3) :: angid(4) = (/'THE','PHI','ALP','OME'/)
4704 character(len=80) :: mapcard !,ucase
4707 ! real(kind=8) :: var,ene
4710 read (inp,'(a)') mapcard
4711 mapcard=ucase(mapcard)
4712 if (index(mapcard,'PHI').gt.0) then
4714 else if (index(mapcard,'THE').gt.0) then
4716 else if (index(mapcard,'ALP').gt.0) then
4718 else if (index(mapcard,'OME').gt.0) then
4721 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
4722 stop 'Error - illegal variable spec in MAP card.'
4724 call readi (mapcard,'RES1',res1(imap),0)
4725 call readi (mapcard,'RES2',res2(imap),0)
4726 if (res1(imap).eq.0) then
4727 res1(imap)=res2(imap)
4728 else if (res2(imap).eq.0) then
4729 res2(imap)=res1(imap)
4731 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
4732 write (iout,'(a)') &
4733 'Error - illegal definition of variable group in MAP.'
4734 stop 'Error - illegal definition of variable group in MAP.'
4736 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
4737 call reada(mapcard,'TO',ang_to(imap),0.0D0)
4738 call readi(mapcard,'NSTEP',nstep(imap),0)
4739 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
4740 write (iout,'(a)') &
4741 'Illegal boundary and/or step size specification in MAP.'
4742 stop 'Illegal boundary and/or step size specification in MAP.'
4746 end subroutine map_read
4747 !-----------------------------------------------------------------------------
4750 use control_data, only: vdisulf
4752 ! implicit real*8 (a-h,o-z)
4753 ! include 'DIMENSIONS'
4754 ! include 'COMMON.IOUNITS'
4755 ! include 'COMMON.GEO'
4756 ! include 'COMMON.CSA'
4757 ! include 'COMMON.BANK'
4758 ! include 'COMMON.CONTROL'
4759 ! character(len=80) :: ucase
4760 character(len=620) :: mcmcard
4762 ! integer :: ntf,ik,iw_pdb
4763 ! real(kind=8) :: var,ene
4765 call card_concat(mcmcard,.true.)
4767 call readi(mcmcard,'NCONF',nconf,50)
4768 call readi(mcmcard,'NADD',nadd,0)
4769 call readi(mcmcard,'JSTART',jstart,1)
4770 call readi(mcmcard,'JEND',jend,1)
4771 call readi(mcmcard,'NSTMAX',nstmax,500000)
4772 call readi(mcmcard,'N0',n0,1)
4773 call readi(mcmcard,'N1',n1,6)
4774 call readi(mcmcard,'N2',n2,4)
4775 call readi(mcmcard,'N3',n3,0)
4776 call readi(mcmcard,'N4',n4,0)
4777 call readi(mcmcard,'N5',n5,0)
4778 call readi(mcmcard,'N6',n6,10)
4779 call readi(mcmcard,'N7',n7,0)
4780 call readi(mcmcard,'N8',n8,0)
4781 call readi(mcmcard,'N9',n9,0)
4782 call readi(mcmcard,'N14',n14,0)
4783 call readi(mcmcard,'N15',n15,0)
4784 call readi(mcmcard,'N16',n16,0)
4785 call readi(mcmcard,'N17',n17,0)
4786 call readi(mcmcard,'N18',n18,0)
4788 vdisulf=(index(mcmcard,'DYNSS').gt.0)
4790 call readi(mcmcard,'NDIFF',ndiff,2)
4791 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
4792 call readi(mcmcard,'IS1',is1,1)
4793 call readi(mcmcard,'IS2',is2,8)
4794 call readi(mcmcard,'NRAN0',nran0,4)
4795 call readi(mcmcard,'NRAN1',nran1,2)
4796 call readi(mcmcard,'IRR',irr,1)
4797 call readi(mcmcard,'NSEED',nseed,20)
4798 call readi(mcmcard,'NTOTAL',ntotal,10000)
4799 call reada(mcmcard,'CUT1',cut1,2.0d0)
4800 call reada(mcmcard,'CUT2',cut2,5.0d0)
4801 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
4802 call readi(mcmcard,'ICMAX',icmax,3)
4803 call readi(mcmcard,'IRESTART',irestart,0)
4804 !!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
4807 call reada(mcmcard,'DELE',dele,20.0d0)
4808 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
4809 call readi(mcmcard,'IREF',iref,0)
4810 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
4811 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
4812 call readi(mcmcard,'NCONF_IN',nconf_in,0)
4813 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
4814 write (iout,*) "NCONF_IN",nconf_in
4816 end subroutine csaread
4817 !-----------------------------------------------------------------------------
4821 use control_data, only: MaxMoveType
4824 ! implicit real*8 (a-h,o-z)
4825 ! include 'DIMENSIONS'
4826 ! include 'COMMON.MCM'
4827 ! include 'COMMON.MCE'
4828 ! include 'COMMON.IOUNITS'
4829 ! character(len=80) :: ucase
4830 character(len=320) :: mcmcard
4833 ! real(kind=8) :: var,ene
4835 call card_concat(mcmcard,.true.)
4836 call readi(mcmcard,'MAXACC',maxacc,100)
4837 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
4838 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
4839 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
4840 call readi(mcmcard,'MAXREPM',maxrepm,200)
4841 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
4842 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
4843 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
4844 call reada(mcmcard,'E_UP',e_up,5.0D0)
4845 call reada(mcmcard,'DELTE',delte,0.1D0)
4846 call readi(mcmcard,'NSWEEP',nsweep,5)
4847 call readi(mcmcard,'NSTEPH',nsteph,0)
4848 call readi(mcmcard,'NSTEPC',nstepc,0)
4849 call reada(mcmcard,'TMIN',tmin,298.0D0)
4850 call reada(mcmcard,'TMAX',tmax,298.0D0)
4851 call readi(mcmcard,'NWINDOW',nwindow,0)
4852 call readi(mcmcard,'PRINT_MC',print_mc,0)
4853 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
4854 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
4855 ent_read=(index(mcmcard,'ENT_READ').gt.0)
4856 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
4857 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
4858 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
4859 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
4860 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
4861 if (nwindow.gt.0) then
4862 allocate(winstart(nwindow)) !!el (maxres)
4863 allocate(winend(nwindow)) !!el
4864 allocate(winlen(nwindow)) !!el
4865 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
4867 winlen(i)=winend(i)-winstart(i)+1
4870 if (tmax.lt.tmin) tmax=tmin
4871 if (tmax.eq.tmin) then
4875 if (nstepc.gt.0 .and. nsteph.gt.0) then
4876 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
4877 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
4879 allocate(sumpro_type(0:MaxMoveType)) !(0:MaxMoveType)
4880 ! Probabilities of different move types
4881 sumpro_type(0)=0.0D0
4882 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
4883 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
4884 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
4885 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
4886 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
4887 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
4888 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
4890 print *,'i',i,' sumprotype',sumpro_type(i)
4891 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
4892 print *,'i',i,' sumprotype',sumpro_type(i)
4895 end subroutine mcmread
4896 !-----------------------------------------------------------------------------
4897 subroutine read_minim
4900 ! implicit real*8 (a-h,o-z)
4901 ! include 'DIMENSIONS'
4902 ! include 'COMMON.MINIM'
4903 ! include 'COMMON.IOUNITS'
4904 ! character(len=80) :: ucase
4905 character(len=320) :: minimcard
4907 ! integer :: ntf,ik,iw_pdb
4908 ! real(kind=8) :: var,ene
4910 call card_concat(minimcard,.true.)
4911 call readi(minimcard,'MAXMIN',maxmin,2000)
4912 call readi(minimcard,'MAXFUN',maxfun,5000)
4913 call readi(minimcard,'MINMIN',minmin,maxmin)
4914 call readi(minimcard,'MINFUN',minfun,maxmin)
4915 call reada(minimcard,'TOLF',tolf,1.0D-2)
4916 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
4917 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
4918 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
4919 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
4920 write (iout,'(/80(1h*)/20x,a/80(1h*))') &
4921 'Options in energy minimization:'
4922 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') &
4923 'MaxMin:',MaxMin,' MaxFun:',MaxFun,&
4924 'MinMin:',MinMin,' MinFun:',MinFun,&
4925 ' TolF:',TolF,' RTolF:',RTolF
4927 end subroutine read_minim
4928 !-----------------------------------------------------------------------------
4929 subroutine openunits
4931 use MD_data, only: usampl
4934 use control_data, only:out1file
4935 use control, only: getenv_loc
4936 ! implicit real*8 (a-h,o-z)
4937 ! include 'DIMENSIONS'
4940 character(len=16) :: form,nodename
4941 integer :: nodelen,ierror,npos
4943 ! include 'COMMON.SETUP'
4944 ! include 'COMMON.IOUNITS'
4945 ! include 'COMMON.MD'
4946 ! include 'COMMON.CONTROL'
4947 integer :: lenpre,lenpot,lentmp !,ilen
4949 character(len=3) :: out1file_text !,ucase
4950 character(len=3) :: ll
4953 ! integer :: ntf,ik,iw_pdb
4954 ! real(kind=8) :: var,ene
4956 ! print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
4957 call getenv_loc("PREFIX",prefix)
4959 call getenv_loc("POT",pot)
4960 call getenv_loc("DIRTMP",tmpdir)
4961 call getenv_loc("CURDIR",curdir)
4962 call getenv_loc("OUT1FILE",out1file_text)
4963 ! print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
4964 out1file_text=ucase(out1file_text)
4965 if (out1file_text(1:1).eq."Y") then
4968 out1file=fg_rank.gt.0
4973 if (lentmp.gt.0) then
4974 write (*,'(80(1h!))')
4975 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
4976 write (*,'(80(1h!))')
4977 write (*,*)"All output files will be on node /tmp directory."
4979 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
4980 if (me.eq.king) then
4981 write (*,*) "The master node is ",nodename
4982 else if (fg_rank.eq.0) then
4983 write (*,*) "I am the CG slave node ",nodename
4985 write (*,*) "I am the FG slave node ",nodename
4988 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
4989 lenpre = lentmp+lenpre+1
4991 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
4992 ! Get the names and open the input files
4993 #if defined(WINIFL) || defined(WINPGI)
4994 open(1,file=pref_orig(:ilen(pref_orig))// &
4995 '.inp',status='old',readonly,shared)
4996 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4997 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4998 ! Get parameter filenames and open the parameter files.
4999 call getenv_loc('BONDPAR',bondname)
5000 open (ibond,file=bondname,status='old',readonly,shared)
5001 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5002 open (ibond_nucl,file=bondname_nucl,status='old',readonly,shared)
5003 call getenv_loc('THETPAR',thetname)
5004 open (ithep,file=thetname,status='old',readonly,shared)
5005 call getenv_loc('ROTPAR',rotname)
5006 open (irotam,file=rotname,status='old',readonly,shared)
5007 call getenv_loc('TORPAR',torname)
5008 open (itorp,file=torname,status='old',readonly,shared)
5009 call getenv_loc('TORDPAR',tordname)
5010 open (itordp,file=tordname,status='old',readonly,shared)
5011 call getenv_loc('FOURIER',fouriername)
5012 open (ifourier,file=fouriername,status='old',readonly,shared)
5013 call getenv_loc('ELEPAR',elename)
5014 open (ielep,file=elename,status='old',readonly,shared)
5015 call getenv_loc('SIDEPAR',sidename)
5016 open (isidep,file=sidename,status='old',readonly,shared)
5018 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5019 open (ithep_nucl,file=thetname_nucl,status='old',readonly,shared)
5020 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5021 open (irotam_nucl,file=rotname_nucl,status='old',readonly,shared)
5022 call getenv_loc('TORPAR_NUCL',torname_nucl)
5023 open (itorp_nucl,file=torname_nucl,status='old',readonly,shared)
5024 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5025 open (itordp_nucl,file=tordname_nucl,status='old',readonly,shared)
5026 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5027 open (isidep_nucl,file=sidename_nucl,status='old',readonly,shared)
5030 #elif (defined CRAY) || (defined AIX)
5031 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
5033 ! print *,"Processor",myrank," opened file 1"
5034 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5035 ! print *,"Processor",myrank," opened file 9"
5036 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5037 ! Get parameter filenames and open the parameter files.
5038 call getenv_loc('BONDPAR',bondname)
5039 open (ibond,file=bondname,status='old',action='read')
5040 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5041 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
5043 ! print *,"Processor",myrank," opened file IBOND"
5044 call getenv_loc('THETPAR',thetname)
5045 open (ithep,file=thetname,status='old',action='read')
5046 ! print *,"Processor",myrank," opened file ITHEP"
5047 call getenv_loc('ROTPAR',rotname)
5048 open (irotam,file=rotname,status='old',action='read')
5049 ! print *,"Processor",myrank," opened file IROTAM"
5050 call getenv_loc('TORPAR',torname)
5051 open (itorp,file=torname,status='old',action='read')
5052 ! print *,"Processor",myrank," opened file ITORP"
5053 call getenv_loc('TORDPAR',tordname)
5054 open (itordp,file=tordname,status='old',action='read')
5055 ! print *,"Processor",myrank," opened file ITORDP"
5056 call getenv_loc('SCCORPAR',sccorname)
5057 open (isccor,file=sccorname,status='old',action='read')
5058 ! print *,"Processor",myrank," opened file ISCCOR"
5059 call getenv_loc('FOURIER',fouriername)
5060 open (ifourier,file=fouriername,status='old',action='read')
5061 ! print *,"Processor",myrank," opened file IFOURIER"
5062 call getenv_loc('ELEPAR',elename)
5063 open (ielep,file=elename,status='old',action='read')
5064 ! print *,"Processor",myrank," opened file IELEP"
5065 call getenv_loc('SIDEPAR',sidename)
5066 open (isidep,file=sidename,status='old',action='read')
5068 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5069 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
5070 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5071 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
5072 call getenv_loc('TORPAR_NUCL',torname_nucl)
5073 open (itorp_nucl,file=torname_nucl,status='old',action='read')
5074 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5075 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
5076 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5077 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
5079 call getenv_loc('LIPTRANPAR',liptranname)
5080 open (iliptranpar,file=liptranname,status='old',action='read')
5081 call getenv_loc('TUBEPAR',tubename)
5082 open (itube,file=tubename,status='old',action='read')
5083 call getenv_loc('IONPAR',ionname)
5084 open (iion,file=ionname,status='old',action='read')
5086 ! print *,"Processor",myrank," opened file ISIDEP"
5087 ! print *,"Processor",myrank," opened parameter files"
5089 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
5090 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5091 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5092 ! Get parameter filenames and open the parameter files.
5093 call getenv_loc('BONDPAR',bondname)
5094 open (ibond,file=bondname,status='old')
5095 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5096 open (ibond_nucl,file=bondname_nucl,status='old')
5098 call getenv_loc('THETPAR',thetname)
5099 open (ithep,file=thetname,status='old')
5100 call getenv_loc('ROTPAR',rotname)
5101 open (irotam,file=rotname,status='old')
5102 call getenv_loc('TORPAR',torname)
5103 open (itorp,file=torname,status='old')
5104 call getenv_loc('TORDPAR',tordname)
5105 open (itordp,file=tordname,status='old')
5106 call getenv_loc('SCCORPAR',sccorname)
5107 open (isccor,file=sccorname,status='old')
5108 call getenv_loc('FOURIER',fouriername)
5109 open (ifourier,file=fouriername,status='old')
5110 call getenv_loc('ELEPAR',elename)
5111 open (ielep,file=elename,status='old')
5112 call getenv_loc('SIDEPAR',sidename)
5113 open (isidep,file=sidename,status='old')
5115 open (ithep_nucl,file=thetname_nucl,status='old')
5116 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5117 open (irotam_nucl,file=rotname_nucl,status='old')
5118 call getenv_loc('TORPAR_NUCL',torname_nucl)
5119 open (itorp_nucl,file=torname_nucl,status='old')
5120 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5121 open (itordp_nucl,file=tordname_nucl,status='old')
5122 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5123 open (isidep_nucl,file=sidename_nucl,status='old')
5125 call getenv_loc('LIPTRANPAR',liptranname)
5126 open (iliptranpar,file=liptranname,status='old')
5127 call getenv_loc('TUBEPAR',tubename)
5128 open (itube,file=tubename,status='old')
5129 call getenv_loc('IONPAR',ionname)
5130 open (iion,file=ionname,status='old')
5132 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
5134 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5135 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5136 ! Get parameter filenames and open the parameter files.
5137 call getenv_loc('BONDPAR',bondname)
5138 open (ibond,file=bondname,status='old',action='read')
5139 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5140 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
5141 call getenv_loc('THETPAR',thetname)
5142 open (ithep,file=thetname,status='old',action='read')
5143 call getenv_loc('ROTPAR',rotname)
5144 open (irotam,file=rotname,status='old',action='read')
5145 call getenv_loc('TORPAR',torname)
5146 open (itorp,file=torname,status='old',action='read')
5147 call getenv_loc('TORDPAR',tordname)
5148 open (itordp,file=tordname,status='old',action='read')
5149 call getenv_loc('SCCORPAR',sccorname)
5150 open (isccor,file=sccorname,status='old',action='read')
5152 call getenv_loc('THETPARPDB',thetname_pdb)
5153 print *,"thetname_pdb ",thetname_pdb
5154 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
5155 print *,ithep_pdb," opened"
5157 call getenv_loc('FOURIER',fouriername)
5158 open (ifourier,file=fouriername,status='old',readonly)
5159 call getenv_loc('ELEPAR',elename)
5160 open (ielep,file=elename,status='old',readonly)
5161 call getenv_loc('SIDEPAR',sidename)
5162 open (isidep,file=sidename,status='old',readonly)
5164 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5165 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
5166 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5167 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
5168 call getenv_loc('TORPAR_NUCL',torname_nucl)
5169 open (itorp_nucl,file=torname_nucl,status='old',action='read')
5170 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5171 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
5172 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5173 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
5174 call getenv_loc('SIDEPAR_SCBASE',sidename_scbase)
5175 open (isidep_scbase,file=sidename_scbase,status='old',action='read')
5176 call getenv_loc('PEPPAR_PEPBASE',pepname_pepbase)
5177 open (isidep_pepbase,file=pepname_pepbase,status='old',action='read')
5178 call getenv_loc('SCPAR_PHOSPH',pepname_scpho)
5179 open (isidep_scpho,file=pepname_scpho,status='old',action='read')
5180 call getenv_loc('PEPPAR_PHOSPH',pepname_peppho)
5181 open (isidep_peppho,file=pepname_peppho,status='old',action='read')
5184 call getenv_loc('LIPTRANPAR',liptranname)
5185 open (iliptranpar,file=liptranname,status='old',action='read')
5186 call getenv_loc('TUBEPAR',tubename)
5187 open (itube,file=tubename,status='old',action='read')
5188 call getenv_loc('IONPAR',ionname)
5189 open (iion,file=ionname,status='old',action='read')
5192 call getenv_loc('ROTPARPDB',rotname_pdb)
5193 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
5196 call getenv_loc('SCPPAR_NUCL',scpname_nucl)
5197 #if defined(WINIFL) || defined(WINPGI)
5198 open (iscpp_nucl,file=scpname_nucl,status='old',readonly,shared)
5199 #elif (defined CRAY) || (defined AIX)
5200 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
5202 open (iscpp_nucl,file=scpname_nucl,status='old')
5204 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
5209 ! 8/9/01 In the newest version SCp interaction constants are read from a file
5210 ! Use -DOLDSCP to use hard-coded constants instead.
5212 call getenv_loc('SCPPAR',scpname)
5213 #if defined(WINIFL) || defined(WINPGI)
5214 open (iscpp,file=scpname,status='old',readonly,shared)
5215 #elif (defined CRAY) || (defined AIX)
5216 open (iscpp,file=scpname,status='old',action='read')
5218 open (iscpp,file=scpname,status='old')
5220 open (iscpp,file=scpname,status='old',action='read')
5223 call getenv_loc('PATTERN',patname)
5224 #if defined(WINIFL) || defined(WINPGI)
5225 open (icbase,file=patname,status='old',readonly,shared)
5226 #elif (defined CRAY) || (defined AIX)
5227 open (icbase,file=patname,status='old',action='read')
5229 open (icbase,file=patname,status='old')
5231 open (icbase,file=patname,status='old',action='read')
5234 ! Open output file only for CG processes
5235 ! print *,"Processor",myrank," fg_rank",fg_rank
5236 if (fg_rank.eq.0) then
5238 if (nodes.eq.1) then
5241 npos = dlog10(dfloat(nodes-1))+1
5243 if (npos.lt.3) npos=3
5244 write (liczba,'(i1)') npos
5245 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba)) &
5247 write (liczba,form) me
5248 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)// &
5249 liczba(:ilen(liczba))
5250 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
5252 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
5254 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)// &
5255 liczba(:ilen(liczba))//'.mol2'
5256 statname=prefix(:lenpre)//'_'//pot(:lenpot)// &
5257 liczba(:ilen(liczba))//'.stat'
5259 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) &
5260 //liczba(:ilen(liczba))//'.stat')
5261 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba)) &
5264 qname=prefix(:lenpre)//'_'//pot(:lenpot)// &
5265 liczba(:ilen(liczba))//'.const'
5270 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
5271 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
5272 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
5273 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
5274 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
5276 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)// &
5278 rest2name=prefix(:ilen(prefix))//'.rst'
5280 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
5283 #if defined(AIX) || defined(PGI)
5284 if (me.eq.king .or. .not. out1file) &
5285 open(iout,file=outname,status='unknown')
5287 if (fg_rank.gt.0) then
5288 write (liczba,'(i3.3)') myrank/nfgtasks
5289 write (ll,'(bz,i3.3)') fg_rank
5290 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
5295 open(igeom,file=intname,status='unknown',position='append')
5296 open(ipdb,file=pdbname,status='unknown')
5297 open(imol2,file=mol2name,status='unknown')
5298 open(istat,file=statname,status='unknown',position='append')
5300 !1out open(iout,file=outname,status='unknown')
5303 if (me.eq.king .or. .not.out1file) &
5304 open(iout,file=outname,status='unknown')
5306 if (fg_rank.gt.0) then
5307 write (liczba,'(i3.3)') myrank/nfgtasks
5308 write (ll,'(bz,i3.3)') fg_rank
5309 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
5314 open(igeom,file=intname,status='unknown',access='append')
5315 open(ipdb,file=pdbname,status='unknown')
5316 open(imol2,file=mol2name,status='unknown')
5317 open(istat,file=statname,status='unknown',access='append')
5319 !1out open(iout,file=outname,status='unknown')
5322 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
5323 csa_seed=prefix(:lenpre)//'.CSA.seed'
5324 csa_history=prefix(:lenpre)//'.CSA.history'
5325 csa_bank=prefix(:lenpre)//'.CSA.bank'
5326 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
5327 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
5328 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
5329 !!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
5330 csa_int=prefix(:lenpre)//'.int'
5331 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
5332 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
5333 csa_in=prefix(:lenpre)//'.CSA.in'
5334 ! print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
5337 write (iout,'(80(1h-))')
5338 write (iout,'(30x,a)') "FILE ASSIGNMENT"
5339 write (iout,'(80(1h-))')
5340 write (iout,*) "Input file : ",&
5341 pref_orig(:ilen(pref_orig))//'.inp'
5342 write (iout,*) "Output file : ",&
5343 outname(:ilen(outname))
5345 write (iout,*) "Sidechain potential file : ",&
5346 sidename(:ilen(sidename))
5348 write (iout,*) "SCp potential file : ",&
5349 scpname(:ilen(scpname))
5351 write (iout,*) "Electrostatic potential file : ",&
5352 elename(:ilen(elename))
5353 write (iout,*) "Cumulant coefficient file : ",&
5354 fouriername(:ilen(fouriername))
5355 write (iout,*) "Torsional parameter file : ",&
5356 torname(:ilen(torname))
5357 write (iout,*) "Double torsional parameter file : ",&
5358 tordname(:ilen(tordname))
5359 write (iout,*) "SCCOR parameter file : ",&
5360 sccorname(:ilen(sccorname))
5361 write (iout,*) "Bond & inertia constant file : ",&
5362 bondname(:ilen(bondname))
5363 write (iout,*) "Bending parameter file : ",&
5364 thetname(:ilen(thetname))
5365 write (iout,*) "Rotamer parameter file : ",&
5366 rotname(:ilen(rotname))
5369 write (iout,*) "Thetpdb parameter file : ",&
5370 thetname_pdb(:ilen(thetname_pdb))
5373 write (iout,*) "Threading database : ",&
5374 patname(:ilen(patname))
5376 write (iout,*)" DIRTMP : ",&
5378 write (iout,'(80(1h-))')
5381 end subroutine openunits
5382 !-----------------------------------------------------------------------------
5385 use geometry_data, only: nres,dc
5387 ! implicit real*8 (a-h,o-z)
5388 ! include 'DIMENSIONS'
5389 ! include 'COMMON.CHAIN'
5390 ! include 'COMMON.IOUNITS'
5391 ! include 'COMMON.MD'
5394 ! real(kind=8) :: var,ene
5396 open(irest2,file=rest2name,status='unknown')
5397 read(irest2,*) totT,EK,potE,totE,t_bath
5400 ! AL 4/17/17: Now reading d_t(0,:) too
5402 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
5405 ! AL 4/17/17: Now reading d_c(0,:) too
5407 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
5410 read (irest2,*) iset
5414 end subroutine readrst
5415 !-----------------------------------------------------------------------------
5416 subroutine read_fragments
5420 use control_data, only:out1file
5423 ! implicit real*8 (a-h,o-z)
5424 ! include 'DIMENSIONS'
5428 ! include 'COMMON.SETUP'
5429 ! include 'COMMON.CHAIN'
5430 ! include 'COMMON.IOUNITS'
5431 ! include 'COMMON.MD'
5432 ! include 'COMMON.CONTROL'
5435 ! real(kind=8) :: var,ene
5437 read(inp,*) nset,nfrag,npair,nfrag_back
5439 !el from module energy
5440 ! if(.not.allocated(mset)) allocate(mset(nset)) !(maxprocs/20)
5441 if(.not.allocated(wfrag_back)) then
5442 allocate(wfrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
5443 allocate(ifrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
5445 allocate(qinfrag(nfrag,nset),wfrag(nfrag,nset)) !(50,maxprocs/20)
5446 allocate(ifrag(2,nfrag,nset)) !(2,50,maxprocs/20)
5448 allocate(qinpair(npair,nset),wpair(npair,nset)) !(100,maxprocs/20)
5449 allocate(ipair(2,npair,nset)) !(2,100,maxprocs/20)
5452 if(me.eq.king.or..not.out1file) &
5453 write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,&
5454 " nfrag_back",nfrag_back
5456 read(inp,*) mset(iset)
5458 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),&
5460 if(me.eq.king.or..not.out1file) &
5461 write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),&
5462 ifrag(2,i,iset), qinfrag(i,iset)
5465 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),&
5467 if(me.eq.king.or..not.out1file) &
5468 write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),&
5469 ipair(2,i,iset), qinpair(i,iset)
5472 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
5473 wfrag_back(3,i,iset),&
5474 ifrag_back(1,i,iset),ifrag_back(2,i,iset)
5475 if(me.eq.king.or..not.out1file) &
5476 write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
5477 wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
5481 end subroutine read_fragments
5482 !-----------------------------------------------------------------------------
5484 !-----------------------------------------------------------------------------
5488 ! implicit real*8 (a-h,o-z)
5489 ! include 'DIMENSIONS'
5490 ! include 'COMMON.CSA'
5491 ! include 'COMMON.BANK'
5492 ! include 'COMMON.IOUNITS'
5494 ! integer :: ntf,ik,iw_pdb
5495 ! real(kind=8) :: var,ene
5497 open(icsa_in,file=csa_in,status="old",err=100)
5498 read(icsa_in,*) nconf
5499 read(icsa_in,*) jstart,jend
5500 read(icsa_in,*) nstmax
5501 read(icsa_in,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
5502 read(icsa_in,*) nran0,nran1,irr
5503 read(icsa_in,*) nseed
5504 read(icsa_in,*) ntotal,cut1,cut2
5505 read(icsa_in,*) estop
5506 read(icsa_in,*) icmax,irestart
5507 read(icsa_in,*) ntbankm,dele,difcut
5508 read(icsa_in,*) iref,rmscut,pnccut
5509 read(icsa_in,*) ndiff
5516 end subroutine csa_read
5517 !-----------------------------------------------------------------------------
5518 subroutine initial_write
5521 ! implicit real*8 (a-h,o-z)
5522 ! include 'DIMENSIONS'
5523 ! include 'COMMON.CSA'
5524 ! include 'COMMON.BANK'
5525 ! include 'COMMON.IOUNITS'
5527 ! integer :: ntf,ik,iw_pdb
5528 ! real(kind=8) :: var,ene
5530 open(icsa_seed,file=csa_seed,status="unknown")
5531 write(icsa_seed,*) "seed"
5533 #if defined(AIX) || defined(PGI)
5534 open(icsa_history,file=csa_history,status="unknown",&
5537 open(icsa_history,file=csa_history,status="unknown",&
5540 write(icsa_history,*) nconf
5541 write(icsa_history,*) jstart,jend
5542 write(icsa_history,*) nstmax
5543 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
5544 write(icsa_history,*) nran0,nran1,irr
5545 write(icsa_history,*) nseed
5546 write(icsa_history,*) ntotal,cut1,cut2
5547 write(icsa_history,*) estop
5548 write(icsa_history,*) icmax,irestart
5549 write(icsa_history,*) ntbankm,dele,difcut
5550 write(icsa_history,*) iref,rmscut,pnccut
5551 write(icsa_history,*) ndiff
5553 write(icsa_history,*)
5556 open(icsa_bank1,file=csa_bank1,status="unknown")
5557 write(icsa_bank1,*) 0
5561 end subroutine initial_write
5562 !-----------------------------------------------------------------------------
5563 subroutine restart_write
5566 ! implicit real*8 (a-h,o-z)
5567 ! include 'DIMENSIONS'
5568 ! include 'COMMON.IOUNITS'
5569 ! include 'COMMON.CSA'
5570 ! include 'COMMON.BANK'
5572 ! integer :: ntf,ik,iw_pdb
5573 ! real(kind=8) :: var,ene
5575 #if defined(AIX) || defined(PGI)
5576 open(icsa_history,file=csa_history,position="append")
5578 open(icsa_history,file=csa_history,access="append")
5580 write(icsa_history,*)
5581 write(icsa_history,*) "This is restart"
5582 write(icsa_history,*)
5583 write(icsa_history,*) nconf
5584 write(icsa_history,*) jstart,jend
5585 write(icsa_history,*) nstmax
5586 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
5587 write(icsa_history,*) nran0,nran1,irr
5588 write(icsa_history,*) nseed
5589 write(icsa_history,*) ntotal,cut1,cut2
5590 write(icsa_history,*) estop
5591 write(icsa_history,*) icmax,irestart
5592 write(icsa_history,*) ntbankm,dele,difcut
5593 write(icsa_history,*) iref,rmscut,pnccut
5594 write(icsa_history,*) ndiff
5595 write(icsa_history,*)
5596 write(icsa_history,*) "irestart is: ", irestart
5598 write(icsa_history,*)
5602 end subroutine restart_write
5603 !-----------------------------------------------------------------------------
5605 !-----------------------------------------------------------------------------
5606 subroutine write_pdb(npdb,titelloc,ee)
5608 ! implicit real*8 (a-h,o-z)
5609 ! include 'DIMENSIONS'
5610 ! include 'COMMON.IOUNITS'
5611 character(len=50) :: titelloc1
5612 character*(*) :: titelloc
5613 character(len=3) :: zahl
5614 character(len=5) :: liczba5
5616 integer :: npdb !,ilen
5620 ! real(kind=8) :: var,ene
5624 if (npdb.lt.1000) then
5625 call numstr(npdb,zahl)
5626 open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
5628 if (npdb.lt.10000) then
5629 write(liczba5,'(i1,i4)') 0,npdb
5631 write(liczba5,'(i5)') npdb
5633 open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
5635 call pdbout(ee,titelloc1,ipdb)
5638 end subroutine write_pdb
5639 !-----------------------------------------------------------------------------
5641 !-----------------------------------------------------------------------------
5642 subroutine write_thread_summary
5643 ! Thread the sequence through a database of known structures
5644 use control_data, only: refstr
5646 use energy_data, only: n_ene_comp
5648 ! implicit real*8 (a-h,o-z)
5649 ! include 'DIMENSIONS'
5651 use MPI_data !include 'COMMON.INFO'
5654 ! include 'COMMON.CONTROL'
5655 ! include 'COMMON.CHAIN'
5656 ! include 'COMMON.DBASE'
5657 ! include 'COMMON.INTERACT'
5658 ! include 'COMMON.VAR'
5659 ! include 'COMMON.THREAD'
5660 ! include 'COMMON.FFIELD'
5661 ! include 'COMMON.SBRIDGE'
5662 ! include 'COMMON.HEADER'
5663 ! include 'COMMON.NAMES'
5664 ! include 'COMMON.IOUNITS'
5665 ! include 'COMMON.TIME1'
5667 integer,dimension(maxthread) :: ip
5668 real(kind=8),dimension(0:n_ene) :: energia
5670 integer :: i,j,ii,jj,ipj,ik,kk,ist
5671 real(kind=8) :: enet,etot,rmsnat,rms,frac,frac_nn
5673 write (iout,'(30x,a/)') &
5674 ' *********** Summary threading statistics ************'
5675 write (iout,'(a)') 'Initial energies:'
5676 write (iout,'(a4,2x,a12,14a14,3a8)') &
5677 'No','seq',(ename(print_order(i)),i=1,nprint_ene),'ETOT',&
5678 'RMSnat','NatCONT','NNCONT','RMS'
5679 ! Energy sort patterns
5684 enet=ener(n_ene-1,ip(i))
5687 if (ener(n_ene-1,ip(j)).lt.enet) then
5689 enet=ener(n_ene-1,ip(j))
5701 ist=nres_base(2,ii)+ipatt(2,i)
5703 energia(i)=ener0(kk,i)
5705 etot=ener0(n_ene_comp+1,i)
5706 rmsnat=ener0(n_ene_comp+2,i)
5707 rms=ener0(n_ene_comp+3,i)
5708 frac=ener0(n_ene_comp+4,i)
5709 frac_nn=ener0(n_ene_comp+5,i)
5712 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5713 i,str_nam(ii),ist+1,&
5714 (energia(print_order(kk)),kk=1,nprint_ene),&
5715 etot,rmsnat,frac,frac_nn,rms
5717 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3)') &
5718 i,str_nam(ii),ist+1,&
5719 (energia(print_order(kk)),kk=1,nprint_ene),etot
5722 write (iout,'(//a)') 'Final energies:'
5723 write (iout,'(a4,2x,a12,17a14,3a8)') &
5724 'No','seq',(ename(print_order(kk)),kk=1,nprint_ene),'ETOT',&
5725 'RMSnat','NatCONT','NNCONT','RMS'
5729 ist=nres_base(2,ii)+ipatt(2,i)
5731 energia(kk)=ener(kk,ik)
5733 etot=ener(n_ene_comp+1,i)
5734 rmsnat=ener(n_ene_comp+2,i)
5735 rms=ener(n_ene_comp+3,i)
5736 frac=ener(n_ene_comp+4,i)
5737 frac_nn=ener(n_ene_comp+5,i)
5738 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5739 i,str_nam(ii),ist+1,&
5740 (energia(print_order(kk)),kk=1,nprint_ene),&
5741 etot,rmsnat,frac,frac_nn,rms
5743 write (iout,'(/a/)') 'IEXAM array:'
5744 write (iout,'(i5)') nexcl
5746 write (iout,'(2i5)') iexam(1,i),iexam(2,i)
5748 write (iout,'(/a,1pe14.4/a,1pe14.4/)') &
5749 'Max. time for threading step ',max_time_for_thread,&
5750 'Average time for threading step: ',ave_time_for_thread
5752 end subroutine write_thread_summary
5753 !-----------------------------------------------------------------------------
5754 subroutine write_stat_thread(ithread,ipattern,ist)
5756 use energy_data, only: n_ene_comp
5758 ! implicit real*8 (a-h,o-z)
5759 ! include "DIMENSIONS"
5760 ! include "COMMON.CONTROL"
5761 ! include "COMMON.IOUNITS"
5762 ! include "COMMON.THREAD"
5763 ! include "COMMON.FFIELD"
5764 ! include "COMMON.DBASE"
5765 ! include "COMMON.NAMES"
5766 real(kind=8),dimension(0:n_ene) :: energia
5768 integer :: ithread,ipattern,ist,i
5769 real(kind=8) :: etot,rmsnat,rms,frac,frac_nn
5771 #if defined(AIX) || defined(PGI)
5772 open(istat,file=statname,position='append')
5774 open(istat,file=statname,access='append')
5777 energia(i)=ener(i,ithread)
5779 etot=ener(n_ene_comp+1,ithread)
5780 rmsnat=ener(n_ene_comp+2,ithread)
5781 rms=ener(n_ene_comp+3,ithread)
5782 frac=ener(n_ene_comp+4,ithread)
5783 frac_nn=ener(n_ene_comp+5,ithread)
5784 write (istat,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5785 ithread,str_nam(ipattern),ist+1,&
5786 (energia(print_order(i)),i=1,nprint_ene),&
5787 etot,rmsnat,frac,frac_nn,rms
5790 end subroutine write_stat_thread
5791 !-----------------------------------------------------------------------------
5793 !-----------------------------------------------------------------------------
5794 end module io_config